Estimation and uncertainty quantification for extreme quantile regions

TL;DR

This paper introduces Bayesian methods for estimating extreme quantile regions and quantifies the uncertainty of these estimates, outperforming existing approaches and demonstrated through simulations and pollution data analysis.

Contribution

The paper develops novel Bayesian schemes for univariate and bivariate extreme quantile region estimation with uncertainty quantification, improving upon existing methods.

Findings

Bayesian methods outperform existing approaches in simulations.

Proposed methods effectively quantify uncertainty in extreme quantile estimates.

Application to pollution data demonstrates practical utility.

Abstract

Estimation of extreme quantile regions, spaces in which future extreme events can occur with a given low probability, even beyond the range of the observed data, is an important task in the analysis of extremes. Existing methods to estimate such regions are available, but do not provide any measures of estimation uncertainty. We develop univariate and bivariate schemes for estimating extreme quantile regions under the Bayesian paradigm that outperforms existing approaches and provides natural measures of quantile region estimate uncertainty. We examine the method's performance in controlled simulation studies. We illustrate the applicability of the proposed method by analysing high bivariate quantiles for pairs of pollutants, conditionally on different temperature gradations, recorded in Milan, Italy.

Figures

Equations

$$\mathcal{Q}=\{\mathbf{x}\in {\mathbb{R}}_{+}^2:f(\mathbf{x})⩽\alpha \},$$

$$F^n(a_{n}x+b_n)\stackrel{n\to \infty }{⟶}\exp \left(-{\left(1+\gamma x\right)}_{+}^{-1/\gamma }\right)=:G(x;\gamma ),$$

$$t\{1-F(a(t)x+U(t))\}\stackrel{t\to \infty }{⟶}{\left(1+\gamma x\right)}^{-1/\gamma },$$

$$F(y)\approx 1-\frac{k}{n}{\left(1+\gamma \frac{y-\mu }{\sigma }\right)}_{+}^{-1/\gamma },n\to \infty ,$$

$$\frac{U(tx)-U(t)}{a(t)}\stackrel{t\to \infty }{⟶}\frac{x^{\gamma }-1}{\gamma },$$

$$Q(p)\approx \mu +\sigma \frac{{\left(\frac{k}{np}\right)}^{\gamma }-1}{\gamma }\text{mbox}asn\to \infty .$$

$$F_{a_{n}s+b_n}(a_{n}x+b_n)\stackrel{n\to \infty }{⟶}1-{\left(1+\gamma \frac{x-s}{\tilde{\sigma }}\right)}_{+}^{-1/\gamma }=:H((x-s)/\tilde{\sigma };\gamma ),$$

$$F(y)$$

$$F^n({\mathbf{a}}_n\mathbf{x}+{\mathbf{b}}_n)\stackrel{n\to \infty }{⟶}G(\mathbf{x}),$$

$$t\left(1-F(t{U}_1(x_1),t{U}_2(x_2))\right)\stackrel{t\to \infty }{⟶}-\log G_{\ast }(\mathbf{x}).$$

$$-\log G_{\ast }(\mathbf{x})=\nu (\{\mathbf{v}\in {\mathbb{R}}_{+}^2:v_1>x_1\text{ or }{v}_2>x_2\})=2{\int }_0^1\left(\frac{w}{x_1}\vee \frac{1-w}{x_2}\right)H(\text{d}w),$$

$$\nu (\{\mathbf{v}\in {\mathbb{R}}_{+}^2:v_1>x_1\text{mbox}or{v}_2>x_2\})={\iint }_{\{v_1>x_1\text{mbox}or{v}_2>x_2\}}g(\mathbf{v})\text{d}\mathbf{v}.$$

$$t{U}_1(t)U_2(t)f\left(t{U}_1(x_1),t{U}_2(x_2)\right)\stackrel{t\to \infty }{⟶}({\gamma }_1{\gamma }_2{)}^{-1}{x}_1^{1-{\gamma }_1}{x}_2^{1-{\gamma }_2}g(\mathbf{x})=:q(\mathbf{x}).$$

$$q_{\ast }(w)=q(w,1-w{)}^{-\frac{1}{1+{\gamma }_1+{\gamma }_2}},q(w,1-w)=\frac{2w^{1-{\gamma }_1}(1-w{)}^{1-{\gamma }_2}h(w)}{{\gamma }_1{\gamma }_2}.$$

$$\nu (\mathcal{S})=2{\int }_{[0,1]}{q}_{\ast }^{-1}(w)h(w)\text{d}w.$$

$${\tilde{\mathcal{Q}}}_n\approx \left\{\left({\mu }_1+{\sigma }_1\frac{{\left(\frac{k\nu (\mathcal{S})x_1}{np}\right)}^{{\gamma }_1}-1}{{\gamma }_1},{\mu }_2+{\sigma }_2\frac{{\left(\frac{k\nu (\mathcal{S})x_2}{np}\right)}^{{\gamma }_2}-1}{{\gamma }_2}\right):(x_1,x_2)\in \mathcal{S}\right\}.$$

$$\mathcal{L}({\mathbf{y}}_{1:n};\mathbf{\theta })=\prod _{i=1}^n\mathcal{L}(y_i;\mathbf{\theta }),$$

