# Semiparametric Estimation of the Shape of the Limiting Bivariate Point Cloud

Reetam Majumder
SECASC
North Carolina State University
Raleigh, NC, 27695
Benjamin A. Shaby
Department of Statistics
Colorado State University
Fort Collins, CO, 80523
Brian J. Reich
Department of Statistics
North Carolina State University
Raleigh, NC, 27695
Daniel Cooley
Department of Statistics
Colorado State University
Fort Collins, CO, 80523

###### Abstract

We propose a model to flexibly estimate joint tail properties by exploiting
the convergence of an appropriately scaled point cloud onto a compact limit
set. Characteristics of the shape of the limit set correspond to key tail
dependence properties. We directly model the shape of the limit set using
Bézier splines, which allow flexible and parsimonious specification of shapes
in two dimensions. We then fit the Bézier splines to data in pseudo-polar
coordinates using Markov chain Monte Carlo sampling, utilizing a limiting
approximation to the conditional likelihood of the radii given angles. By
imposing appropriate constraints on the parameters of the Bézier splines, we
guarantee that each posterior sample is a valid limit set boundary, allowing
direct posterior analysis of any quantity derived from the shape of the curve.
Furthermore, we obtain interpretable inference on the asymptotic dependence
class by using mixture priors with point masses on the corner of the unit box.
Finally, we apply our model to bivariate datasets of extremes of variables
related to fire risk and air pollution.

Keywords: Bayesian inference, Bézier curves, extreme values, gauge function,
limit set

##  1 Introduction

Multivariate tail risk calculations require knowledge of the strength of
dependence in the joint tail of the relevant distribution. Here, we propose a
model to flexibly estimate joint tail characteristics in a way that coherently
links several existing measures of tail dependence. To do this, we describe
tail dependence of a multivariate distribution through its associated _gauge
function_ (Balkema et al., 2010; Balkema and Nolde, 2010; Nolde, 2014). The
homogeneity property of the gauge function allows us to recover the entire
gauge function from its unit level set, which bounds the support of the
appropriately scaled data points in the limit (Nolde and Wadsworth, 2022;
Wadsworth and Campbell, 2024). We represent the unit level set of the gauge
function using a semiparemtric model, specified such that the required
constraints on such functions are automatically satisfied. In this way, we
obtain a posterior sample of gauge functions, with each member of the sample
being a valid gauge function not requiring any _post hoc_ adjustments such as
re-scaling or truncation.

Efforts to exploit the limit set representation of multivariate extreme values
(Davis et al., 1988; Kinoshita and Resnick, 1991; Balkema et al., 2010;
Balkema and Nolde, 2010) have appeared only recently. Wadsworth and Campbell
(2024) decompose the data into pseudo-polar coordinates and use a limiting
argument to approximate the distribution of the radii with a truncated Gamma
distribution whose parameters depend on the gauge function. They first
transform the data to unit exponential margins, and assuming a parametric form
for the gauge function, perform maximum likelihood estimation with the
truncated Gamma likelihood. They extend this approach using mixtures of
parametric forms, but need to perform post-hoc re-scaling of the mixtures to
satisfy the required properties of valid gauge functions.

In contrast to Wadsworth and Campbell (2024), whose primary focus is
estimating probabilities of sets in the joint tail region, Simpson and Tawn
(2022) focus on inference for the limit set boundary itself, taking a flexible
semiparametric approach. They estimate the sample limit set by approximating
the limiting upper endpoint of the distribution of radii with an estimated
high quantile, as a function of angle. To do this, they fit a generalized
Pareto distribution, whose scale parameter varies by angle, to the large
radii. The radii are calculated by decomposing the bivariate data points
transformed to unit exponential margins with a rank transformation. As the
result is not a valid limit set, they perform a subsequent scaling and
truncation procedure based on a Hill estimator (Hill, 1975) to force their
estimate to satisfy the required conditions.

Like Simpson and Tawn (2022), our focus is flexible estimation of the limit
set boundary, though our methodology is quite different. Here, we directly
model the boundary of the limiting scaled point cloud, which is prescribed by
the unit level set of the gauge function, as a Bézier spline. Bézier splines
are constituted of Bézier curves, with points on the curve that can be
represented as Bernstein polynomials (for reviews, see Hazewinkel, 2012;
Farouki, 2012). Similar semiparametric approaches have been used previously in
multivariate extremes to characterize the Pickands dependence function for
extremal dependence (Marcon et al., 2014, 2017; Vettori et al., 2018), the
angular density in the context of multivariate regular variation (Hanson et
al., 2017), and the angular dependence function (Murphy-Barltrop et al.,
2023), which we will explore below as a direct byproduct of the gauge
function. Bézier splines are convenient here because they allow parsimonious
specification of shapes in $\mathbb{R}^{2}$ which are defined by a small
number of control points. Placing appropriate constraints on the control
points can ensure that the resultant shapes satisfy the conditions required of
limit set boundaries. To estimate the parameters of the Bézier spline, we use
the result from Wadsworth and Campbell (2024) which says that, given a gauge
function evaluated at the data points, the distribution of the large radial
components decays like a Gamma distribution whose rate parameter depends on
the gauge function. We then use standard Markov chain Monte Carlo machinery to
sample from the posterior distribution.

Our approach has several advantages. First, we model the shape of the limiting
point cloud in a way that automatically results in a valid limit set, without
the need for _post hoc_ fixes. Second, our model allows the boundary of the
limit set to exactly touch the corners of the unit box; this in particular
gives a clean interpretation of the distinction between asymptotic
independence (AI) and asymptotic dependence (AD) classes, since this
distinction corresponds to whether or not the boundary touches the upper right
corner. Third, our approach produces a posterior sample of valid limit set
curves which yields a realistic picture of the state of knowledge about the
joint tail region given the data. In addition, we also note that our work
builds on a growing literature on Bayesian approaches to multivariate extreme
value analysis (e.g. Boldi and Davison, 2007; Sabourin et al., 2013; Sabourin
and Naveau, 2014; de Carvalho et al., 2022; Padoan and Rizzelli, 2022, to name
a few).

The rest of the paper is arranged as follows. Section 2 introduces the
limiting scaled point cloud and how it can be used to model the tail behavior
of bivariate data. Section 3 develops the modeling of the limit set boundary
using Bézier splines. Section 4 contains a simulation study demonstrating our
approach. Section 5 contains two applications—the Santa Ana Winds dataset, and
ozone concentration data for the contiguous US—where we use Bézier splines to
model the tail dependence in the data. Section 6 concludes.

##  2 The Limiting Scaled Point Cloud

Consider a collection of $n$ independent random vectors in
$\mathbb{R}^{2}_{+}$, $\bm{X}_{1},\ldots,\bm{X}_{n}$, each having joint
density $f_{\bm{X}}$, with standard exponential margins. At times, it will be
convenient to transform the components of $\bm{X}=(X_{1},X_{2})^{\text{T}}$
into pseudo-polar coordinates $(R,\bm{W})$, as $R=X_{1}+X_{2}$, and
$\bm{W}=\bm{X}/R$. Note that for $\bm{W}=(W_{1},W_{2}),W_{2}=1-W_{1}$.

Now define the scaled point cloud as the collection of points divided by
$\log{n}$, $\\{\bm{X}_{1}/\log{n},\ldots,\bm{X}_{n}/\log{n}\\}$. If we assume
that $\lim_{t\rightarrow\infty}-\log f_{\bm{X}}(t\bm{x})/t=g(\bm{x})$,
$\bm{x}\in\mathbb{R}^{2}_{+}$, for some continuous function $g$, then the
scaled point cloud converges onto a compact limit set

$G=\\{\bm{x}\in\mathbb{R}^{2}:g(\bm{x})\leq 1\\}$

as $n\rightarrow\infty$ (Davis et al., 1988; Kinoshita and Resnick, 1991;
Nolde, 2014; Nolde and Wadsworth, 2022). The function $g$ is called the _gauge
function_ associated with the density $f_{\bm{X}}$. Denote the boundary of $G$
as $\partial G$.

Every gauge function $g$ is homogeneous of order one, with
$g(c\bm{x})=cg(\bm{x})$ for any $c>0$ (Nolde, 2014). We will exploit this
property by modeling the limit set boundary $\partial G$ directly and using
its associated gauge function, induced by homogeneity, for estimation (see
Section 3.2). Any valid limit set $G$ must satisfy the following constraints
on its shape:

  1. 1.

$G$ is _star-shaped_ , meaning that for any $t\in(0,1)$, if $\bm{x}$ is in
$G$, then $t\bm{x}$ is also in $G$.

  2. 2.

The supremum of the boundary $\partial G$ is 1 in each component direction.
That is, $\partial G$ touches, but does not cross, the upper and right-hand
sides of the unit box.

Figure 1: Schematic of $\eta$. The blue curve is the unit level set of gauge
function $g(\bm{x})$, which forms the limit set boundary $\partial G$, and the
proportional distance to the red point from the origin is the tail dependence
coefficient $\eta$. While the red point is always on the diagonal, the
intersection of the shaded red region and the blue curve does not necessarily
occur on the diagonal.

We seek a flexible way of representing the boundary $\partial G$ of the limit
set $G$ that satisfies conditions 1 and 2 and can be estimated from iid
samples of the random vector $\bm{X}$. The shape of the limit set contains
useful information about the extremal dependence of the distribution of the
data. Nolde and Wadsworth (2022) linked particular features of the shape of
$G$ with various indices of joint tail dependence in the literature. The
residual tail dependence coefficient (Ledford and Tawn, 1996), the angular
dependence function (Wadsworth and Tawn, 2013), components of the conditional
extremes model (Heffernan and Tawn, 2004), and the index $\tau_{1}(\delta)$
(Simpson et al., 2020) all have direct connections to the shape of $G$. Our
primary focus is on the residual tail dependence coefficient, $\eta\in(0,1]$,
which is defined by assuming that, for $\bm{X}$ in exponential margins, its
survivor function satisfies

$P(X_{1}>x,X_{2}>x)=\mathcal{L}(e^{x})e^{-x/\eta}$

as $x\rightarrow\infty$, for some function $\mathcal{L}$ that is slowly
varying at infinity (Ledford and Tawn, 1996). Then the coefficient $\eta$
describes the strength of dependence in the joint tail, with $\eta\in(1/2,1)$
indicating positive tail dependence but AI, and $\eta=1$ indicating AD,
assuming $\mathcal{L}(x)\nrightarrow 0$. The dependence class (AI vs. AD) is
defined by the limiting conditional probability $\chi$, where

$\lim_{x\rightarrow\infty}\frac{P(X_{1}>x,X_{2}>x)}{P(X_{1}>x)}>0,$

with $\chi=0$ characterizing AI, and $\chi>0$ characterizing AD.

The residual tail dependence coefficient, $\eta$, can be calculated (Nolde,
2014; Nolde and Wadsworth, 2022) from shape of the limit set as

$\eta=\min\\{r:r\times[1,\infty]^{2}\cap G=\emptyset\\}.$

This is illustrated schematically in Figure 1, where one can think of sliding
the shaded box down the ray with slope 1 until it first touches the boundary
$\partial G$. The radius corresponding to this first point of intersection is
$\eta$. A corollary is that (assuming as above that
$\mathcal{L}(x)\nrightarrow 0$) when $\bm{X}$ is AD, $\eta=1$, so $\partial G$
necessarily touches the upper right-hand corner of the unit box. Conversely,
when $\bm{X}$ is AI, $\eta<1$, so $\partial G$ does not touch the upper right-
hand corner, and is referred to as _blunt_.

##  3 Modeling the Shape Using Bézier Splines

###  3.1 A Bézier spline representation of the limit set boundary

Figure 2: Examples of Bézier curves of orders 1, 2, and 3. The red control
points (end points) $\mathbf{p}_{0}$ and $\mathbf{p}_{m}$ always lie on the
curve, while the blue control points usually do not.

Bézier curves (e.g. Hazewinkel, 2012; Farouki, 2012) are a class of parametric
functions that can be used as building blocks to represent complex shapes.
Bézier curves are defined by a set of control points $\mathbf{p}_{0}$ to
$\mathbf{p}_{m}$, where $m$ is the order of the curve. Figure 2 plots examples
of Bézier curves of orders 1–3. The end points (red) define the beginning and
end of the curve; intermediate control points (blue) of each curve control its
shape but generally do not lie on the curve. A quadratic Bézier curve, for
example, traces the path:

$B(t)=(1-t)[(1-t)\mathbf{p}_{0}+t\mathbf{p}_{1}]+t[(1-t)\mathbf{p}_{1}+t\mathbf{p}_{2}],$

for $0\leq t\leq 1$. Rearranging this equation simplifies it to:

$B(t)=(1-t)^{2}\mathbf{p}_{0}+2t(1-t)\mathbf{p}_{1}+t^{2}\mathbf{p}_{2}.$

A useful property is that if the three points are co-linear, then a quadratic
Bézier curve simplifies to a linear Bézier curve. Several Bézier curves can in
turn be linked together at the end points to form a Bézier spline. The end
points of each Bézier curve within the spline now function as knots for the
spline. Splines comprised of quadratic Bézier curves are particularly useful
since analytical solutions for quadratic equations are straightforward to
obtain. In addition, increasing the order to cubic splines would make it
difficult to constrain the shapes to the unit box, and would prevent the
shapes from having the sharp corners required to represent AD limit set
boundaries.

Figure 3: Examples of unit level sets of gauge functions that can be expressed
using Bézier splines comprised of 3 quadratic Bézier curves. The red points
are the end points of each curve, while the blue points are intermediate
points controlling the shapes of the curves.

Because they are parsimoniously parameterized and straightforward to
constrain, Bézier splines are convenient for modeling the boundary $\partial
G$ of the limit set $G$. We specify $\partial G$ as a Bézier spline comprised
of three quadratic Bézier curves $g_{B}=\\{B_{1}(t),B_{2}(t),B_{3}(t)\\}$,
where $B_{1}(t):=B(t;\mathbf{p}_{0},\mathbf{p}_{1},\mathbf{p}_{2})$,
$B_{2}(t):=B(t;\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4})$, and
$B_{3}(t):=B(t;\mathbf{p}_{4},\mathbf{p}_{5},\mathbf{p}_{6})$, for
$\mathbf{p}_{i}\in\mathbb{R}^{2}$, $i=0,1,\ldots,6$. The three curves trace
the paths:

$\displaystyle
B_{1}(t)=(1-t)^{2}\mathbf{p}_{0}+2t(1-t)\mathbf{p}_{1}+t^{2}\mathbf{p}_{2},$
$\displaystyle
B_{2}(t)=(1-t)^{2}\mathbf{p}_{2}+2t(1-t)\mathbf{p}_{3}+t^{2}\mathbf{p}_{4},$
$\displaystyle
B_{3}(t)=(1-t)^{2}\mathbf{p}_{4}+2t(1-t)\mathbf{p}_{5}+t^{2}\mathbf{p}_{6},$

for $0\leq t\leq 1$. We denote the point $\mathbf{p}_{i}:=(p_{i,1},p_{i,2})$,
$0\leq p_{i,1},p_{i,2}\leq 1$, and place two sets of constraints on the curves
in order to elicit valid gauge functions which satisfy conditions 1 and 2. The
first set of constraints ensure that the Bézier spline touches all four edges
of the unit square:

$\displaystyle p_{0,1}$ $\displaystyle=p_{6,2}=0,$ $\displaystyle p_{2,2}$
$\displaystyle=p_{4,1}=1.$

The second set of constraints are sufficient conditions to ensure that the
star-shaped property holds for the spline:

$\displaystyle p_{1,1}$ $\displaystyle\leq p_{2,1},$ $\displaystyle
m(\mathbf{0},\mathbf{p}_{1})$ $\displaystyle\geq
m(\mathbf{0},\mathbf{p}_{2}),$ $\displaystyle m(\mathbf{0},\mathbf{p}_{4})$
$\displaystyle\geq m(\mathbf{0},\mathbf{p}_{5}),$ $\displaystyle p_{4,2}$
$\displaystyle\geq p_{5,2},$ $\displaystyle p_{3,1}$ $\displaystyle=p_{3,2},$
$\displaystyle p_{3,1}$ $\displaystyle\geq\min(p_{2,1},p_{4,2})$