$$\mathcal{L}(y_i;\mathbf{\theta })\propto \{\begin{array}{ll}{G}^{k/n}(t;\mathbf{\theta }),& \text{mbox}if{y}_i⩽t,\\ \frac{\partial }{\partial y}{G}^{k/n}(y;\mathbf{\theta }){|}_{y=y_i},& \text{mbox}if{y}_i>t,\end{array}$$

$$\frac{\partial }{\partial y}{G}^{k/n}(y;\mathbf{\theta }){|}_{y=y_i}=G^{k/n}((y_i-\mu )/\sigma ;\gamma ){\left(1+\frac{\gamma (y_i-\mu )}{\sigma }\right)}^{-1/\gamma -1}\frac{1}{\sigma }\frac{k}{n}.$$

$${\pi }^{(j)}=\min \left(\frac{\mathcal{L}({\mathbf{y}}_{1:n};{\mathbf{\theta }}^{\prime })\Pi ({\mathbf{\theta }}^{\prime })}{\mathcal{L}({\mathbf{y}}_{1:n};{\mathbf{\theta }}^{(j)})\Pi ({\mathbf{\theta }}^{(j)})},1\right),$$

$$\Sigma^{(j+1)}=\left\{\begin{array}[]{ll}(1+[\tau^{(j)}]^{2}/j)\mathrm{I}_{3},&j\leqslant 100\\ \frac{1}{j-1}\sum_{k=1}^{j}({\bm{\theta}}^{(k)}-\bar{{\bm{\theta}}}^{(j)})({\bm{\theta}}^{(k)}-\bar{{\bm{\theta}}}^{(j)})^{\top}+([\tau^{(j)}]^{2}/j)\mathrm{I}_{3},&j>100,\end{array}\right.$$

$$\log {\tau }^{(j+1)}=\log {\tau }^{(j)}+c({\pi }^{(j)}-{\pi }^{\ast }),$$

$$z_i:=z_i(\cdot )\equiv z_i(\cdot ;\mathbf{\theta })=\frac{k}{n}{\left(1+{\gamma }_i\frac{\cdot -{\mu }_i}{{\sigma }_i}\right)}_{+}^{-1/{\gamma }_i},i=1,2,$$

$$F(\mathbf{y})\approx G(\mathbf{y};{\mathbf{\vartheta }}^{\ast }),$$

$$F(\mathbf{y})\approx \exp \left(-L(\mathbf{z};\mathbf{\vartheta })\right)=:\tilde{G}(\mathbf{z};\mathbf{\vartheta }),\mathbf{y}⩾\mathbf{t},$$

$$A_{\kappa }(v;{\mathbf{\beta }}_{\kappa })=\frac{1}{\kappa +1}\sum _{j=0}^{\kappa }{\beta }_j\text{Be}(v;j+1,\kappa -j+1),$$

$$h_{\kappa -1}(w;{\mathbf{\eta }}_{\kappa })=\sum _{j=0}^{\kappa -2}({\eta }_{j+1}-{\eta }_j)\text{Be}(w;j+1,\kappa -j-1),$$

$$\mathcal{L}(\mathbf{\vartheta })=\prod _{i=1}^n\mathcal{L}({\mathbf{y}}_i;\mathbf{\vartheta }),$$

$$\mathcal{L}(\mathbf{y};\mathbf{\vartheta })\propto \{\begin{array}{ll}\tilde{G}(\mathbf{z};\mathbf{\vartheta }){|}_{\mathbf{z}=\mathbf{z}(\mathbf{t})},& \text{mbox}if\mathbf{y}⩽\mathbf{t}\\ -\frac{\partial z_1(y_1)}{\partial y_1}\tilde{G}(z_1,z_2;\mathbf{\vartheta })L^{(z_1)}(z_1,z_2;\mathbf{\vartheta }){|}_{z_1=z_1(y_1),z_2=z_2(t_2)},& \text{mbox}if{y}_1>t_1,y_2⩽t_2\\ -\frac{\partial z_2(y_2)}{\partial y_2}\tilde{G}(z_1,z_2;\mathbf{\vartheta })L^{(z_2)}(z_1,z_2;\mathbf{\vartheta }){|}_{z_1=z_1(t_1),z_2=z_2(y_2)},& \text{mbox}if{y}_1⩽t_1,y_2>t_2\\ \mathcal{J}\cdot \tilde{G}(\mathbf{y};\mathbf{\vartheta })(L^{(z_1)}(\mathbf{z};\mathbf{\vartheta })L^{(z_2)}(\mathbf{z};\mathbf{\vartheta })-L^{(z_1,z_2)}(\mathbf{z};\mathbf{\vartheta }){|}_{\mathbf{z}=\mathbf{z}(\mathbf{y})},& \text{mbox}if\mathbf{y}>\mathbf{t},\end{array}$$

$$\mathcal{J}=\frac{\partial z_1(y_1)}{\partial y_1}\cdot \frac{\partial z_2(y_2)}{\partial y_2},$$

Full text

Keywords: Air pollution; Bayesian nonparametrics; Bernstein polynomials; Extremal dependence; Extreme quantile regions; Max-stable distributions.

1 Introduction

Estimating quantiles and how they can change depending on influential predictors (i.e. quantile regression) is a recurring problem in many applied fields including medicine, survival analysis, economics, finance and environmental science (e.g. Yu et al., 2003). Assessing high quantiles, associated with fixed, low occurrence probabilities that lie beyond the range of $n$ existing data observations, is of crucial importance in risk management. The solution to such a problem is not obvious when none of the observed data points exceed this event. It is likely that the exceedance probability, $p$, is smaller than $1/n$. This is an extreme value problem.

More precisely, consider a random variable $X$ with distribution function $F$ defined on $\mathbb{R}_{+}:=(0,\infty)$. For $0<p<1$, let $Q(p):=F^{\leftarrow}(1-p)$ be the $(1-p)$-th quantile of $F$, where $F^{\leftarrow}$ is the left-continuous inverse function of $F$, i.e. $F^{\leftarrow}(x):=\inf\{y:F(y)\geqslant x\}$. Let $X_{1},\ldots,X_{n}$ be independent random variables with distribution $F$. Interest is in estimating the quantile $Q(p)$ when the exceedance probability $p$ is very small. We refer to $Q$ as an extreme quantile. The extreme-value approach (e.g. Beirlant et al., 2004; de Haan and Ferreira, 2006) assumes that as the sample size $n$ grows to infinity a suitable asymptotic probabilistic model can be used to approximate the desired quantile. In the asymptotic setting the exceedance probability depends on $n$, i.e. $p=p_{n}$ with $p\to 0$ as $n\to\infty$. Hence, for a sufficiently large number of observations, extreme-value theory can be used to compute an approximation, $Q_{n}$, of the extreme quantile.

A more challenging problem is estimating, for some fixed probability, an extreme bivariate region – a subset of two-dimensional Euclidean space – in which a future event would fall when none of the two-dimensional observations fall in such a region, and are likely to lie far from it (e.g. de Haan and Ferreira, 2006, Ch. 6). In some applications there may not be any specific shape for the critical region, whereas in others it may be well defined. One illustration of the former is air pollution monitoring, where the critical region is a set of combinations of two pollutants’ concentrations where at least one of them is at a high level. Here, the shape of the critical region depends on the type of pollutants and the intensity of their dependence.

Einmahl et al. (2013) proposed a simple and practically useful method for defining critical regions which are generated through the level sets of a probability density function $f$, under suitable conditions (see Section 2.2). Precisely, a critical region is the level set of $f$ given by

[TABLE]

with $\alpha>0$, such that $\mathbb{P}({\mathcal{Q}})=p$ for some very small $p\in(0,1)$. We refer to ${\mathcal{Q}}$ as an extreme quantile region. In particular ${\mathcal{Q}}^{\complement}=\{{\bm{x}}\in\mathbb{R}^{2}_{+}:\,f({\bm{x}})>\alpha\}$ is the set with smallest area such that $\mathbb{P}({\mathcal{Q}}^{\complement})=1-p$. See also Cooley et al. (2019) for an alternative approach. Let ${\bm{X}}_{1},\ldots,{\bm{X}}_{n}$ be independent two-dimensional random vectors with distribution function $F$. Similar to the univariate case, the extreme-value approach suggests adopting an appropriate asymptotic probabilistic model as $n$ approaches infinity. From this, a set ${\mathcal{Q}}_{n}$, that is not fixed but depends on $n$, may be derived such that $\mathbb{P}({\mathcal{Q}}_{n})=p_{n}$ with $p\equiv p_{n}$ and $p\to 0$ as $n\to\infty$, which provides a sufficiently close approximation to the extreme quantile region. In particular, in both univariate and bivariate cases we assume the mild condition that $np\to c\in[0,\infty)$ as $n\to\infty$.

Beyond a definition and a method of estimation of such quantiles, a critical aspect in practice is to provide some quantification of the uncertainty around the estimates. Several estimators of extreme quantiles already exist in the univariate setting (e.g. Beirlant et al., 2004; de Haan and Ferreira, 2006; Wang and Li, 2015). However, quantifying the uncertainty of these estimators can be difficult as their asymptotic variance depends on the parameters of second-order conditions which can be problematic to estimate in some settings. Estimators of extreme quantile regions also exist (Cai et al., 2011; Einmahl et al., 2013) but to the best of our knowledge measures of their uncertainty are not available. A readily applicable uncertainty quantification for extreme quantiles is the methodological contribution of this article.

In this paper we develop a Bayesian approach for inferring both univariate extreme quantiles as well as extreme quantile regions. In the univariate case, we define a parametric Bayesian method for the extreme quantiles by exploiting the well known and widely used censored-likelihood approach (e.g. Prescott and Walden, 1983; Smith, 1994; Ledford and Tawn, 1996; Huser et al., 2016; Bienvenüe and Robert, 2017), based on the likelihood function of the univariate Generalized Extreme Value (GEV) family of distributions (e.g. Sisson and Coles, 2003). Inference is performed using an adaptive random-walk Metropolis-Hastings algorithm (Garthwaite et al., 2016). In the bivariate case we combine the univariate approach for the estimation of the GEV marginal parameters with the non-parametric Bayesian approach for the estimation of the extremal dependence proposed by Marcon et al. (2016). As a result we obtain a semi-parametric Bayesian inferential method based on the censored-likelihood corresponding to a suitable two-dimensional extreme-value distribution. Through such a censored-likelihood it is possible to simultaneously estimate the dependence structure of the extreme-value distribution together with the parameters of its margins, which are in turn members of the GEV class. The components of such a pseudo-posterior distribution are combined to produce a pseudo-posterior distribution for the extreme quantile regions. Accordingly this approach allows for the direct estimation of extreme quantile intervals and quantile regions, as well as clear measures of their uncertainty.

There are many settings in environmental science where the analysis of univariate extreme quantiles is of particular concern. These include evaluating tropical cyclone wind speed conditional on certain climate variables (Jagger and Elsner, 2009), precipitation conditional on global climate model projections (Friederichs and Thorarinsdottir, 2012), both precipitation and wind power given outputs from numerical weather prediction models (Bremnes, 2004a, b), and ozone and particulate matter concentrations conditional on meteorological variables (Porter et al., 2015). However it is clear that some environmental variables depend on others: e.g. wind speed is dependent on wind gust or air pressure (e.g. Marcon et al., 2017), and ozone is dependent on particular matter (e.g. Falk et al., 2019). Similarly, when designing offshore structures, engineers need estimates of extreme quantiles of oceanographic variables. Due the strong dependence between the involved variables e.g. wave heights and wind speeds, bivariate extreme value tools (in the simplest case) are sometimes necessary (Naess and Karpa, 2015). As a consequence, the joint study of two (or more) dependent variables will produce a more comprehensive and accurate analysis because it is based on more available information.

The goal of this article is to provide an inferential framework (i.e. a method and software to implement it) that can provide estimates of extreme bivariate quantile regions, which can then be used to inform part of a larger study. For example, reliability engineers designing offshore platforms can use this framework as both a visual and quantitative tool to produce extreme quantile curves (such as those in Figure 6) that can be informative for reliability design. Here, we particularly focus on the problem of air pollution.

According to the World Health Organisations (WHO) and the Global Ambient Air Quality Database (2018 update, https://www.who.int/airpollution/en/), air pollution kills an estimated 7 million people worldwide each year, $4.2$ million of which as a result of exposure to ambient (outdoor) air pollution. Furthermore, $25\%$ of all heart diseases deaths are attributable to air pollution. Based on WHO guidelines (World Health Organization, 2006) the European policy on air-quality standards (Guerreiro et al., 2016) regulates emissions of the pollutants particulate matter (PM10), ozone (O3), nitrogen oxide (NO), nitrogen dioxide (NO2) and sulphur dioxide (SO2), aiming to reduce the negative impact that they have on human health, the environment and climate. For example, short-term pollutant concentrations (obtained through the empirical quantiles of pollutants time series) that should not be exceeded in order to protect against to air pollution peaks are: a daily average concentration of $50\,\mu$g/m3 ($90.4$ percentile) for PM10, a daily maximum concentration of $120\,\mu$g/m3 ($93.2$ percentile) for O3, an hourly average concentration of $200\,\mu$g/m3 ($99.8$ percentile) for NO2 and a daily average concentration of $125\,\mu$g/m3 ($99.2$ percentile) for SO2. We refer to such concentrations as the pollutant limit thresholds. In contrast, long-term pollutant thresholds are set based on yearly mean concentrations, rather than using quantiles. See Guerreiro et al. (2016) for details. The above thresholds are determined for each individual pollutant, although it is well known that some pollutants are dependent on each other (e.g. Dahlhaus, 2000; Clapp and Jenkin, 2001; Heffernan and Tawn, 2004; World Health Organization, 2006).

With our methodology (and software) we are able to estimate extreme quantiles while taking into account the dependence between pairs of pollutants, and thereby be informed about potentially dangerous pollutant concentration combinations. This will be of particular interest to air pollution emission regulators and public health analysts. Estimating extreme pollutant concentrations conditional on e.g. meteorological variables can also provide information on whether some predictors are important for understanding the evolution of air pollution at high concentration levels. Here, we investigate this explicitly through an analysis of the behaviour air pollution data recorded over the last $\sim 20$ years in Milan, Italy.

The remainder of this article is organised as follows. In Section 2 we briefly review the extreme-value approach for approximating the extreme quantiles in both the univariate (Section 2.1) and bivariate (Section 2.2) cases. In Section 3 we introduce our Bayesian semi-parametric approach for inferring extreme quantiles and extreme quantile regions. Section 4 provides an extensive simulation study examining the performance of the proposed method. The article concludes with an analysis of the extremes of air pollution recorded in Milan (Section 5) followed by a Discussion.

2 Estimating extreme quantiles

2.1 The univariate case

Let $F$ be a distribution on $\mathbb{R}_{+}$, and assume that $F$ is in the domain of attraction of the GEV family of distributions, $F\in{\mathcal{D}}(G)$. Then, for $n=1,2,\ldots$ there are norming constants $a_{n}>0$ and $b_{n}$ such that

[TABLE]

where $\gamma>0$ is the extreme-value index that describes the heaviness of the tail of the distribution $G$ and $(a)_{+}=\max(a,0)$ (see e.g. de Haan and Ferreira, 2006, Ch. 1).

There are several extreme-value based methods for modelling extreme quantiles; see e.g. de Haan and Ferreira (2006, Ch. 1, 4) and Wang and Li (2015) for a compendium. Here, we briefly review one which is useful for our purposes. By de Haan and Ferreira (2006, Theorem 1.1.6) the result in (2.1) is also equivalent to the result

[TABLE]

where $U(t)=F^{\leftarrow}(1-1/t)$ for $t>1$ and $a(\cdot)$ is a suitable function (de Haan and Ferreira, 2006, Theorem 1.1.6). Let $k\equiv k_{n}$ and assume that $k\to\infty$ and $k/n\to 0$ as $n\to\infty$. Set $t=n/k$ and $y=a(n/k)x+U(n/k)$. Then, from (2.2) we obtain

[TABLE]

where $\mu,\sigma>0$ are location and scale parameters parameters. By few manipulations and de Haan and Ferreira (2006, Theorem 1.1.8) we have that $a(t)\approx\sigma$ and $U(t)\approx\mu$ as $t\to\infty$. Result (2.1) is also equivalent to

[TABLE]

see de Haan and Ferreira (Theorem 1.1.6 2006, for details). Since $Q(p)\equiv U(1/p)$, from (2.4) we obtain, using arguments in de Haan and Ferreira (2006, Ch. 3.1),

[TABLE]

A suitable adjustment of the GEV distribution allows us to also derive the approximate quantile function in (2.5). Specifically, for some threshold $s>0$ and for $x\geqslant s$, let $F_{s}(x):=\mathbb{P}(X\leqslant x|X>s)$. Using de Haan and Ferreira (2006, Theorem 1.1.6), for $a_{n}=a(n)$ and $b_{n}=U(n)$, for all $x\geqslant s$ we can write

[TABLE]

where $\tilde{\sigma}=1+\gamma s$ and $H(\cdot;\gamma)$ is the Generalized Pareto (GP) family of distributions. Set $y=a_{n}x+b_{n}$ and $t:=a_{n}s+b_{n}$, then for large $t$ and any $y\geqslant t$ we obtain

[TABLE]

where $\bar{\sigma}=\sigma+\gamma(t-\mu)$. Writing again $F(y)=1-p$ and $t=n/k$ in (2.7), and by noting that $-\log(1-p)\approx p$ for $p\to 0$, then the quantile function in (2.5) can be obtained by inverting the expression in (2.7) with respect to $p$.

2.2 The bivariate case

Let ${\bm{X}}=(X_{1},X_{2})$ be a random vector with joint distribution function $F$ on $\mathbb{R}^{2}_{+}$ with margins $F_{j}$, $j=1,2$, and probability density function $f$. Assuming that $F\in\mathcal{D}(G)$, then there are sequences of norming constants ${\bm{a}}_{n}>\bm{0}$ and ${\bm{b}}_{n}$ such that

[TABLE]

where $G$ is the so-called bivariate max-stable distribution (de Haan and Ferreira, 2006, Ch 6), whose margins $G(x_{j};\gamma_{j})$ are members of the GEV family (2.1) with tail indices $\gamma_{j}>0$, $j=1,2$. In the following we denote $G_{*}({\bm{x}}):=G(x_{1}^{\gamma_{1}},x_{2}^{\gamma_{2}}),\quad{\bm{x}}\in\mathbb{R}_{+}^{2},$ as an extreme-value distribution with unit-Fréchet margins: $G_{*}(x_{1},\infty)=\exp(-1/x_{1})$ and $G_{*}(\infty,x_{2})=\exp(-1/x_{2})$ for every $x_{1},x_{2}>0$.

The convergence result in (2.8) implies convergence at both marginal and dependence levels. For marginal convergence we have that (2.1), (2.2) and (2.3) hold for $F_{j}$, with $a_{n,j}$, $b_{n,j}$, $a_{j}(t)$, $U_{j}(t)$ and $y_{j}=a_{j}(n/k)x_{j}+U_{j}(n/k)$ with $j=1,2$. Furthermore, $a_{j}(t)\approx\sigma_{j}>0$ and $U_{j}(t)\approx\mu_{j}>0$ as $t\to\infty$. For convergence of the dependence structure, for every $(x,y)\in(0,\infty]^{2}\setminus\{(\infty,\infty)\}$, from de Haan and Ferreira (2006, Ch 6.1.2) we have that

[TABLE]

The result in (2.9) implies the existence of a measure $\nu$, named the exponent measure (see de Haan and Ferreira, 2006, Ch 6 for details) such that for every ${\bm{x}}\in\mathbb{R}_{+}^{2}$,

[TABLE]

where $H$ is a probability measure on $[0,1]$ satisfying the mean constraint $\int_{0}^{1}wH(\textrm{d}w)=\int_{0}^{1}(1-w)H(\textrm{d}w)=1/2.$ In the following we denote both the probability measure and its distribution function by $H$, with the difference being determined by the context.

The extreme quantile regions method introduced in Einmahl et al. (2013) requires also some further assumptions at the density level. The density $f$ is assumed to be decreasing in each coordinate, outside of $(0,M]^{2}$ for some $M>0$, and bounded away from zero on $(0,M]^{2}$. There exist a nonnegative Lebesgue integrable function $g$ such that for every ${\bm{x}}\in\mathbb{R}_{+}^{2}$,

[TABLE]

and on $\mathbb{R}^{2}_{+}$,

[TABLE]

We refer to $g$ and $q$ as the density of the exponent measure the basic density function, respectively. By the above assumption on $g$ it follows that $H$ has a continuous density $h:=\partial H/\partial w$ on $(0,1)$ with no atoms at [math] and $1$, i.e. $H(\{0\})=H(\{1\})=0$. By the homogeneity of $g$ we have that $h(w)=2^{-1}g(w,1-w)$, where $r=x_{1}+x_{2}$ and $w=x_{1}/r$. We refer to $H$ and $h$ as the angular measure and density, respectively.

The extreme-value approach for modelling extreme quantile regions ${\mathcal{Q}}$ in (1.1), works with the set ${\mathcal{Q}}_{n}=\left\{{\bm{x}}\in\mathbb{R}^{2}\colon f({\bm{x}})\leqslant\alpha\right\}$, where $\alpha=\alpha_{n}$ is not fixed but depends on the sample size $n$ such that $\mathbb{P}({\mathcal{Q}}_{n})=p$ with $p=p_{n}\to 0$ as $n\to\infty$. In particular, Einmahl et al. (2013) suggest focusing on the fixed set ${\mathcal{S}}=\{{\bm{x}}\in\mathbb{R}^{2}_{+}:q({\bm{x}})\leqslant 1\}=\left\{{\bm{x}}\in\mathbb{R}^{2}_{+}\colon r\geqslant q_{*}^{-1}(w),w\in[0,1]\right\}$ that we call the basic set, where

[TABLE]

and where $q(w,1-w)$ is obtained from the relation $q(rw,r(1-w))=q(w,1-w)/r^{(1+\gamma_{1}+\gamma_{2})}$. Abusing terminology refer to $q_{*}$ as the angular basic density function. According to the exponent measure, the size of ${\mathcal{S}}$ is

[TABLE]

The idea is to then inflate the basic set ${\mathcal{S}}$ into an extreme set, ${\widetilde{\mathcal{Q}}_{n}}$, depending on $n$, so that $\mathbb{P}({\widetilde{\mathcal{Q}}_{n}})\approx p$ for $n\to\infty$ and such that ${\widetilde{\mathcal{Q}}_{n}}$ is a good approximation of ${\mathcal{Q}}_{n}$, i.e. so that $\mathbb{P}({\mathcal{Q}}_{n}\triangle{\widetilde{\mathcal{Q}}_{n}})/p\to 0$ as $n\to\infty$, where $B\triangle D=B\backslash D\,\cup\,D\backslash B$, for two nonempty sets $B$ and $D$. Here we consider a slightly different definition of the set ${\widetilde{\mathcal{Q}}_{n}}$ to that given in Einmahl et al. (2013). Specifically, by exploiting the fact that $U_{j}(n/k)\approx\mu_{j}$ and $a_{j}(n/k)\approx\sigma_{j}$, $j=1,2$ as $n\to\infty$, then for large $n$, we have that ${\widetilde{\mathcal{Q}}_{n}}$ can be approximated as

[TABLE]

The approximation in (2.13) is consistent with the formula in (2.5) used to approximate the univariate extreme quantiles. In Section 3.2 we show how the bivariate max-stable distribution in (2.8) can be used to estimate the extreme quantile region in (2.13).

3 Inference

3.1 The univariate case

We describe an approximate Bayesian framework for estimating the extreme quantile in (2.5) for small $p$ (where the meaning of “small $p$” is given in Section B.1). In particular, we explore the Bayesian paradigm using a censored likelihood (Sisson and Coles, 2003), based on (2.7). Specifically, let $X_{1},\ldots,X_{n}$ be independent and identically distributed random variables with distribution function $F$ on $\mathbb{R}_{+}$, where $F\in{\mathcal{D}}(G)$. First we define a high threshold $t$. Let $X_{1,n}\leqslant X_{2,n}\leqslant\cdots\leqslant X_{n,n}$ be the $n$ order statistics and $F_{n}$ be the empirical distribution function. The threshold may then be defined as $T=X_{n-k,n}$ for large $k$ such that $1-F_{n}(X_{n-k,n})=k/n$ is close to zero, for instance $k/n=0.10,0.05,0.01$. Next, let ${\bm{y}}_{1:n}=(y_{1},\ldots,y_{n})$ be a realisation of $(X_{1},\ldots,X_{n})$ and $t$ be the corresponding threshold. Then on the basis of the approximation in (2.7) we define the censored likelihood function

[TABLE]

where each contribution to the likelihood depends on the domain where an observation $y_{i}$ falls. That is,

[TABLE]

where $G^{k/n}(y;{\bm{\theta}})\equiv G^{k/n}((y-\mu)/\sigma;\gamma)$ with ${\bm{\theta}}=(\mu,\sigma,\gamma)^{\top}$, $G^{k/n}((y-\mu)/\sigma;\gamma)$ is given in (2.7) and

[TABLE]

Assuming a prior $\Pi({\bm{\theta}})$ for ${\bm{\theta}}\in{\bm{\Theta}}\subseteq\mathbb{R}_{+}^{3}$, we draw samples from the resulting posterior distribution using an adaptive Metropolis-Hastings algorithm. Specifically, we directly apply the adaptive (Gaussian) random-walk Metropolis-Hastings (RWMH) algorithm discussed in Garthwaite et al. (2016). The current state of the chain ${\bm{\theta}}^{(j)}$ at time $j$ is updated by proposing a draw ${\bm{\theta}}^{\prime}\sim h({\bm{\theta}}|{\bm{\theta}}^{(j)})=\phi_{3}({\bm{\theta}}^{(j)},\tau^{(j)}\Sigma^{(j)})$ where $\phi_{d}({\bm{a}},A)$ denotes a $d$-dimensional Gaussian density function with mean ${\bm{a}}$ and covariance matrix $A$. Because $h({\bm{\theta}}^{\prime}|{\bm{\theta}})=h({\bm{\theta}}|{\bm{\theta}}^{\prime})$ is symmetric the acceptance probability of setting ${\bm{\theta}}^{(j+1)}={\bm{\theta}}^{\prime}$ reduces to

[TABLE]

otherwise ${\bm{\theta}}^{(j+1)}={\bm{\theta}}^{(j)}$. Following Haario et al. (2001), the proposal covariance matrix $\Sigma^{(j)}$ is specified as

[TABLE]

where $\mathrm{I}_{d}$ is the $d$-dimensional identity matrix, $\bar{{\bm{\theta}}}^{(j)}=j^{-1}({\bm{\theta}}^{(1)}+\cdots+{\bm{\theta}}^{(j)})$, and $\tau^{(j)}>0$ is a scaling parameter that affects the acceptance rate of proposal parameter values. Following Garthwaite et al. (2016) we adaptively update $\tau$ using a Robbins-Monro process so that

[TABLE]

where $c=(2\pi)^{1/2}\exp(\zeta_{0}^{2}/2)/(2\zeta_{0})$ is a steplength constant, $\zeta_{0}=-1/\Phi(\pi^{*}/2)$, and where $\Phi$ is the univariate standard Gaussian distribution function. The parameter $\pi^{*}$ is the desired overall sampler acceptance probability, here specified as $\pi^{*}=0.234$ following Roberts et al. (1997). This algorithm is summarised in Step 1 of Algorithm 1. Proposed Gaussian updates for the GEV scale parameter are performed on $\log\sigma$. See Garthwaite et al. (2016) for further details.

3.2 The bivariate case

Inference for the extreme quantile region ${\widetilde{\mathcal{Q}}_{n}}$ in (2.13) requires estimation of the two sets of marginal parameters ${\bm{\theta}}_{i}=(\mu_{i},\sigma_{j},\gamma_{i})$, $i=1,2$, together with the basic set ${\mathcal{S}}$ and its measure $\nu({\mathcal{S}})$. In particular, the estimation of the basic set ${\mathcal{S}}$ its measure $\nu({\mathcal{S}})$ in (2.12) requires estimation of the angular density $h$. We extend the Bayesian procedure of the univariate case (Section 3.1) to the bivariate setting to simultaneously estimate the marginal parameters and the angular density. This framework utilises the censored likelihood based on the bivariate max-stable distribution in (2.8).

Specifically, on the basis of the marginal domain of attraction (see 2.1) and the approximation in (2.3) for each marginal distribution, we define the transformations

[TABLE]

where ${\bm{\theta}}=({\bm{\theta}}_{1},{\bm{\theta}}_{2})$. From the bivariate domain of attraction (2.8) and the max-stability property (de Haan and Ferreira, 2006, Ch 6), for large ${\bm{y}}$, where ${\bm{y}}={\bm{a}}_{n}{\bm{x}}+{\bm{b}}_{n}$, we have

[TABLE]

where ${\bm{\vartheta}}^{*}=({\bm{\theta}}{{}^{*}},h)$, ${\bm{\theta}}^{*}=({\bm{\theta}}_{1}^{*},{\bm{\theta}}_{2}^{*})$, ${\bm{\theta}}_{i}^{*}=(\mu_{i}^{*},\sigma_{i}^{*},\gamma_{i})$, $i=1,2$, and where $h$ is the angular density. Replicating the arguments involving the Pareto distribution for deriving the approximations in (2.6) and (2.7) for each margin, then for ${\bm{y}}\geqslant{\bm{t}}$, where ${\bm{t}}=(t_{1},t_{2})$ is a large threshold, we obtain the further approximation

[TABLE]

where ${\bm{z}}=(z_{1},z_{2})$ with $z_{i}=z_{i}(y_{i};{\bm{\theta}})$ given by (3.4) and where ${\bm{\vartheta}}=({\bm{\theta}},h)$.

Following Marcon et al. (2016) we model the angular density using Bernstein polynomials. Writing $L({\bm{z}})=(z_{1}+z_{2})A(v)$ with $v=z_{2}/(z_{1}+z_{2})$, we model Pickands dependence function $A(v)$ through a Bernstein polynomial of degree $\kappa=0,1,\ldots$ as

[TABLE]

where ${\bm{\beta}}_{\kappa}=(\beta_{1},\ldots,\beta_{\kappa})^{\top}$ is a parameter vector satisfying suitable conditions so that $A_{\kappa}$ in (3.6) defines a valid Pickands dependence function (Marcon et al., 2016, Section 3.1), and $\textrm{Be}(\cdot;a,b)$ is the beta density function with parameters $a>0$ and $b>0$. Modelling Pickands dependence function with a polynomial in Bernstein form is equivalent to modelling the angular distribution with a Bernstein polynomial. The corresponding angular density in Bernstein form is then

[TABLE]

where $w\in[0,1]$ and the elements of ${\bm{\eta}}_{\kappa}=(\eta_{0},\ldots,\eta_{\kappa-2})^{\top}$ must satisfy suitable conditions so that $h_{\kappa-1}$ in (3.7) is a valid angular density (Marcon et al., 2016, Section 3.1). The vectors of coefficients ${\bm{\beta}}_{\kappa}$ and ${\bm{\eta}}_{\kappa}$ are related via a one-to-one relationship (Marcon et al., 2016, Proposition 3.2).

Let $({\bm{y}}_{1},\ldots,{\bm{y}}_{n})$ be a sample of independent bivariate observations from a distribution $F$ on $\mathbb{R}^{2}_{+}$, for which $F\in\mathcal{D}(G)$. Then, exploiting the approximation $\widetilde{G}$ in (3.5) we construct the censored likelihood function

[TABLE]

where each likelihood contribution depends on the domain where ${\bm{y}}_{i}$ falls. Specifically,

[TABLE]

where

[TABLE]

[TABLE]

and where $v=z_{2}/(z_{1}+z_{2})$, $z_{i}(\cdot;{\bm{\theta}})$, $i=1,2$ (3.4), $A^{\prime}_{\kappa}$ and $A^{\prime\prime}_{\kappa}$ are the first and second derivatives of $A_{\kappa}$ with respect to $v$, and where we re-express the full vector of unknown parameters as ${\bm{\vartheta}}=({\bm{\theta}}_{1},{\bm{\theta}}_{2},\kappa,{\bm{\beta}}_{\kappa})$. We infer the marginal GEV parameters and the extremal dependence structure through a semiparametric Bayesian approach. In particular, we combine the inferential scheme for each marginal parameter set ${\bm{\theta}}_{1},{\bm{\theta}}_{2}$ described above with the trans-dimensional MCMC scheme for inferring the dependence structure over the unknown number of elements in ${\bm{\beta}}_{\kappa}$ suggested by Marcon et al. (2016) (see also Antoniano-Villalobos and Walker 2013; Fan and Sisson 2011). The latter takes into account that at each MCMC iteration the dimension of ${\bm{\beta}}_{\kappa}$ changes with $\kappa$ and that the size of ${\bm{\beta}}_{\kappa}$ is potentially infinite (see Marcon et al., 2016, Section 3.2 for details).

While it is more convenient to derive the likelihood function using the representation of extremal dependence through Pickands dependence function (Marcon et al., 2016), it is simpler to define a prior distribution on the angular distribution. In this case, the prior distribution for ${\bm{\beta}}_{\kappa}$ can be deduced from the prior distribution on ${\bm{\eta}}_{\kappa}$ by exploiting the one-to-one relationship between the two (Marcon et al., 2016). Accordingly, the MCMC algorithm may be implemented using the likelihood function parametrised through Pickands dependence function and the prior distribution for this model that is induced from the prior on the angular distribution.

Recall that for all $\kappa\geqslant 3$, ${\bm{\eta}}_{\kappa}$ must satisfy some suitable constraints so that $h_{\kappa-1}$ in (3.7) is a valid angular density. Specifically, the elements of ${\bm{\eta}}_{\kappa}$ must be a nondecreasing sequence in $[0,1]$ such that their sum is equal to $\kappa/2$ (Marcon et al., 2016, Proposition 3.1). Then, following Marcon et al. (2016), for any fixed $\kappa\geqslant 3$, the prior on $({\bm{\eta}}_{\kappa},\kappa)$, with ${\bm{\eta}}_{\kappa}=(\eta_{0},\ldots,\eta_{\kappa-1})$, is defined as

[TABLE]

where $\eta_{0}=p_{0}$ and $\eta_{\kappa-1}=1-p_{1}$ represent the point masses at the endpoints of the simplex ($H(\{0\})$ and $H(\{1\})$), $\Pi(\kappa)=\mathrm{NegBin}(\kappa-3|m_{NB},\sigma_{NB})$ with mean $m_{NB}>0$ and variance $\sigma_{NB}>0$, $\Pi(p_{0})=\textrm{Unif}(0,1/2)$ and $\Pi(p_{1}|\kappa,p_{0})=\textrm{Unif}(a,b)$ with $a=\max\{0,(\kappa-1)p_{0}-\kappa/2+1\}$ and $b=(p_{0}-\kappa/2-1)/(\kappa-1)$. The conditional prior distribution on the remaining parameters is set as

[TABLE]

where the domain $\mathcal{D}_{\bm{\eta}}$ is a suitable set such that the elements of ${\bm{\eta}}_{\kappa}$ satisfy the required conditions above specified (see Marcon et al., 2016, Section 3.2 for details). Notice that, although the theory in Section 2.2 assumes models for which the condition $H(\{0\})=H(\{1\})=0$ holds, our proposed inferential method is capable of estimating point masses at the vertices [math] and $1$, if this exists.

Within the trans-dimensional MCMC update, the pair $({\bm{\eta}}_{\kappa^{(j)}}^{(j)},\kappa^{(j)})$ at time $j$ is updated through the proposal distribution

[TABLE]

where $q_{\bm{\eta}}({\bm{\eta}}_{\kappa}|\kappa)=\Pi({\bm{\eta}}_{\kappa}|\kappa)$ is the conditional prior implied by (3.9), and $q_{\kappa}(\kappa|\kappa^{(j)})$ is defined such that if $\kappa^{(j)}=3$, it places mass on $\kappa=4$ with probability $1$ and if $\kappa^{(j)}>3$ it places mass on $\kappa^{(j)}-1$ and $\kappa^{(j)}+1$ with equal probability. Using $q_{\bm{\eta}}({\bm{\eta}}_{\kappa}|\kappa)=\Pi({\bm{\eta}}_{\kappa}|\kappa)$ means that these terms cancel in the between-model acceptance probability, whether implemented under the ${\bm{\eta}}_{\kappa}$ or ${\bm{\beta}}_{\kappa}$ parameterisation.

The full MCMC sampler is summarised in Algorithm 1, where ${\mathcal{L}}\left({\bm{\theta}}_{1},{\bm{\theta}}_{2},\kappa,{\bm{\beta}}_{\kappa}\right)$ indicates the bivariate censored likelihood (3.8). Separate Robbins-Monro RWMH updates are implemented for each set of marginal parameters ${\bm{\theta}}_{i},i=1,2$, and the above scheme for the dependence parameter updates. See Marcon et al. (2016, Section 3.2 ) for further details.

4 Simulation experiments

4.1 Univariate

We generate $n=1500$ observations from each of three distributions: Fréchet with location, scale and shape parameters equals to $\psi_{0}=3$, $\varsigma_{0}=1$ and $\xi_{0}=1/3$, respectively; $\mbox{Half-}t$ with scale and degrees of freedom equal to $\sigma_{0}=1$ and $\nu_{0}=1/3$; Inverse Gamma with shape parameters equal to $\eta_{0}=1/2$ and $\lambda_{0}=1$. We recall that $\text{Fr}(y;\psi,\varsigma,\xi)=\exp(-((y-\psi)/\varsigma)^{-\xi})$, with $y,\xi>0$, is the Fréchet family of distributions with scale $\varsigma>0$ and location $\psi<y$ parameters. These distributions are in the domain of attraction of a GEV distribution with tail indices $\gamma_{0}=3,3$ and $2$, respectively (Beirlant et al., 2004, Ch. 2). The univariate likelihood function (3.1) is used, censoring observations below the $90$-th empirical quantile. We specify a prior distribution for ${\bm{\theta}}$ as a product of uniform prior distributions on the real line for $\mu$, $\log(\sigma)$ and $\gamma$, i.e. $\Pi({\bm{\theta}})\equiv\Pi(\mu,\sigma,\gamma):=\Pi(\mu)\Pi(\log(\sigma))\Pi(\gamma)$. This improper prior distribution, with $\Pi(\mu,\sigma,\gamma)\propto 1/\sigma$ with $\sigma>0$, leads to a proper posterior distribution (Northrop and Attalides, 2016). We run the MCMC sampler for $M=50,000$ iterations for each dataset.

Each row in Figure 1 shows the estimation results for each dataset: Fréchet (top), Half-$t$ (middle) and inverse Gamma (bottom). The columns on the left present trace plots of the scaling parameter $\tau^{2}$ (initialised at $\tau^{(0)}=1$) and the sampler average acceptance probability. Through the Robbins-Monro process both quantities converge rapidly, with the sampler acceptance rate effectively achieving the target (“optimal”) acceptance rate of $\pi^{*}=0.234$ (solid horizontal line) after no more than $m=30,000$ iterations, which we remove as sampler burn-in. The centre-right panels illustrate histogram and kernel density estimates of the posterior distribution of the tail index $\gamma$, with dashed and solid vertical lines indicating the posterior mean and true value, respectively. The crosses along the horizontal axis show the lower and upper bounds of the estimated $95\%$ credibility interval. In each case the posterior for $\gamma$ puts most of its mass where the true value lies.

Finally, the rightmost panels display the estimated posterior densities of quantiles corresponding to the exceedance probabilities $p=1/750$ (light grey), $1/1500$ and $1/3000$ (dark grey), derived from (2.5). Vertical dashed and solid lines again represent the estimated posterior mean and true quantile values with the latter that are on the log scale: $19.86$, $21.94$, $24.02$ for the Fréchet; $18.73$, $20.81$, $22.89$ for the Half-$t$; $13.48$, $14.87$, $16.25$ for the Inverse Gamma. The corresponding estimated central $95\%$ credible intervals are: $(17.67,22.73)$, $(19.41,25.36)$, $(21.16,27.95)$; $(17.22,22.38)$, $(18.93,25.00)$, $(20.68,27.60)$ and $(12.58,16.58)$, $(13.86,18.58)$, $(15.14,20.57)$, respectively. The points on the $x$-axis are the upper $5\%$ of the observed dataset. As the sample size is $n=1500$, the largest value is a realisation of an event occurring with probability $1/1500$, corresponding to the second investigated quantile. For example, on the log scale, the largest observation from the Fréchet distribution (top row of Figure 1) is $24.36$. This means that the probability of observing an event taking the value $24.36$ or greater is less than $1/750$, and is more likely to be an event with probability closer to $1/3000$.

In all cases, the true tail indices and quantiles are always included inside the estimated central $95\%$ credible intervals, confirming the accuracy of our proposed method for estimating such quantities. In Section 2 of the Supplementary Material we investigate the sensitivity of the above results to variations in the experimental setup. In particular we consider a higher censoring threshold (at the 95-th empirical quantile) and alternative prior specifications for $\mu,\sigma,\gamma$. Performance is qualitatively similar to the above in each case.

4.2 Bivariate

Estimating and quantifying bivariate extreme quantile regions is more challenging than the univariate setting. Here we examine data simulated from two distributions defined on $\mathbb{R}^{2}_{+}$: the positive bivariate Cauchy and a positive truncated bivariate-$t$ density. The former has been previously considered in the literature (e.g. Einmahl et al., 2013; Cai et al., 2011), and we consider the second one as an alternative flexible model. Two further examples (the so-called Asymmetric and Clover densities) are investigated in the Supplementary Material (Section 3). As for the univariate setting, we simulate $n=1500$ observations from each density and marginally censor all observations that fall below the corresponding $90$-th marginal empirical quantile. We are again interested in estimating quantile regions for events with probability $p=1/750$, $1/1500$ and $1/3000$, corresponding to regions that contain very little or no observed data. Details about the selected distributions are:

• •

Bivariate Cauchy distribution on $\mathbb{R}_{+}^{2}$: The Cauchy probability density function is

[TABLE]

The Cauchy distribution is very heavy-tailed with tail indices $\gamma_{1}=\gamma_{2}=1$, and angular density $h(w)=2^{-1}\left(w^{2}+(1-w)^{2}\right)^{-3/2}$. The angular basic density function is $q_{*}(w)=\left(w^{2}+(1-w)^{2}\right)^{1/2}$ and the associated basic set is

[TABLE]

• •

Truncated bivariate-$t$ distribution on $\mathbb{R}^{2}_{+}$: We consider a truncated two-dimensional Student-$t$ distribution with unit-variances, correlation $0\leqslant\rho\leqslant 1$ and $\nu>0$ degrees of freedom. The bivariate-$t$ probability density function is

[TABLE]

where $t_{d,\nu}$ denotes the $d$-dimensional centred Student-$t$ probability density function with correlation $0\leqslant\rho\leqslant 1$ and $\nu>0$ degrees of freedom, $T_{2,\nu}$ is the $d$-dimensional centred Student-$t$ distribution function, and $\bm{0}=(0,0)^{\top}$. The tail indices are $\gamma_{1}=\gamma_{2}=1/\nu$ and vary with the degree of freedom and the angular density is

[TABLE]

The angular basic density is equal to

[TABLE]

and the associated basic set is

[TABLE]

Details concerning the derivation of the angular density are available in Appendix A.1. It is well known that a bivariate-$t$ distribution defined on $\mathbb{R}^{2}$ is in the domain of attraction of the Extremal-$t$ model (see e.g. Beranger and Padoan, 2015). The angular distribution corresponding to this model places mass on all the subsets of the unit simplex, which in the bivariate case corresponds to the interval $(0,1)$ and the vertices $\{0\}$ and $\{1\}$. However, when considering a truncated bivariate-$t$ distribution defined on $\mathbb{R}^{2}_{+}$, it can be shown that the corresponding angular distribution does not put mass at the endpoints, i.e. $H(\{0\})=H(\{1\})=0$ (see Appendix A.1). Note that when $\nu=1$ and $\rho=0$ the truncated bivariate-$t$ distribution reduces to the bivariate Cauchy on $\mathbb{R}^{2}_{+}$.

As before, we specify a uniform prior distribution $\Pi(\mu_{i},\sigma_{i},\gamma_{i})$ for each margin $i=1,2$. As the dependence structures for the above models are known to be smooth (see solid lines in the top left panels of Figures 2–3), it is expected they can be well modelled through relatively low degree polynomials. Hence we set the prior distribution for the polynomial degree as $\Pi(\kappa)=\mathrm{NegBin}(m_{NB}=3.2,\sigma_{NB}=4.48)$. Even though the selected models do not permit mass in the corners of the simplex (i.e. $H(\{0\})=H(\{1\})=0$), in the analysis we still allow for the possibility of non-zero point masses at the endpoints of the simplex by specifying $\Pi(p_{0})=\Pi(p_{1})=\textrm{Unif}(0,0.1)$. This prior is slightly different to those in Marcon et al. (2016) to better represent our prior belief that $p_{0}$ and $p_{1}$ are likely to be small for these data. See Marcon et al. (2016) for an analysis of alternative prior specifications.

We run Algorithm 1 for $M=50,000$ iterations and determine the burn-in period $m$ by visual inspection of trace plots of the marginal scaling parameters $\tau_{1},\tau_{2}$ and the overall acceptance rates of the marginal ($\pi_{i}^{*}=0.234,i=1,2$) and dependence proposals.

For each draw from the posterior distribution, we can construct the angular density $h_{\kappa-1}(w;{\bm{\eta}})$ using (3.7), combine this with the marginal shape indices $\gamma_{1},\gamma_{2}$ into $q(w)$ (2.11), and compute the angular basic density $q_{*}(w)$ (see Section 2.2 for details). Hence, for each fixed $w\in(0,1)$, we obtain samples from the posterior of the angular basic density $q_{*}(w)$. Its inverse $q_{*}^{-1}(w)$ and pointwise central credible intervals are shown in the top-left panels of Figures 2–3.

The boundary of the basic set ${\mathcal{S}}$ corresponds to those points $(x,y)\in(0,\infty)^{2}$ such that $x+y=q_{*}^{-1}(w)$, i.e. the points $\left(wq_{*}^{-1}(w),(1-w)q_{*}^{-1}(w)\right)$. For fixed $w\in(0,1)$, the posterior distribution of this boundary is available, from which we may calculate the posterior mean and 90% central credible intervals. Computed over all $w\in(0,1)$ the (pointwise) posterior mean and credible intervals for the basic set are illustrated in the top-centre panels of Figures 2–3.

To construct the extreme quantile regions for a given small probability $p$ (top-right panel and bottom row of Figures 2–3), consider a point $w\in(0,1)$. As before, we may construct the posterior distribution of the points $(x,y)=\left(wq_{*}^{-1}(w),(1-w)q_{*}^{-1}(w)\right)\in\mathbb{R}^{2}_{+}$ with angle $w$ and radial component $q_{*}^{-1}(w)$, as an estimate of a point on the boundary of ${\mathcal{S}}$. For each $(x,y)$ point in this posterior, we may compute the quantile $\widetilde{Q}_{n}$ associated with probability $p$ via (2.13), leading to a posterior distribution approximating the $1-p$ bivariate quantile level at a particular (estimated) point in ${\mathcal{S}}$. We then use the marginal $0.05$ and $0.95$ quantiles of this bivariate posterior to define an approximate $90\%$ credible region, and the marginal posterior means to produce a mean estimate. This procedure is repeated for other $w\in(0,1)$.

The top right panel in Figures 2–3 provides a comparison between the true extreme quantile region (solid line) and the posterior mean estimate (dotted line), for $p=1/750,1/1500$ and $1/3000$, with the observed data (points) overlaid. The bottom panels in each figure illustrate each extreme quantile region separately, but with the estimated 90% pointwise credible regions and the point estimates of the quantile regions given in Einmahl et al. (2013) (dotted lines) are included for comparison. We refer to the latter as the EdHK estimator.

Figure 2 illustrates the results for the Bivariate Cauchy distribution. The estimated (pointwise posterior mean) inverse angular basic density $q_{*}^{-1}(w)$ (dotted line; top-right panel) describes the behaviour of the true function (solid line) well, and is fully included in the (pointwise) estimted $90\%$ credible regions. The curvature in the dependence structure around $w=0.5$ is not as pronounced as the true $q_{*}^{-1}(w)$, which consequently impacts on the estimated basic set and extreme quantile regions. The bottom panels show that the mean extreme quantile levels (dotted) line are close to the true levels (solid line) and provide a similar fit to the EdHK point estimate (dashed) for each probability $p=1/750,1/1500$ and $1/3000$. In addition all three curves consistently appear in the centre of the $90\%$ credible regions.

Finally, Figure 3 presents analysis for the truncated bivariate-$t$ distribution. The dependence structure $q_{*}^{-1}(w)$ is very accurately estimated within the interior of the simplex, although towards the endpoints the estimated $q_{*}^{-1}(w)$ appears to approach [math] whereas the true value is approximately $1$. Producing quantile regions that seem drawn to the origin for low values of either component, may at first appear erroneous. However, inspection of the simulated dataset (top-right panel) reveals that there are no extreme observations that lie along the axes, thereby explaining the behaviour of our estimator. On balance, the estimated posterior means follow the general shape of the true quantile regions and provide similar results than the competitor EdH estimates (dashed lines). The true quantile regions are amply included in the (pointwise) estimated $90\%$ credible regions.

Overall, the proposed methodology is able to accurately estimate both marginal and dependence structure of process extremes. As a point estimator, the posterior mean of the extreme quantile region performs well in recovering the true region, while being responsive to the dependence structure within the data itself, as with the bivariate-$t$ data analysis (Figure 3). It also seems to perform more consistently than the EdHK estimator (e.g. Figures 5 and 6 of the Supplementary Material). The credible regions both provide some measure of the uncertainty inherent with low probability events, while also allowing for judgements regarding whether exceptionally high observations can still be considered to belong to events with particular probabilities (e.g. whether the single large outlier in Figure 3 can be considered a $1/1500$ or $1/3000$ probability event).

5 Analysis of extreme air pollution levels in Milan

Understanding the extremes of air pollutants is of critical importance, especially in the context of climate change (De Sario et al., 2013). In recent work, Martins et al. (2017) used standard univariate extreme value theory to compare the air quality between the two largest urban regions in Brazil, Heffernan and Tawn (2004) and Beranger and Padoan (2015) estimated the extremal dependence between pollutants in Leeds, U.K., Coles and Pan (1996) studied the extreme temporal behaviour of airborne NO2 particles in Milan between 1984–1994, and Vettori et al. (2019) performed a spatial analysis of air pollution in Los Angeles, U.S.A.

Here we study extreme air pollution levels recorded in Milan, Italy, over the winter period October 31st–February 27/28th, between December 31st 2001 and December 30th 2017. We examine the daily mean level of PM10 and daily maximum levels of NO, NO2 and SO2. These data were also considered by Falk et al. (2019, Table 3) although there the objective was estimating joint exceedance probabilities (with excesses in all variables).

Here, the aim is to estimate bivariate extreme quantile regions given by a high concentration of the two pollutants and with a small probability of occurrence. This analysis is important for understanding the long term impact of air quality on health. When modelling univariate extremes it is common to express marginal parameters via a regression model which may be functions of space or other covariates (Padoan et al., 2010). For air pollution in particular, evidence is overwhelming that interactions between the pollution and temperatures cannot be disregarded (De Sario et al., 2013; Cheng and Kan, 2012; Katsouyanni et al., 1993; Roberts, 2004). The leftmost panels of Figure 4 show scatterplots of each pollutant against the daily maximum temperature. All four indicate a possible quadratic relationship between pollutant and maximum temperature, and so we write the $i$-th marginal mean $\mu_{i}=\beta_{0,i}+\beta_{1,i}z+\beta_{2,i}z^{2}$ as a quadratic function of the maximum temperature $z$. Other covariates could have been included if they were available, and regressions on $\sigma_{i}$ and $\gamma_{i}$ constructed, although we did not pursue that here. For observations that fall below the threshold (black points in left panels of Figure 4), as these observation are censored it could be argued that the level of the covariates at the threshold should be used when evaluating the likelihood contribution. This would beneficially reduce computational costs. However, as the covariates are still available for censored observations, they still provide valuable information to estimate the regression coefficients.

Using mean residual life plots (e.g. Beirlant et al., 2004, Section 1.2.2) we set the marginal $90$-th empirical quantile as the threshold (points above the threshold are blue in Figure 4), which results in $t=642.9$ ($n=1796$) for NO, $t=139.5$ ($n=1799$) for NO2, $t=108$ ($n=1779$) for PM10 and $t=45.2$ ($n=1809$) for SO2. These thresholds are comparable with those in Falk et al. (2019, Table 3). We specify all marginal prior distributions $\Pi(\mu_{i},\sigma_{i},\gamma_{i})\propto 1/\sigma_{i}$, for $i=1,\ldots,4$ as in Section 4.1, and dependence parameter prior distributions as in Section 3.2 with $\Pi(\kappa)=\mbox{NegBin}(\kappa=3|m_{NB}=6,\sigma_{NB}=8)$ to allow for higher degree polynomial modelling of $h_{\kappa-1}(w;{\bm{\eta}})$ if required. For the univariate analysis we implement an MCMC sampler with 300k iterations and retain the final 50k iterations for inference; for the multivariate analysis we retain the final 20k iterations from a chain of length 50k. The univariate sampler ran in under 1 hour on a single-core 2.2 GHz Intel Core i7 processor on a MacBook Pro; the multivariate sampler took around 2.5 hours.

The image plots in Figure 4 illustrate the (univariate) posterior distribution of the quantiles associated with the probabilities $p=1/600,1/1200$ and $1/2400$ (left to right) as a function of the maximum daily temperature, which correspond to pollutant levels that would be expected to be reached once every 5, 10 and 20 winters. Restricting our interpretation to the range of observed temperature, we note that extreme PM10 quantiles (third row) appear higher for lower temperatures rather than high. Indeed, the main sources of PM10 pollution include combustion engines (both diesel and petrol) and combustion for energy production in households. Accordingly, when maximum temperatures are low, one may expect household heating systems to work at higher capacity and an increase in the use of cars rather than less energy consumptive transport methods such as walking or cycling.

Extreme quantile behaviour for NO2 (second row) as function of maximum temperature appears similar to PM10 although the larger range of the quantiles ($y$-axis) reduces the visual curvature. This behaviour can be explained by the fact that nitrogen dioxide, mainly emitted by power generation, industrial and traffic sources, is an important constituent of particulate matter. Sulfur dioxide (fourth row) forms either naturally (decomposition or combustion of organic matter) or due to human activity (smelting sulfur-containing mineral ores). Accordingly, it is not surprising to observe decreasing extreme SO2 quantiles across temperature levels (fourth row) as maximum temperature moves from [math] to $20^{\circ}\mathrm{C}$. The decreasing quantile behaviour as temperature gets cold might be due to a reduction in human activity. The posterior mean of the tail indices \gamma_{\textrm{NO}},\gamma_{\textrm{NO{}{2}}},\gamma_{\textrm{PM{}{10}}} and \gamma_{\textrm{SO{}_{2}}} and their corresponding $90\%$ credible interval are $0.050$ $(0.003,0.138)$, $0.186$ $(0.055,0.337)$, $0.054$ $(0.004,0.148)$ and $0.064$ $(0.005,0.167)$ respectively.

To date, the link between extreme levels of multiple pollutants and human health has not been well explored, although a multi-pollutant approach to assess the health risks associated with air pollution has been emerging (e.g. Dominici et al., 2010; Wesson et al., 2010). In the following we restrict our attention to the three pollutants with the largest tail indices: NO2, SO2 and PM10.

Figure 5 shows information regarding extremal dependence between the three pollutant pairs: NO2/SO2, NO2/PM10 and SO2/PM10. The estimated basic sets demonstrate weak dependence between the pairs NO2/SO2 and SO2/PM10 and stronger dependence between NO2 and PM10. Further, from the pairwise scatterplots it is apparent that neither the coldest (blue dots) nor the warmest (red dots) daily maximum temperatures appear to induce the largest pollutant levels. As discussed above, NO2 is a key constituent of PM10 and so it is realistic to expect the observed strong dependence between these pollutants. Similarly, the sources of PM10 (combustion engines) and SO2 (natural or smelting mineral ores) differ, explaining the independence between the extremes of these pollutants to some extent.

Figure 6 illustrates extreme quantiles regions corresponding to events that would expect to occur once every 5, 10 and 20 winters (left to right panels) when the daily maximum temperature is fixed to the minimum, median and maximum observed daily maximum temperatures (blue to red shading). These values are $(-6.3^{\circ}\textrm{C},8.6^{\circ}\textrm{C},22.2^{\circ}\textrm{C})$, $(-5.1^{\circ}\textrm{C},8.6^{\circ}\textrm{C},22.2^{\circ}\textrm{C})$ and $(-5.1^{\circ}\textrm{C},8.7^{\circ}\textrm{C},22.2^{\circ}\textrm{C})$ for the pairs NO2/SO2, NO2/PM10 and SO2/PM10 respectively. The top panels of Figure 6 exhibit a small quadratic influence of maximum daily temperature on the joint levels of NO2 and SO2, and similarly for SO2 and PM10 (bottom row). It is apparent that extreme levels of SO2 are reduced with warmer weather. This phenomena can be loosely understood by the fact that sulfur dioxide is an aerosol which cools the planet by reflecting some of the sun’s energy back into space. As such, large levels of SO2 are more likely to be associated with cold temperatures. Overall the estimated extreme quantile regions capture the behaviour of the data well, with the NO2/PM10 pair exhibiting the strongest level of dependence. As expected, as the probability $p$ of an event decreases the quantile regions become larger and reach higher pollutant levels. Despite having observations for 17 consecutive winters, our method is able to provide quantile levels and credible regions for events with probability of occurrence $p$ at arbitrarily low levels.

In terms of the pollutant thresholds suggested by European emission regulation, the top-left panel of Figure 6 shows that when the the daily maximum temperature is fixed to the maximum, we expect that NO2 individually exceeds approximatively $200$ $\mu$g/m3 (i.e. the short-term limit threshold) once every $5$ winters ($p=1/600$) on average. However, this event is as equally probable as the joint event that both NO2 and PM10 each exceed approximatively $200$ $\mu$g/m3. Such a concentration for PM10 is four times higher than the individual limit threshold for this pollutant (see Section B.1). Similarly, for the bottom-left panel of Figure 6, when the daily maximum temperature is fixed to the median temperature, we expect that PM10 individually exceeds approximatively $200$ $\mu$g/m3 on average once every $5$ winters. However, this event is as equally probable as the joint event that both SO2 and PM10 each exceed approximatively $200$ $\mu$g/m3. This concentration for SO2 is almost twice the individual limit threshold for this pollutant (see Section B.1). These results indicate both the strong dependence between harmful pollutants at extreme levels in Milan, and that limit threshold alerts for poor air quality in this city are likely to be issued for multiple pollutants simultaneously.

6 Discussion

We have presented a new method for estimating extreme quantile regions that is responsive to varying levels of extremal dependence, and comes with natural measures of uncertainty, both for model parameters and extreme quantile regions, under the Bayesian paradigm. The method was able to outperform the existing (EdHK) approach of Einmahl et al. (2013) which does not provide any measure of uncertainty. This methodology provides a useful general tool for long-term analysis of the air pollution at the extreme level. We explored this in Section 5 for assessing and quantifying the health risks associated with multiple extreme pollutant exposures (Dominici et al., 2010; Wesson et al., 2010).

The methods developed in this paper have been incorporated into accessible functions in the R Package ExtremalDep. This package is freely available from the CRAN repository, see https://CRAN.R-project.org/package=ExtremalDep. The R code for the simulations and the multivariate analysis of real air pollution data is available online at the address https://www.borisberanger.com/zip/BPS_2019.zip.

Acknowledgements

The authors are grateful to Andrea Krajina for sharing the code for the frequentist estimation of bivariate extreme quantile regions and her valuable suggestions and help. SAP is supported by the Bocconi Institute for Data Science and Analytics (BIDSA), Italy and PRIN 2015 research grant Modern Bayesian Nonparametric Methods. SAS and BB are supported by the Australian Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS; CE140100049) and the Australian Research Council Discovery Project scheme (FT170100079). The authors are also indebted to the Associate Editor and two anonymous reviewers for their careful reading of the manuscript and their constructive remarks

Appendix A Appendix

A.1 The Extremal-$t$ model with restriction to the positive reals

Consider a student-$t$ distribution restricted to $(0,\infty)$. Using Beirlant et al. (2004, p.59) the norming constants required in (2.1) are

[TABLE]

Applying the conditional tail dependence function framework of Nikoloulopoulos et al. (2009), the exponent function can be written as:

[TABLE]

where $(Z_{1},Z_{2})^{\top}$ follows a centred bivariate-$t$ distribution on $(0,\infty)^{2}$ with unit variance, correlation $\rho$ and degree of freedom $\nu$. The conditional distribution of $Z_{2}|Z_{1}=z_{1}$ is a truncated $t$ distribution on $(0,\infty)$ with mean $\rho z_{1}$, variance $(\nu+z_{1}^{2})(1-\rho^{2})/(\nu+1)$ and $\nu+1$ degrees of freedom. As a consequence, we obtain

[TABLE]

and

[TABLE]

Identical calculations can be applied to the second term in the exponent function. Transforming the margins to unit-Fréchet margins allows expression of the exponent function as

[TABLE]

It is easy to verify that $\lim_{y\to 0}\partial/\partial xV(x,y)=0$ and $\lim_{x\to 0}\partial/\partial yV(x,y)=0$ which implies $H(\{0\})=H(\{1\})=0$. Finally, note that due to the form of $V(x,y)$, taking the double partial derivative with respect to $x$ and $y$ is equivalent to the double partial derivative of the exponent function of the Extremal-$t$ model multiplied by a scaling term. Hence the angular density on the interior of the $2$-dimensional unit simplex is as given.

Appendix B Supplementary Material for “Estimation and uncertainty

quantification for extreme quantile regions”

B.1 Introduction

Note that all references below of the form ($\cdot^{*}$) refer to equation ($\cdot$) in the main paper. We recall that “EdHK” is the abbreviation for Einmahl-de Haan-Krajina with reference to the results developed in Einmahl et al. (2013)

B.2 Complement to Section 4.1

We first replicate the simulation experiment in Section 4.1 of the main article, but this time censoring observations below the $95$-th empirical quantile. We simulate $n=1500$ observations from the distributions: Fréchet with location, scale and shape parameters equal to $\psi_{0}=3$, $\varsigma_{0}=1$ and $\xi_{0}=1/3$, respectively; $\mbox{Half-}t$ with scale and degrees of freedom equal to $\sigma_{0}=1$ and $\nu_{0}=1/3$; Inv-Gamma with shape parameters equal to $\eta_{0}=1/2$ and $\lambda_{0}=1$. The results are illustrated in Figure 7, which is organised the same as Figure 1 in the main text.

Each row shows the results for each dataset: Fréchet (top), Half-$t$ (middle) and inverse Gamma (bottom). The left columns report the trace plots of the scaling parameter $\tau^{2}$ and the sampler average acceptance probability. The centre-right panels show histogram and kernel density estimates of the posterior distribution of the tail index $\gamma$, with dashed and solid vertical lines indicating the posterior mean and true value, respectively. Finally, the rightmost panels show the estimated posterior densities of quantiles corresponding to the exceedance probabilities $p=1/750$ (light grey), $1/1500$ and $1/3000$ (dark grey), derived from (2.4*∗*). Vertical dashed and solid lines represent the posterior mean and true quantile values.

The conclusions are similar to those obtained by censoring observations below the $90$-th empirical quantile; see the main text for details. Once again, the true tail indices and quantiles always fall within the $95\%$ credible intervals obtained from the estimated posterior distributions. Using a higher threshold leads to the following differences with the results described in the main article: for the Fréchet distribution (top row) the posterior means are all greater, and for the inverse-Gamma (bottom row) the posterior means are all greater. The posterior means are almost unchanged for the Half-$t$ distribution (middle row).

We additionally perform a sensitivity analysis to examine the sensitivity of the posterior to the prior distribution. For each univariate distribution considered in the main article (i.e. Fréchet, $\mbox{Half-}t$ and Inv-Gamma) we simulate $n=1500$ observations and then censor these below the $90$-th empirical quantile. We consider the following three prior distributions:

• A)

$\Pi(\mu,\sigma,\gamma):=\Pi(\mu)\cdot\Pi(\log(\sigma))\cdot\Pi(\gamma)$, where the prior distributions for $\mu$, $\log(\sigma)$ and $\gamma$ are uniform on the real line;

• B)

$\Pi(\mu,\sigma,\gamma):=\mathcal{N}\cdot(\mu;0,2^{2})\cdot\mathrm{log}\mathcal{N}(\sigma|0,(1/3)^{2})\cdot\mathcal{N}(\gamma;0,(3/2)^{2})$, where the $\mathcal{N}(x;a,b^{2})$ and $\mathrm{log}\mathcal{N}(x;a,b^{2})$ denote normal and lognormal distribution with mean variance parameters $a$ and $b^{2}$, respectively;

• C)

$\Pi(\mu,\sigma,\gamma):=\mathcal{N}(\mu;0,5^{2})\cdot\mathrm{log}\mathcal{N}(\sigma;0,5^{2})\cdot\mathcal{N}(\gamma;0,6^{2})$.

Prior A) is the improper prior used in the main text, prior B) is strongly informative, and prior C) is proper but weakly informative. Each MCMC sampler is then run for $M=50,000$ iterations for each simulated dataset. The results are shown in Figures 8–10, with columns (left-to-right) displaying the results under prior distributions A)–C).

The strongly informative prior distribution B) leads to estimated posterior distributions that are heavily influenced by the prior. As a result, the point and interval estimates of the tail index and quantiles are similarly driven by these beliefs, and so these do not correspond to the true values. Of course, if the practitioner truly believes this prior specification, then the resulting predictive inference is strictly correct in the usual Bayesian sense.

Under priors A) and C) the strength of prior information is much weaker, and as a result the posterior and resulting inferences are driven more by the data. The resulting point and interval estimates naturally accord more strongly with the true values in this case.

Appendix C Complement to Section 4.2

We present two extra simulation experiments, involving estimating bivariate extreme quantile regions from the so-called Asymmetric and Clover densities. We simulate $n=1500$ observations from each of the following bivariate densities:

• •

Asymmetric distribution on $\mathbb{R}_{+}^{2}$: The probability density of the asymmetric distribution is

[TABLE]

with $c\approx 0.58$. This density is heavy-tailed but its upper tail is less heavy than that of the bivariate Cauchy. Its tail indices are $\gamma_{1}=4/5$ and $\gamma_{2}=3/5$, and the angular density is given by

[TABLE]

where $c_{1}\approx 0.589$ and $c_{2}\approx 0.593$. The angular basic density is

[TABLE]

and the associated basic set is

[TABLE]