where $m(\mathbf{p},\mathbf{p^{\prime}})$ denotes the slope of the line
connecting the points $\mathbf{p}$ and $\mathbf{p^{\prime}}$, and
$\mathbf{0}=(0,0)$ is the origin. The final condition prevents unrealistic
cases where $p_{2,1}$ and $p_{4,2}$ are both 1, but $p_{3,1}<1$. Thus, we
arrive at a model for the limit set boundary $\partial G$, indexed by the 9
univariate parameters,

$\bm{\theta}_{g}=(p_{0,2},p_{1,1},p_{1,2},p_{2,1},p_{3,1},p_{4,2},p_{5,1},p_{5,2},p_{6,1})^{\text{T}}.$

Figure 3 plots Bézier splines under these constraints, each representing a
gauge function with different dependence properties. Top row plots correspond
to AI scenarios, whereas plots in the bottom row correspond to AD scenarios.
The four red control points are the knots of the spline. The three blue
control points affect the shape, and the spline passes through them only if
they are co-linear with the preceding and proceeding control points. In the
general case, there are 9 coordinates, each admitting a uniform support, which
need to be estimated to fully specify a valid gauge function. Richer models
can be achieved using more control points; however, it would come with
increased computational cost and additional constraints to ensure that
Conditions 1 and 2 hold. The quadratic Bézier spline with 4 knots therefore
constitutes a parsimonious representation for $\partial G$ which is still
flexible enough to capture multiple dependence regimes and mimic most of the
common parametric models.

###  3.2 Statistical inference for the limit set boundary

With a model defined for the limit set boundary $\partial G$, we turn to the
question of how to estimate the shape from iid copies of the random vector
$\bm{X}$ in standard exponential margins. After transforming $\bm{X}$ to
pseudo-polar coordinates $(R,\bm{W})$, a convenient form (Wadsworth and
Campbell, 2024) for the conditional density of a large radius $R$, given the
angle $\bm{W}$, is

$f_{R\,|\,W}(r\,|\,\bm{w})\propto r^{d-1}\exp\\{-rg(\bm{w})[1+o(1)]\\},\quad
r\rightarrow\infty,$

where $d$ is the dimension of $\bm{X}$ (we have only considered $d=2$ here).
For likelihood-based inference, Wadsworth and Campbell (2024) show that the
$o(1)$ term can be moved outside the exponent in most cases and therefore
ignored; they consequently consider the approximation adequate for radii
larger than a threshold $r_{0}(\bm{w})$. This yields a truncated Gamma
likelihood:

$R\,|\,\bm{W}=\bm{w},R>r_{0}(\bm{w}),\bm{\theta}_{g}\sim\text{truncGamma}\left(\alpha,g_{\bm{\theta}_{g}}(\bm{w})\right).$
(1)

For most common bivariate copulas, $\alpha=2$. However, the quality of the
approximation tends to vary (for further details, see Wadsworth and Campbell,
2024), and $\alpha$ is usually treated as a parameter to be estimated. Thus,
given a gauge function $g_{\bm{\theta}_{g}}$ an approximate likelihood for the
large radii given the angles is

$L(\bm{\theta}_{g},\alpha;(r_{1},\bm{w}_{1}),\ldots,(r_{n_{0}},\bm{w}_{n_{0}}))=\prod_{i=1}^{n_{0}}\frac{g_{\bm{\theta}_{g}}(\bm{w})^{\alpha}}{\Gamma(\alpha)}\frac{r_{i}^{\alpha-1}e^{-r_{i}g_{\bm{\theta}_{g}}(\bm{w}_{i})}}{1-F(r_{0}(\bm{w}_{i});\alpha,g_{\bm{\theta}_{g}}(\bm{w}_{i}))},$
(2)

where $n_{0}$ is the number of points exceeding the threshold $r_{0}(w)$, and
$F(\,\cdot\,;\alpha,g_{\bm{\theta}_{g}}(\bm{w}_{i}))$ is the CDF of a Gamma
distribution with shape parameter $\alpha$ and rate parameter
$g_{\bm{\theta}_{g}}(\bm{w}_{i})$.

Figure 4: Schematic of how to calculate $g_{\bm{\theta}_{G}}$ at a data point
$(\bm{x})$, given a boundary curve $\partial G$. The value of the gauge
function is the distance from the origin to $\bm{x}$, relative to the distance
from the origin of the intersection of the ray connecting $\bm{x}$ with the
origin and the boundary $\partial G$.

To calculate the gauge function at each data point, as required in the
likelihood (2), we exploit the homogeneity property of $g$. This gives us that
the value of the gauge function evaluated at a point $\bm{x}$ is the distance
from the origin to $\bm{x}$, relative to the distance from the origin of the
intersection of the ray connecting $\bm{x}$ with the origin and the boundary
$\partial G$. In the schematic in Figure 4, the intersection with $\partial G$
is denoted as $\bm{x}_{\partial G}$, so that

$g_{\bm{\theta}_{G}}(\bm{x})=\frac{\|\bm{x}\|}{\|\bm{x}_{\partial G}\|}.$ (3)

We also need to select a threshold $r_{0}(w)$, as a function of angle.
Wadsworth and Campbell (2024) and Simpson and Tawn (2022) both chose
thresholds as functions of the angle, first as empirical quantiles of moving
windows of angle and then using smooth semi parametric quantile regression. We
employ a much simpler approach, and choose a high quantile in each threshold
marginal component. We have found that this very basic strategy results in
estimation performance at least comparably as good as more complicated
alternatives. In addition, choosing marginal thresholds has two key
advantages. First, it is simple to implement and requires no intricate tuning.
Second, it permits, in principle, transformation to standard exponential
margins within a hierarchical model, whereas thresholds that depend jointly on
both components do not. With this in mind, we choose a value $\tau\in(0,1)$,
and then set each marginal threshold at the $\tau^{\text{th}}$ marginal
empirical quantile $q_{\tau,X_{1}}$ for $X_{1}$ and $q_{\tau,X_{2}}$ for
$X_{2}$. In pseudo-polar coordinates, this gives a radial threshold of

$r_{0}(\bm{w})=\begin{cases}\frac{q_{\tau,X_{2}}}{1-w},&\quad
w\in\left[0,\frac{q_{\tau,X_{1}}}{q_{\tau,X_{1}}+q_{\tau,X_{2}}}\right]\\\
\frac{q_{\tau,X_{1}}}{w},&\quad
w\in\left(\frac{q_{\tau,X_{1}}}{q_{\tau,X_{1}}+q_{\tau,X_{2}}},1\right]\end{cases}$

A sensitivity study comparing the effect of the choice of threshold on
estimation of $\eta$ is provided in the Supplementary Material (Majumder et
al., 2024a, Section A).

###  3.3 Prior distributions for control points

Our model for the limit set boundary $\partial G$ is indexed by the 9
univariate parameters, viz.
$\bm{\theta}_{g}=(p_{0,2},p_{1,1},p_{1,2},\allowbreak
p_{2,1},p_{3,1},p_{4,2},p_{5,1},p_{5,2},p_{6,1})^{\text{T}}$. To inform prior
selection for these control points, we examine the boundary $\partial G$ of
four parametric copula models, and learn the conditions on the control points
which allow the Bézier curve to mimic their shapes. The four copulas that we
consider are the Gaussian, inverted logistic, logistic, and asymmetric
logistic. The first two are AI, while the final two are AD; analytical
expressions of dependence measures for these models are provided in Table 1,
replicated from Simpson and Tawn (2022). Since the limit set boundary
$\partial G$ can take on a variety of shapes, including the AD case where it
touches the upper right-hand corner of the unit box (see e.g., Nolde and
Wadsworth, 2022, Section 4), we let the coordinates of the control points
(i.e., the members of $\bm{\theta}_{g}$) vary freely in $[0,1]$ for
flexibility, which permits the possibility of them being exactly equal to 0 or
1 to accommodate AD as seen in the logistic and asymmetric logistic copulas
and the very weak dependence as seen in the inverted logistic copula. We now
outline the support for the distributions of the control points which will
allow the Bézier splines to mimic our four copulas of interest, and then
specify priors for all parameters to be used in the remainder of this study.

The limit set boundary for an AD copula is obtained whenever $p_{2,1}=1$ or
$p_{4,1}=1$; additionally, the logistic copula is implied by
$\mathbf{p}_{2}=\mathbf{p}_{3}=\mathbf{p}_{4}=(1,1)$ (Figure 3, bottom-right).
This is equivalent to collapsing the second curve of the Bézier spline to a
single point, and can be incorporated into our model by having a semi-
continuous prior distribution for $p_{2,1},p_{3,1},\mbox{ and }p_{4,2}$ with
support over (0,1] which includes a point mass at 1. Similarly, collapsing the
second curve and having a semi-continuous prior on $p_{0,2}$ and $p_{6,1}$ can
incorporate an approximation of the limit set boundary for an asymmetric
logistic copula (Figure 3, bottom left). Finally, approximating the limit set
boundary for an inverted logistic copula using a Bézier spline would require
the first and third curves to collapse onto the $x=0$ and $y=0$ lines
respectively (Figure 3, top-right). To accommodate this case, we set semi-
continuous priors with point masses at 0 for $p_{1,1},p_{2,1},p_{4,2},\mbox{
and }p_{5,2}$. To flexibly accommodate the wide range of limit set boundary
shapes possible within our framework, we set priors as follows:

$\displaystyle\alpha$ $\displaystyle\sim\mbox{LogNormal}(1,1),$ $\displaystyle
p_{1,2},p_{5,1}$
$\displaystyle\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\mbox{Uniform}(0,1).$

The LogNormal prior on $\alpha$ has its density concentrated near $\alpha=2$.
The remaining control points have priors that are the mixture of a standard
Uniform distribution and at least one point mass (to allow important geometric
features of $\partial G$ with positive probability). They have the following
forms:

$\displaystyle p_{1,1}$
$\displaystyle\stackrel{{\scriptstyle\text{iid}}}{{\sim}}0.1\cdot\mathbb{I}(p_{1,1}=0)+0.9\cdot\mbox{Uniform}(0,1),$
$\displaystyle p_{2,1}$
$\displaystyle\stackrel{{\scriptstyle\text{iid}}}{{\sim}}0.1\cdot\mathbb{I}(p_{2,1}=0)+0.8\cdot\mbox{Uniform}(0,1)+0.1\cdot\mathbb{I}(p_{2,1}=1),$
$\displaystyle p_{3,1}\,|\,p_{2,1},\,p_{4,2}$ $\displaystyle\sim
0.6\cdot\mbox{Uniform}(0,1)+0.4\cdot\mathbb{I}(p_{3,1}=1,\max(p_{2,1},p_{4,2})=1).$

The points $p_{5,2}$ and $p_{4,2}$ are distributed identically to $p_{1,1}$
and $p_{2,1}$, respectively. The point mass probabilities were chosen on the
basis of a sensitivity study whose aim was to have good discrimination between
the AD and AI cases based on posterior probabilities of $\eta=1$, while also
simultaneously being able to provide unbiased, consistent estimates of $\eta$
when $\eta<1$. In particular, the prior on $p_{3,1}$ ensures that it is 1 only
if AD is implied by either $p_{2,1}$ or $p_{4,2}$. These prior assumptions can
accommodate the logistic, inverted logistic, and asymmetric logistic copulas,
as well as intermediate forms such as the Gaussian copula which do not require
point masses.

###  3.4 Additional bivariate extremal dependence measures

Alongside the tail dependence coefficient $\eta$ (Ledford and Tawn, 1996)
discussed in Section 2, we consider two additional indices of tail dependence
which can be derived from the gauge function. The first of these is the
angular dependence function $\lambda(\omega)$ which considers different
scalings for the two components of $\bm{X}$. In particular, Wadsworth and Tawn
(2013) considered asymptotic probabilities of the following form:

$P(X_{1}>\omega
x,X_{2}>(1-\omega)x)=\mathcal{L}_{\omega}(e^{x})e^{-x\lambda(\omega)},$

for $\omega\in[0,1]$, $\lambda(\omega)\in(0,1]$, and some function
$\mathcal{L}_{\omega}$ that is slowly varying at infinity. The function
$\lambda(\omega)$ therefore captures both extremal dependence regimes, with AD
corresponding to the pointwise lower bound of $\lambda(\omega)$. Evaluation of
$\lambda(\omega)$ for rays $\omega$ near 0 and 1 corresponds to regions where
one variable is larger than the other. In particular, $\lambda(\omega)$ is a
generalization of $\eta$, with $\eta=1/\\{2\lambda(1/2)\\}$. Murphy-Barltrop
et al. (2023) found that global estimators (such as the Simpson-Tawn
estimator) which simultaneously estimate $\lambda(\omega)$ for all values of
$\omega$ tend to provide better estimates compared to pointwise estimators
(such as the Hill estimator). Since the Bézier spline estimator is also a
global estimator of $\lambda(\omega)$, examining the fit of $\lambda(\omega)$
is a useful measure to compare the two global estimators.

We also investigate the dependence measure $\tau_{1}(\delta)$ (Simpson et al.,
2020), given by:

$P(X_{1}>x,X_{2}\leq\delta
x)=\mathcal{L}_{\delta}(e^{x})e^{-x/\tau_{1}(\delta)},$

for some function $\mathcal{L}_{\delta}$ that is slowly varying at infinity,
with $\delta\in[0,1]$. $\tau_{1}(\delta)$ is monotonically increasing in
$\delta$, with $\tau_{1}(1)=1$. This dependence measure characterizes the
probability of $X_{1}$ being large while $X_{2}$ is of a smaller order.
Specifically, if there exists a $\delta^{*}<1$ such that
$\tau_{1}(\delta^{*})=1$, it implies that $X_{1}$ can be large while $X_{2}$
is small (with $\delta$ determining just how small). If no such $\delta^{*}$
exists, then $X_{1}$ can be large only if $X_{2}$ is also large. We can define
$\tau_{2}(\delta)$ analogously, and both
$\tau_{1}(\delta),\tau_{2}(\delta)\in(0,1]$.

Table 1 provides the analytical expressions for these measures for the four
dependence copulas that we consider in our study. Like with $\eta$, these
dependence measures can be exactly deduced from limit set boundaries, and
hence from gauge functions. Since the Bézier splines are quadratic
polynomials, we can easily calculate $\tau_{1}(\delta)$ and $\lambda(\omega)$
for any estimated limit set boundary, simply by finding the intersections of
polynomials and lines, which have closed-form solutions. While we do not
present results for $\tau_{2}(\delta)$ in our study, it can be calculated in
the same manner as $\tau_{1}(\delta)$. We refer the reader to Simpson and Tawn
(2022) for a detailed discussion on how each of these measures can be obtained
from gauge functions in more general settings.

Copula | $g(\bm{x})=g(x_{1},x_{2})$
---|---
Gaussian | $\\{x_{1}+x_{2}-2\rho(x_{1}x_{2})^{1/2}\\}/(1-\rho^{2})$
Logistic | $\gamma^{-1}\max(x_{1},x_{2})+(1-\gamma^{-1})\min(x_{1},x_{2})$
Inv-Logistic | $(x_{1}^{1/\gamma}+x_{2}^{1/\gamma})^{\gamma}$
Asy-Logistic | $\begin{cases}\gamma^{-1}\max(x_{1},x_{2})+(1-\gamma^{-1})\min(x_{1},x_{2}),&\mbox{ if }\gamma<1\\\ x_{1}+x_{2},&\mbox{ if }\gamma=1\end{cases}$
(a) Gauge function $g$ for the four bivariate copulas.