• •

Clover distribution on $\mathbb{R}_{+}^{2}$: The probability density of the Clover distribution is

[TABLE]

This density is heavy-tailed and its joint upper tail is heavier than that of the the bivariate Cauchy. Its tail indices are $\gamma_{1}=1$ and $\gamma_{2}=5/4$, and the angular density is

[TABLE]

where $c\approx 0.208$. The angular basic density is

[TABLE]

and the associated basic set is

[TABLE]

For each margin we censor all observations that fall below the corresponding $90$-th marginal empirical quantile. We specify the prior distribution $\Pi(\mu_{i},\sigma_{i},\gamma_{i})\propto 1/\sigma$ for each margin $i=1,2$ and a $\Pi(\kappa)=\mathrm{NegBin}(m_{NB}=3.2,\sigma_{NB}=4.48)$ prior distribution for the polynomial degree. We allow for the possibility of non-zero point masses at the endpoints of the simplex by specifying the uniform prior distributions $\Pi(p_{0})=\Pi(p_{1})=\textrm{Unif}(0,0.1)$. We run Algorithm 1 in the main text for $M=50,000$ iterations and determine the burn-in period $m$ by visual inspection of trace plots of the marginal scaling parameters $\tau_{1},\tau_{2}$ and the overall acceptance rates of the marginal ($\pi_{i}^{*}=0.234,i=1,2$) and dependence proposals. Then, we estimate the bivariate quantile regions corresponding to the exceedance probabilities $p=1/750$, $1/1500$ and $1/3000$.