Copula | $\eta$ | $\lambda(\omega)$ | $\tau_{1}(\delta)=\tau_{2}(\delta)$
---|---|---|---
Gaussian | $(1+\rho)/2$ | $\begin{cases}\max(\omega,1-\omega),&\mbox{ if }t_{\omega}\leq\rho^{2}\\\ \frac{1-2\rho\sqrt{\omega(1-\omega)}}{1-\rho^{2}},&\mbox{ if }t_{\omega}\geq\rho^{2}\end{cases}$ | $\begin{cases}1,&\mbox{ if }\delta\geq\rho^{2}\\\ \frac{1-\rho^{2}}{1+\delta-2\rho\sqrt{\delta}},&\mbox{ if }\delta\leq\rho^{2}\end{cases}$
Logistic | 1 | $\max(\omega,1-\omega)$ | $(\gamma^{-1}+1-\gamma^{-1}\delta)^{-1}$
Inv-Logistic | $2^{-\gamma}$ | $\\{\omega^{1/\gamma}+(1-\omega)^{1/\gamma}\\}^{\gamma}$ | 1
Asy-Logistic | 1 | $\max(\omega,1-\omega)$ | 1
(b) Dependence measures for the four bivariate copulas. Here,
$t_{\omega}=\min(\omega,1-\omega)/\max(\omega,1-\omega)$.

Table 1: Gauge function $g$, and a summary of dependence measures for the four
bivariate copulas used in our study. Table has been reproduced from Simpson
and Tawn (2022).

##  4 Simulation Study

###  4.1 Study Setup

We demonstrate the appropriateness of using Bézier splines to model the limit
set boundary corresponding to the gauge function by means of a simulation
study comparing its performance with the Simpson-Tawn estimator (Simpson and
Tawn, 2022). We consider four bivariate copulas to generate dependent data
with exponential marginal distributions: the Gaussian, the logistic, the
inverted logistic, and the asymmetric logistic. The Gaussian copula is
parameterized by its correlation $\rho\in[0,1)$, while the dependence
parameter for the remaining three copulas is $\gamma\in(0,1)$. Table 1 lists
the gauge functions associated with each copula, as well as the set of
corresponding extremal dependence coefficients
$\\{\eta,\lambda(\omega),\tau_{1}(\delta)\\}$. The four copulas cover a range
of AD behavior; the logistic and asymmetric logistic copulas, in particular,
are AD. Further details for these four copulas can be found in Nolde and
Wadsworth (2022). For each copula, we consider five parameter settings:
$\rho,\gamma=\\{0.3,0.4,0.5,0.6,0.7\\}$. Smaller values of $\gamma$ correspond
to stronger tail dependence, whereas larger values of $\rho$ lead to stronger
tail dependence.

For each copula and parameter combination, we generate 100 datasets of
$n=5,000$ data points each. The datasets are converted to pseudo-polar
coordinates $(R,\bm{W})$, and the $n_{0}$ points that are above the
$\tau=0.75$ quantile marginal threshold for at least one variable are used to
model the gauge function for each dataset. The $n_{0}$ radii are assumed to
approximately distributed according to a truncated Gamma distribution with a
common shape parameter $\alpha$ and a rate parameter equal to an appropriate
gauge function evaluated at the data point. We use Metropolis updates for all
parameters and run 11,000 MCMC iterations for each dataset, discarding the
first 1,000 as burn-in. All Metropolis updates are tuned to give an acceptance
probability of 0.4, and posterior convergence is diagnosed based on the visual
inspection of trace plots. We compare our estimated limit set boundaries with
the Simpson-Tawn estimator in terms of how well they estimate the set of
dependence coefficients outlined in Table 1. The Simpson-Tawn estimator for
the datasets was evaluated using the default settings recommended by the
authors in Simpson and Tawn (2022). We use the root mean square error (RMSE)
to compare estimates of the scalar $\eta$, and the root mean integrated square
error (RMISE) to compare estimates of the functions $\lambda(\omega)$ and
$\tau_{1}(\delta)$. The methodology is implemented in R (R Core Team, 2023);
the code is available on GitHub through the BezELS package (Majumder et al.,
2024b).

###  4.2 Parameter estimates

(a) Estimated gauge functions in pseudo-polar space.

(b) Estimated gauge functions in Euclidean space.

Figure 5: Limit set boundaries based on Bézier splines (blue) and
corresponding functions for the data-generating model (red) in pseudo-polar
and Euclidean space for the Gaussian, logistic, inverted logistic, and
asymmetric logistic copulas with dependence parameters set to $0.5$

.

Each panel in Figure 5(b) plots the limit set boundaries elicited by the
estimated Bézier splines based on the posterior distribution from a single
dataset with the dependence parameter set to 0.5. Plots in Figure 5(a) display
the dependence modeling in pseudo-polar coordinates. The dashed grey line
corresponds to the threshold $r_{0}(\bm{w})$ for angles $\bm{w}$. Each
estimated limit set boundary is represented as a functional boxplot (Hyndman
and Shang, 2010; Sun and Genton, 2011), a visual representation of functional
data analogous to a classical boxplot. Each functional quantile depicted in
the functional boxplot is a function contained in the sample; in this case,
the sample consists of limit set boundaries based on Bézier splines, in
pseudo-polar coordinates, drawn from the posterior. The curves are ordered
according to a notion of modified band depth (López-Pintado and Romo, 2009).
The median is shown in dark blue, and the limit set boundary corresponding to
the data-generating model is shown in red. The envelope represents the $50\%$
central region, analogous to the box of a classical boxplot. The outer light
blue lines of the functional boxplot correspond to the whiskers of a classical
boxplot. Finally, the vertical lines indicate the maximum envelope of the
functions except outliers. Plots in Figure 5(b) display the dependence models
in Euclidean coordinates, with the median Bézier spline in dark blue and the
limit set boundary corresponding to the data-generating model in red. They are
overlaid on Bézier splines evaluated from 500 draws from the posterior
distribution, plotted in gray. The boxplots on the top and right margins
correspond to the posterior distributions of $p_{2,1}$ and $p_{4,2}$
respectively, which serve as a visual indicator of asymptotic dependence in
the data. Specifically, if the median of either boxplot is 1, the posterior of
$\eta$ is 1. The Bézier splines are able to adequately represent the geometric
form of all four dependence copulas.

Figure 6 shows boxplots of posterior median of $\eta$ for the four dependence
copulas based on the Bézier spline (blue) and Simpson-Tawn estimators (green).
Analytical values of $\eta$ obtained using expressions in Table 1 are shown as
red dots in each plot, and the coverage of equi-tailed 95% intervals are noted
in plain text below each boxplot. Coverage for the Simpson-Tawn estimator is
based on 100 bootstrapped samples from each of the 100 datasets. We plot the
median instead of the mean for the Bézier spline estimates since the posterior
distributions are often highly asymmetric due to point-mass prior
distributions. The Bézier spline estimator has low bias and nominal coverage
for estimates of $\eta$. The Simpson-Tawn estimator has nominal or near-
nominal coverage in most cases except for the asymmetric logistic copula when
the dependence parameter $\gamma$ is 0.5 or higher. It also has noticeably
higher bias than the Bézier spline estimator for both AD copulas, and shows a
sharp decline in coverage as the strength of dependence drops for the
asymmetric logistic copula.

Figure 6: Sampling distribution of the posterior medians of $\eta$ based on
the Bézier spline estimate (blue) for the four dependence copulas, alongside
estimates using the Simpson-Tawn estimator (green). The red dots indicate the
true values, and coverage of equi-tailed 95% intervals are noted below each
boxplot.

(a) Sampling distribution of the posterior medians of $\lambda(0.40)$.

(b) Sampling distribution of the posterior medians of $\tau_{1}(0.25)$.

Figure 7: Sampling distribution of the posterior medians of $\lambda(0.40)$
and $\tau_{1}(0.25)$ based on the Bézier spline (blue) and Simpson-Tawn
(green) estimators for four dependence copulas, with dependence parameters set
to 0.5. The red lines indicate the true values, and coverage of equi-tailed
95% intervals are noted below each boxplot.

We evaluate $\lambda(\omega)$ and $\tau_{1}(\delta)$ based on the Bézier
spline and Simpson-Tawn estimators for
$\omega,\delta=0.01,0.02,\allowbreak\ldots,0.99$. Figure 7(a) shows boxplots
for the posterior medians of $\lambda(0.40)$ based on the two estimators, with
the corresponding dependence parameter set to 0.5. Figure 7(b) similarly plots
the distribution of $\tau_{1}(0.25)$ for the two estimators. Both estimators
are better at estimating $\lambda(\omega)$ than $\tau_{1}(\delta)$, and the
Bézier spline estimator tends to have better coverage than the Simpson-Tawn
estimator.

Table 2 summarizes the RMSE ratio for estimates of $\eta$ and RMISE ratios for
estimates of $\lambda(\omega)$ and $\tau_{1}(\delta)$ based on the Bézier
spline and Simpson-Tawn estimators. Most of the values are greater than 1,
indicating that dependence measures based on Bézier spline estimates of the
gauge function have comparable or better RMSE/RMISE than those based on the
Simpson-Tawn estimator. However, for the inverted logistic copula, the
Simpson-Tawn estimator outperforms the Bézier spline estimator when estimating
of $\tau_{1}(\delta)$. Our experiments indicate that the error arises when the
posterior $p_{4,2}>\delta$; this leads to $\tau_{1}(\delta)$ estimates to be
less than 1, whereas the theoretical value for the inverted logistic copula is
1 for all $\delta$. Our approach, however, is still able to correctly estimate
$\tau(\delta)=1$ in all the inverted logistic scenarios for all but extremely
small values of $\delta$.

Table 2: RMSE ratios for estimates of $\eta$ and RMISE ratios for estimates of $\lambda(\omega)$ and $\tau_{1}(\delta)$ based on the Bézier spline $(\hat{\eta},\hat{\lambda},\mbox{ and }\hat{\tau})$ and Simpson-Tawn $(\tilde{\eta},\tilde{\lambda},\mbox{ and }\tilde{\tau})$ estimators over simulated datasets for four copulas and five dependence levels. Measure | Copula | Dependence parameter value
---|---|---
|  | 0.3 | 0.4 | 0.5 | 0.6 | 0.7
$\frac{RMSE(\tilde{\eta})}{RMSE(\hat{\eta})}$ | Gaussian | 1.02 | 1.11 | 1.04 | 1.06 | 0.81
Logistic | 1.38 | 0.81 | 1.36 | 1.17 | 1.07
Inv-Logistic | 0.70 | 0.99 | 1.11 | 1.03 | 1.01
Asy-Logistic | 2.39 | 2.01 | 2.23 | 2.37 | 2.39
$\frac{RMISE(\tilde{\lambda})}{RMISE(\hat{\lambda})}$ | Gaussian | 1.19 | 1.18 | 1.05 | 0.94 | 0.75
Logistic | 1.35 | 0.61 | 1.10 | 1.04 | 1.28
Inv-Logistic | 0.82 | 1.10 | 1.20 | 1.10 | 1.04
Asy-Logistic | 6.17 | 6.65 | 10.91 | 7.43 | 5.67
$\frac{RMISE(\tilde{\tau})}{RMISE(\hat{\tau})}$ | Gaussian | 0.85 | 1.00 | 1.08 | 1.17 | 1.18
Logistic | 1.00 | 0.95 | 0.89 | 0.83 | 1.00
Inv-Logistic | 0.23 | 0.27 | 0.17 | 0.28 | 0.17
Asy-Logistic | 2.38 | 2.96 | 1.77 | 1.72 | 1.23

Table 3 provides the number of datasets (out of 100) in each scenario where
the posterior medians of $\eta$ is estimated to be 1. The values in the
parentheses are corresponding point estimates of $\eta$ based on the Simpson-
Tawn estimator. For the Bézier spline estimators, the values were always
estimated to be near 0 for the AI copulas, and always high for the AD copulas.
While there were some cases where the Bézier spline estimates the posterior of
$\eta$ as 1 when the dependence is high in an AI copula, both methods are good
at estimating AI behaviour correctly. On the other hand, the Simpson-Tawn
estimates show a decline in their ability to estimate the correct value of
$\eta=1$ for the AD copulas when dependence is low. This is especially
noticeable for the asymmetric logistic copula, and has been documented by
Simpson and Tawn (2022) as well. The Bézier spline estimator is much better at
predicting AD correctly across all scenarios. We conclude that the Bézier
splines are adept at representing limit set boundaries associated with common
parametric copula models, and are also flexible enough to represent a wider
variety of edge cases. In all cases, the true value of $\eta$ was well-
estimated from the posterior distribution, and the model is particularly adept
at identifying AD.

Table 3: Number of datasets (out of 100) where the posterior median of $\eta$ is 1 for each scenario. Values in parenthesis correspond to the Simpson-Tawn estimator. | Dependence parameter value
---|---
| 0.3 | 0.4 | 0.5 | 0.6 | 0.7
Gaussian | 00 (00) | 00 (00) | 00 (00) | 00 (00) | 06 (00)
Logistic | 98 (91) | 93 (79) | 96 (77) | 91 (74) | 82 (63)
Inv-Logistic | 02 (00) | 00 (00) | 00 (00) | 00 (00) | 00 (00)
Asy-Logistic | 96 (30) | 80 (25) | 82 (22) | 83 (14) | 84 (00)

###  4.3 Additional simulation studies

Two additional simulation studies are presented in the Supplementary Material
(Majumder et al., 2024a). In both cases, the results are compared based on
RMSE/RMISE values of bivariate extremal dependence coefficients, as well as
how often it estimates $\eta=1$. The first study considers how our relatively
straightforward method of selecting the quantile threshold affects the
estimation process. This is carried out by comparing it against an ‘oracle’
threshold which is an asymptotic approximation to the true conditional
quantile $q_{\tau}(w)$ and requires knowledge of the true gauge function to
compute. We are unable to find any meaningful improvement in the estimates
when we consider an oracle threshold, which indicates that our choice of
threshold is adequate for the scenarios we have considered. The second study
repeats the simulation study presented in this section, but for a sample size
$n=600$. This is carried out to ensure that our approach is still valid for
small data sizes like the one that arises in our second application, presented
in the following section. Our results indicate that despite slightly higher
RMISE values and slightly lower coverage, the Bézier splines can still capture
the shape of the true limit set boundary with low bias. Comparisons with the
Simpson-Tawn estimator provide results that are quite similar to the ones
presented in this section.

##  5 Applications

###  5.1 Analysis of the Santa Ana Winds data

Figure 8: Santa Ana wind speeds and dryness measured at the March Air Force
Base station. Data above the $0.75$ marginal quantile threshold are in red.

We apply our method to the Santa Ana winds and dryness data (Cooley et al.,
2019). The Santa Ana winds are a multivariate meteorological regime that has
been implicated as a major driver of large wildfires in southern California
(Billmire et al., 2014). Wildfires are related to several conditions like
temperature, humidity, wind speed, and fuel supply (Littell et al., 2018).
Historically, the autumn months of September, October, and November have had a
higher number of wildfires compared to the winter months, and are associated
with warm temperatures, low humidity, and high winds. Cooley et al. (2019)
surmised that the data exhibits AD and used the framework of regular variation
to estimate probabilities associated with two different risk regions. The
regular variation structure employed by them, however, cannot capture the
nuance of AI and models all levels of asymptotic independence (including exact
independence) within a single degenerate model.

We consider daily dryness (%) and wind speed (m/s) data collected at the March
Air Reserve Base station in Riverside County from the HadISD dataset (Dunn et
al., 2012). The dryness is defined in this case as the $100-$RH, where RH is
the relative humidity measured as a percentage. The bivariate time series
represents a measure of the daily risk of fire. The station’s variables have
appeared to be associated with known Santa Ana events. The data consists of
3,902 days for the months of September–November from 1972–2015; we assume
temporal stationarity as in Cooley et al. (2019). Figure 8 plots the data both
in its original scale as well as the rank transformed scale. The data shows
tail dependence, noticeable in the rank transformed data with a cluster of
values in the upper right corner. Our goal is to study the tail dependence
between the two variables by estimating a gauge function for the data after a
further transformation to unit exponential margins.

We analyze this data at a threshold of $\tau=0.75$, providing us with 1,529
points that are above the threshold in at least one margin. The conditional
distribution of the radii are assumed to be a truncated Gamma distribution for
this analysis. We run 2 MCMC chains of 11,000 iterations each, discarding the
first 1,000 from each chain as burn-in. The priors and the remainder of the
MCMC settings are identical to the simulation study.