The complete description on how to derive the estimated posterior densities for the angular basic density $q_{*}(w)$, the basic set ${\mathcal{S}}$ and the quantile $\widetilde{Q}_{n}$ and their respective posterior means and credibility bands is given in Section 4.2 of the main text.

The results for the asymmetric distribution (Figure 11) indicate a generally well captured dependence structure (top-left panel), albeit with an estimated dip down towards [math] at the boundaries of the simplex. This corresponds to a basic set where values of one component that are close to [math] also tend to drag the other component towards [math]. The top right panel illustrates that the quantile regions follow the behaviour of the data quite well, particularly the largest observations. The mean quantile estimates (bottom panels) appear to closely match the true levels, while the EdHK estimates do not capture either the shape or the magnitude of the region well, and do not convincingly lie within the 90% credible quantile regions.

The structure of the clover distribution (Figure 12) is more challenging to model. Our procedure, despite having some difficulty capturing the detailed fluctuations for the given observed dataset size ($n=1500$), appears to estimate the general trend of the dependence behaviour fairly well. Indeed, the inverse of the true angular basic density $q_{*}^{-1}(w)$ and the true basic set ${\mathcal{S}}$ are almost fully contained within the relatively narrow 90% credible regions. The EdHK estimator performs reasonably poorly, giving higher quantile estimates than our mean quantile estimate, with both overestimating the true quantile level, although the overestimation is only slight for our mean estimator. This can be explained by the presence of an observation above the $(2999/3000)$ quantile and two observations above the $(1499/1500)$ quantile.