Figure 9: Limit set boundaries based on Bézier splines in pseudo-polar space
(left) and Euclidean space (right) for assessing tail dependence between Santa
Ana windspeed ($X_{1}$) and dryness ($X_{2}$). Median curves are plotted in
dark blue.

Figure 10: PP plot (left) and exponential QQ plot (right) for the truncated
Gamma model fitted to the Santa Ana winds data. The black points correspond to
the fit for the median Bézier spline. Gray lines are based on 100 random draws
from the posterior.

Figure 9 plots the estimated limit set boundaries, with the median curve
plotted in blue. The plot on the left is in pseudo-polar coordinates and
depicts a functional boxplot of the estimated limit set boundary, while the
plot on the right is in Euclidean coordinates. The posterior median of $\eta$
is estimated to be 1, and $\mbox{P}(\eta=1\,|\,\bm{X})=0.60$, suggesting AD
between wind speed and dryness. Figure 10 evaluates the goodness-of-fit for
the truncated Gamma model in terms of PP plots as well as QQ plots, where each
point is transformed to unit exponential using the shape and rate parameters
of the assumed truncated Gamma distribution. In both cases, the black points
are obtained from Gamma parameters implied by the median curve of the limit
set boundary, while the gray lines correspond to 100 random draws from the
posterior. Both plots suggest that the model fits the data adequately.

To test sensitivity to the threshold level, we repeated the experiment at an
additional threshold levels of $\tau=0.70,0.80,\mbox{ and }0.90$. In all 3
additional cases, $\mbox{P}(\eta=1\,|\,\bm{X})>0.50$, with the smallest value
(of 0.56) occurring for $\tau=0.90$ and the largest value (of 0.64) for
$\tau=0.70$. Taken together, we conclude that dryness and windspeed for the
Santa Ana dataset is asymptotically dependent with high probability,
consistent with the results of Cooley et al. (2019).

###  5.2 Analysis of Ozone concentration data

Figure 11: PP plot (left) and exponential QQ plot (right) for the truncated
Gamma model fitted to the Ozone concentration data. The gray lines correspond
to the fit for the median Bézier splines at 100 random locations.

For our second study, we consider air pollution measurement data across the US
from the Community Multiscale Air Quality (CMAQ) model (Binkowski and Roselle,
2003; Wyat Appel et al., 2007, 2008) as well as EPA Air Quality System (AQS)
(US EPA, 2017) data for the contiguous US. While CMAQ is a numerical model
available across the entire country at a 12km resolution, the AQS dataset
consists of observations monitored at $1,376$ stations across the US. Among
them, only 519 stations had over 600 observations, which is what we use for
this analysis. The full dataset has previously been used by Gong et al. (2021)
to develop a combined data product for 12 air pollutants. When fusing data
products, it is important to calibrate the model data to ground truth. In our
application, we will verify how strong the dependence is between the AQS and
CMAQ datasets for ozone, one of the 12 pollutants made available by both
datasets.

Our data consists of daily ozone readings for the months of July–September
from 2010–2014, resulting in a bivariate time series of CMAQ and AQS data for
up to 610 days at each station. The sample correlations between the AQS and
CMAQ data for the 519 stations range from 0.29–0.86 with a median of 0.69,
suggesting a high level of agreement in the bulk of the distribution. To
assess tail dependence, we fit a gauge function with truncated Gamma
likelihood for data censored at $\tau=0.75$ threshold independently at every
station. We run 2 MCMC chains for 11,000 iterations each for each station’s
data, discarding the first 1,000 as burn-in.

Figure 12: Posterior median of $\eta$ at 519 AQS monitoring stations in the
US.

The posterior median of $\eta$ has an average value of 0.81 across the 519
locations, and is 1 (AD) for 79 of those stations. This suggests that the CMAQ
data product can adequately represent the tail behavior of observational
ambient ozone. Figure 12 plots the posterior probability of AD based on the
truncated Gamma model. While we are unable to discern any spatial pattern for
high or low posterior values of $\eta$ from the map, we do note that several
of the low values are in urban areas with high population densities.

Finally, To study the sensitivity of our posterior to $\tau$, we repeated our
analysis at two additional threshold levels of $\tau=0.70\mbox{ and }0.80$.
Estimates from both these cases were quite similar to the baseline case of
$\tau=0.75$, with correlations $>0.90$ for both the posterior median of $\eta$
and $\mbox{P}(\eta=1)$. There were 76 and 80 locations respectively where the
posterior median of $\eta$ was 1, and 62 of those locations were shared with
the baseline case. Thus, our results are not very sensitive to the choice of
threshold for small data sizes.

##  6 Discussion

Key aspects of tail dependence in multivariate distributions can be described
through their corresponding gauge functions. In this study, we propose a semi-
parametric method for estimating gauge functions by modeling their unit level
sets as Bézier splines comprised of three quadratic Bézier curves. The splines
can represent the gauge function unit level sets, and hence limit set
boundaries, of varying shapes, and are parsimoniously parameterized by a small
number of control points. The quadratic specification makes it straightforward
to obtain analytical solutions for the shape of the limit set, and constraints
on the control points ensure that the resultant shapes are valid limit set
boundaries. Bayesian estimation of the Bézier splines requires only standard
MCMC techniques and allows important cases on the edge of the parameter space
to be represented by employing mixture priors with point masses. We
demonstrate the efficacy of our model using numerical studies as well as two
real data applications involving fire weather in California, and ambient air
pollution from ozone across the US.

We have only considered bivariate random vectors here, but the modeling
strategy can scale to three dimensions by using Bézier surfaces, with the
control points of constituent Bézier curves set in $\mathbb{R}^{3}$ instead of
$\mathbb{R}^{2}$. It will, however, require more complicated constraints to
ensure that the star shaped property holds, and dimensions greater than three
appear to be infeasible. In addition, it appears possible to extend our
modeling framework to include negative dependence by transforming to Laplace
margins rather than exponential margins. This has been previously suggested
(Simpson and Tawn, 2022; Wadsworth and Campbell, 2024) and recently
implemented in the context of radially stable Pareto distributions
(Papastathopoulos et al., 2023), and within a semi-parametric angular-radial
model (Mackay and Jonathan, 2023). Implementing our Bézier model in Laplace
margins (in two dimensions) would require somewhere between two to four times
the number of control points as we have now, specifying appropriate
constraints on their support, and more sophisticated sampling algorithms to
ensure convergence. Finally, while three Bézier curves are sufficient to
ensure that the boundary of our estimated limit set in the two dimensional
case exactly touches the corners of the unit box, the curves themselves don’t
necessarily need to be quadratic. It is possible to use a nonparametric
Bayesian framework to construct curves that have an arbitrary number of
intermediate control points (i.e., control points excluding the start and end
points). Though this would be computationally more expensive, the resulting
Bézier spline is likely to have broader support over the space of all limit
set boundaries.

## Funding

This work was supported by grants from the Southeast National Synthesis
Wildfire and the United States Geological Survey’s National Climate Adaptation
Science Center (G21AC10045), and the National Science Foundation (DMS-2001433,
DMS-2152887).

## References

  * Balkema et al. (2010) Balkema, A. A., Embrechts, P., and Nolde, N. (2010), “Meta densities and the shape of their sample clouds,” J. Multivariate Anal., 101, 1738–1754.
  * Balkema and Nolde (2010) Balkema, G. and Nolde, N. (2010), “Asymptotic independence for unimodal densities,” Adv. in Appl. Probab., 42, 411–432.
  * Billmire et al. (2014) Billmire, M., French, N. H. F., Loboda, T., Owen, R. C., and Tyner, M. (2014), “Santa Ana winds and predictors of wildfire progression in southern California,” Int. J. Wildland Fire, 23, 1119–1129.
  * Binkowski and Roselle (2003) Binkowski, F. S. and Roselle, S. J. (2003), “Models-3 Community Multiscale Air Quality (CMAQ) model aerosol component 1. Model description,” J. Geophys. Res. - Atmos., 108.
  * Boldi and Davison (2007) Boldi, M.-O. and Davison, A. C. (2007), “A mixture model for multivariate extremes,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69, 217–229.
  * Cooley et al. (2019) Cooley, D., Hunter, B. D., and Smith, R. L. (2019), “Univariate and multivariate extremes for the environmental sciences,” in Handbook of environmental and ecological statistics, Chapman and Hall/CRC, pp. 153–180.
  * Davis et al. (1988) Davis, R. A., Mulrow, E., and Resnick, S. I. (1988), “Almost sure limit sets of random samples in ${\bf R}^{d}$,” Adv. in Appl. Probab., 20, 573–599.
  * de Carvalho et al. (2022) de Carvalho, M., Kumukova, A., and dos Reis, G. (2022), “Regression-type analysis for multivariate extreme values,” Extremes, 25, 595–622.
  * Dunn et al. (2012) Dunn, R. J. H., Willett, K. M., Thorne, P. W., Woolley, E. V., Durre, I., Dai, A., Parker, D. E., and Vose, R. S. (2012), “HadISD: a quality-controlled global synoptic report database for selected variables at long-term stations from 1973–2011,” Clim. Past, 8, 1649–1679.
  * Farouki (2012) Farouki, R. T. (2012), “The Bernstein polynomial basis: a centennial retrospective,” Comput. Aided Geom. D., 29, 379–419.
  * Gong et al. (2021) Gong, W., Reich, B. J., and Chang, H. H. (2021), “Multivariate spatial prediction of air pollutant concentrations with INLA,” Environ. Res. Commun., 3, 101002.
  * Hanson et al. (2017) Hanson, T. E., de Carvalho, M., and Chen, Y. (2017), “Bernstein polynomial angular densities of multivariate extreme value distributions,” Statistics & Probability Letters, 128, 60–66.
  * Hazewinkel (2012) Hazewinkel, M. (2012), Encyclopaedia of Mathematics, no. Vol. 1 in Encyclopaedia of Mathematics, Springer Netherlands.
  * Heffernan and Tawn (2004) Heffernan, J. E. and Tawn, J. A. (2004), “A conditional approach for multivariate extreme values,” J. R. Stat. Soc. Ser. B Stat. Methodol., 66, 497–546, with discussions and reply by the authors.
  * Hill (1975) Hill, B. M. (1975), “A Simple General Approach to Inference About the Tail of a Distribution,” Ann. Stat., 3, 1163–1174.
  * Hyndman and Shang (2010) Hyndman, R. J. and Shang, H. L. (2010), “Rainbow Plots, Bagplots, and Boxplots for Functional Data,” J. Computat. Graph. Stat., 19, 29–45.
  * Kinoshita and Resnick (1991) Kinoshita, K. and Resnick, S. I. (1991), “Convergence of scaled random samples in ${\bf R}^{d}$,” Ann. Probab., 19, 1640–1663.
  * Ledford and Tawn (1996) Ledford, A. W. and Tawn, J. A. (1996), “Statistics for near independence in multivariate extreme values,” Biometrika, 83, 169–187.
  * Littell et al. (2018) Littell, J. S., McKenzie, D., Wan, H. Y., and Cushman, S. A. (2018), “Climate Change and Future Wildfire in the Western United States: An Ecological Approach to Nonstationarity,” Earth’s Future, 6, 1097–1111.
  * López-Pintado and Romo (2009) López-Pintado, S. and Romo, J. (2009), “On the Concept of Depth for Functional Data,” J. Am. Stat. Assoc., 104, 718–734.
  * Mackay and Jonathan (2023) Mackay, E. and Jonathan, P. (2023), “Modelling multivariate extremes through angular-radial decomposition of the density function,” arXiv preprint, arXiv:2310.12711.
  * Majumder et al. (2024a) Majumder, R., Reich, B. J., Shaby, B. A., and Cooley, D. S. (2024a), “Supplement to “Semiparametric Estimation of the Shape of the Limiting Multivariate Point Cloud”,” .
  * Majumder et al. (2024b) Majumder, R., Shaby, B. A., Reich, B. J., and Cooley, D. S. (2024b), BezELS: Bezier splines for Estimating Limit Sets, r package version 0.1.0.
  * Marcon et al. (2017) Marcon, G., Padoan, S., Naveau, P., Muliere, P., and Segers, J. (2017), “Multivariate nonparametric estimation of the Pickands dependence function using Bernstein polynomials,” J. Stat. Plan. Infer., 183, 1–17.
  * Marcon et al. (2014) Marcon, G., Padoan, S. A., Naveau, P., and Muliere, P. (2014), “Nonparametric estimation of the dependence among multivariate rainfall maxima,” in METMA VII-GRASPA 14.
  * Murphy-Barltrop et al. (2023) Murphy-Barltrop, C., Wadsworth, J., and Eastoe, E. (2023), “Improving estimation for asymptotically independent bivariate extremes via global estimators for the angular dependence function,” arXiv preprint, arXiv:2303.13237.
  * Nolde (2014) Nolde, N. (2014), “Geometric interpretation of the residual dependence coefficient,” J. Multivariate Anal., 123, 85–95.
  * Nolde and Wadsworth (2022) Nolde, N. and Wadsworth, J. L. (2022), “Linking representations for multivariate extremes via a limit set,” Adv. Appl. Probab., 54, 688–717.
  * Padoan and Rizzelli (2022) Padoan, S. A. and Rizzelli, S. (2022), “Consistency of Bayesian inference for multivariate max-stable distributions,” The Annals of Statistics, 50, 1490–1518.
  * Papastathopoulos et al. (2023) Papastathopoulos, I., De Monte, L., Campbell, R., and Rue, H. (2023), “Statistical inference for radially-stable generalized Pareto distributions and return level-sets in geometric extremes,” arXiv preprint, arXiv:2310.06130.
  * R Core Team (2023) R Core Team (2023), R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.
  * Sabourin and Naveau (2014) Sabourin, A. and Naveau, P. (2014), “Bayesian Dirichlet mixture model for multivariate extremes: a re-parametrization,” Computational Statistics & Data Analysis, 71, 542–567.
  * Sabourin et al. (2013) Sabourin, A., Naveau, P., and Fougères, A.-L. (2013), “Bayesian model averaging for multivariate extremes,” Extremes, 16, 325–350.
  * Simpson and Tawn (2022) Simpson, E. S. and Tawn, J. A. (2022), “Estimating the limiting shape of bivariate scaled sample clouds for self-consistent inference of extremal dependence properties,” arXiv preprint, arXiv:2207.02626.
  * Simpson et al. (2020) Simpson, E. S., Wadsworth, J. L., and Tawn, J. A. (2020), “Determining the dependence structure of multivariate extremes,” Biometrika, 107, 513–532.
  * Sun and Genton (2011) Sun, Y. and Genton, M. G. (2011), “Functional Boxplots,” J. Computat. Graph. Stat., 20, 316–334.
  * US EPA (2017) US EPA (2017), “Air Quality System Database,” Accessed on 3 August, 2017\.
  * Vettori et al. (2018) Vettori, S., Huser, R., and Genton, M. G. (2018), “A comparison of dependence function estimators in multivariate extremes,” Stat. Comput., 28, 525–538.
  * Wadsworth and Campbell (2024) Wadsworth, J. L. and Campbell, R. (2024), “Statistical inference for multivariate extremes via a geometric approach,” Journal of the Royal Statistical Society Series B: Statistical Methodology, qkae030.
  * Wadsworth and Tawn (2013) Wadsworth, J. L. and Tawn, J. A. (2013), “A new representation for multivariate tail probabilities,” Bernoulli, 19, 2689–2714.
  * Wyat Appel et al. (2008) Wyat Appel, K., Bhave, P. V., Gilliland, A. B., Sarwar, G., and Roselle, S. J. (2008), “Evaluation of the community multiscale air quality (CMAQ) model version 4.5: Sensitivities impacting model performance; Part II—particulate matter,” Atmos. Environ., 42, 6057–6066.
  * Wyat Appel et al. (2007) Wyat Appel, K., Gilliland, A. B., Sarwar, G., and Gilliam, R. C. (2007), “Evaluation of the Community Multiscale Air Quality (CMAQ) model version 4.5: Sensitivities impacting model performance: Part I—Ozone,” Atmos. Environ., 41, 9603–9615.