Inference for bivariate extremes via a semi-parametric angular-radial model

The modelling of multivariate extremes is an active area of research, with applications spanning many domains, including meteorology [6], metocean engineering [24, 62], financial modelling [3] and flood risk assessment [13]. Typically, approaches in this research field are comprised of two steps: first, modelling the extremes of individual variables and transforming to common margins, followed by modelling of the dependence between the extremes of different variables. We refer to this dependence as the extremal dependence structure henceforth.

This article discusses inference for multivariate extremes using an angular-radial model for the probability density function, illustrated using examples in two dimensions. To place the proposed model in context, we first provide a brief synopsis of the existing literature for multivariate extremes. Given a random vector $(X,Y)\in\mathbb{R}^{2}$ with marginal distributions functions $F_{X}$ and $F_{Y}$, the strength of dependence in the upper tail of $(X,Y)$ can be quantified in terms of the tail dependence coefficient, $\chi\in[0,1]$, defined as the limiting probability

when this limit exists [23]. When $\chi>0$, the components of $(X,Y)$ are said to be asymptotically dependent (AD) in the upper tail, and when $\chi=0$, they are said to be asymptotically independent (AI). Much of the focus of recent work in multivariate extreme value theory has been related to developing a general framework for modelling joint extremes of $(X,Y)$ which is applicable to both AD and AI cases, and can be used to evaluate joint tail behaviour in the region where at least one variable is large.

To discuss the approaches proposed to date and their associated limitations, it is helpful to categorise them in terms of whether they assume heavy- or light-tailed margins, and whether they consider the distribution or density function. Classical multivariate extreme value theory assumes heavy-tailed margins, and is based on the framework of multivariate regular variation [MRV, 53]. It addresses the case where $\chi>0$, and has been widely studied – see Beirlant et al., [2], de Haan and Ferreira, [11] and Resnick, [54] for reviews. Under some regularity conditions, equivalent asymptotic descriptions of joint extremal behaviour can be obtained from either the density or distribution function [12].

In the MRV framework, any distribution with $\chi=0$ has the same asymptotic representation. To address this issue, [30, 31] proposed a method to characterise joint extremes for both AI and AD distributions in the region where both variables are large. Model forms for the Ledford-Tawn representation were proposed by Ramos and Ledford, [50]. The resulting framework also assumes heavy-tailed margins and is referred to as hidden regular variation [HRV, 52]. However, for AI distributions, a description of extremal behaviour in the region where both variables are large may not be the most useful, since extremes of both variables are unlikely to occur simultaneously. Moreover, for AI distributions with certain regularity conditions, the asymptotic representation in this framework is governed only by the properties of the distribution along the line $y=x$ [36]. To provide a more useful representation for AI distributions, applicable in the region where either variable is large, Wadsworth and Tawn, [64] introduced an asymptotic model for the joint survivor function on standard exponential margins. In contrast to the MRV framework, the resulting model provides a useful description of AI distributions, but all AD distributions have the same representation.

More recently, there has been interest in modelling the limiting shapes of scaled sample clouds, or limit sets. The study of limit sets has been around since the 1960s [14, 8], and recent works from Nolde, [42] and Nolde and Wadsworth, [43] have shown that these sets are directly linked to several representations for multivariate extremes. For a given distribution, the limit set is obtained by evaluating the asymptotic behaviour of the joint density function on light tailed margins. Many recent approaches have focused on estimation of the limit set in order to approximate extremal dependence properties; see, for instance, Simpson and Tawn, [55], Wadsworth and Campbell, [63], Majumder et al., [38] and Papastathopoulos et al., [47]. However, the limit set itself does not provide a full description of the asymptotic joint density or distribution, so is less useful from a practical modelling perspective.

To understand the limitations of the methods discussed above, it is instructive to provide an illustration of the joint distribution and density functions on heavy- and light-tailed margins for AI and AD random vectors. All the methods discussed above have equivalent representations in angular-radial coordinates, so without loss of generality, we consider the angular-radial dependence. The first step for most methods for modelling multivariate extremes is to transform variables to common margins. Define

so that $(X_{P},Y_{P})$ and $(X_{E},Y_{E})$ have standard Pareto and exponential margins, respectively. Note that $(X_{E},Y_{E})=(\log(X_{P}),\log(Y_{P}))$, and that the dependence structure or copula of $(X,Y)$ remains unchanged by the marginal transformation [56]. Furthermore, the joint survivor function $\bar{F}_{P}(x,y)=\Pr(X_{P}>x,Y_{P}>y)$ is related to the joint survivor function of $(X_{E},Y_{E})$ by $\bar{F}_{E}(x,y)=\bar{F}_{P}(\exp(x),\exp(y))$. Moreover, if $(X_{E},Y_{E})$ has joint density function $f_{E}(x,y)$, then $(X_{P},Y_{P})$ has joint density $f_{P}(\exp(x),\exp(y))=\exp(-r)f_{E}(x,y)$, where $r=x+y$..

Figure 1.1 shows the joint survivor and density functions for the AD Joe copula, as defined in the Supplementary Material, on standard Pareto and exponential margins. Rays of constant angle on each margin are also shown. On Pareto margins, with the axes shown on a logarithmic scale, lines of constant angle asymptote to lines with unit gradient. As such, the MRV framework provides a description of joint tail behaviour in the region close to the line $\log(y)=\log(x)$, i.e., where $X_{P}$ and $Y_{P}$ are of similar magnitudes. In this region, the contours of the joint density and survivor functions asymptote to a curve of constant shape, which describes the joint extremal behaviour in this region. In contrast, the angular-radial description appears different on exponential margins. For the joint survivor function, the contours of constant probability appear to asymptote towards the line $\max(x,y)=c$ for some constant $c$. Wadsworth and Tawn, [64] showed that is the case for all AD distributions. Informally, this is because for a distribution to be AD, the probability mass must be concentrated close to the line $y=x$, so when the density is integrated to obtain the survivor function, the dominant contribution comes from this region. In contrast to the joint survivor function, the angular-radial description of the joint density is not the same for all AD distributions on exponential margins.

Figure 1.2 shows a similar set of plots for the AI Gaussian copula, also defined in the Supplementary Material. In this case, contours for both the joint density and joint survivor function are curved on both sets of margins. The angular-radial model on Pareto margins describes the section of the curves close to the line $y=x$, which asymptote to straight lines as $x+y\to\infty$. Therefore, the HRV description of the asymptotic behaviour is effectively a straight line approximation to a curve, and is only applicable in the region close to the line $y=x$; see Mackay and Jonathan, [36] for details. In contrast, the angular-radial description of both the density and survivor functions on exponential margins provides a more useful description of asymptotic behaviour. That is, the representation on exponential margins is valid for the full angular range, whereas the representation on Pareto margins is only valid in the joint exceedance region where we are unlikely to observe the largest values of either variable for variables which are AI.

In some applications it is useful to describe the extremal behaviour of a random vector for both large and small values of certain variables; see Section 1.2. In this case, it is more useful to work on symmetric two-sided margins, rather than one-sided margins. Figure 1.3 shows the joint survivor and density functions for a Gaussian copula on standard Laplace margins. The angular-radial variation of the joint survivor function is useful in the first quadrant of the plane, but is less useful in the other quadrants. In the second and fourth quadrants, the contours of the joint survivor function asymptote to the corresponding marginal levels, providing no information about the asymptotic behaviour of the distribution in this region. In contrast, the joint density function provides useful asymptotic information in all regions of the plane.

This motivates an intuitively-appealing angular-radial description of the joint density function, referred to as the semi-parametric angular-radial (SPAR) model [32], which we consider in detail in this article. A similar model was recently proposed by Papastathopoulos et al., [47], although the application was only considered for standard Laplace margins. However, the SPAR framework can be applied on any type of margin. Mackay and Jonathan, [36] showed that on heavy-tailed margins, SPAR is consistent with the MRV/HRV frameworks, and on light-tailed margins, SPAR is consistent with limit set theory. However, the SPAR framework is more general than limit set theory, as it provides an explicit model for the density in extreme regions of the variable space. Moreover, there are distributions which have degenerate limit sets in some regions, for which there is still a useful SPAR representation.

In the SPAR framework, variables are transformed to angular-radial coordinates, and it is assumed that the conditional radial distribution is in the domain of attraction of an extreme value distribution. This implies the radial tail conditioned on angle can be approximated by a non-stationary generalised Pareto (GP) distribution. The SPAR approach generalises the model proposed by Wadsworth et al., [65], in which angular and radial components are assumed to be independent. In the Wadsworth et al., [65] model, the margins and angular-radial coordinate system are selected so that the assumption of independent angular and radial components is satisfied. The SPAR framework removes this requirement, providing a more flexible representation for multivariate extremes.

While a strong theoretical foundation for the SPAR model is provided in Mackay and Jonathan, [36], inference for this model has not yet been demonstrated. Inference via this framework would offer advantages over many existing approaches, and a fitted SPAR model could be used to estimate extreme quantities commonly applied in practice, such as risk measures [41] and joint tail probabilities [27].

The SPAR model reframes multivariate extreme value modelling as non-stationary peaks over threshold (POT) modelling with angular dependence. Many approaches have been proposed for non-stationary POT inference e.g. Randell et al., [51], Youngman, [71], Zanini et al., [73]. In this paper, we introduce an ‘off-the-shelf’ inference framework for the SPAR model. This framework, which utilises generalised additive models [GAMs; 68] for capturing the relationship between radial and angular components, offers a high degree of flexibility and can capture a wide variety of extremal dependence structures, as demonstrated in Sections 6 and 7. Our approach offers utility across a wide range of applications and provides a convenient, practical framework for performing inference on multivariate extremes. Moreover, our inference framework is ready to use by practitioners; open-source software for fitting the SPAR model is available at https://github.com/callumbarltrop/SPAR. For ease of discussion and illustration, we restrict attention to the bivariate setting throughout, noting that the SPAR model is not limited to this setting.

To demonstrate the practical applicability of our proposed inference framework, we consider three bivariate metocean time series made up of zero-up-crossing period, $T_{z}$, and significant wave height, $H_{s}$, observations. We label these data sets as A, B and C, with each data set corresponding to a location off the coast of North America. data sets A and B were previously considered in a benchmarking exercise for environmental contours [19]. Observations were recorded on an hourly basis over 40, 31 and 42 year time periods for data sets A, B and C, resulting in $n_{A}=320740$, $n_{B}=241815$ and $n_{C}=328247$ observations, respectively, once missing observations are taken into account. Exploratory analysis indicates the joint time series are approximately stationary over the observation period. Understanding the joint extremes of metocean variables is important in the field of ocean engineering for assessing the reliability of offshore structures. Wave loading on structures is dependent on both wave height and period, and the largest loads on a structure may not necessarily occur with the largest wave heights. Resonances in a structure may result in the largest responses occurring with either short- or long-period waves, meaning it is necessary to characterise the joint distribution in both of these ranges. These data sets are illustrated in Figure 1.4.

Metocean data sets of this type can often exhibit complex dependence structures, for which many multivariate models fail to account. For example, data set B exhibits clear asymmetry in its dependence structure. Moreover, as demonstrated in Haselsteiner et al., [19], many existing approaches for modelling metocean data sets perform poorly in practice, often misrepresenting the joint tail behaviour or not offering sufficient flexibility to capture the complex data structures. These shortcomings can have drastic consequences if fitted models are used to inform the design bases for offshore structures, as is common in practice.

This paper is structured as follows. In Section 2, we briefly introduce the SPAR model and outline our assumptions. In Section 3, we introduce a technique to estimate the density of the angular component. In Section 4, we introduce a framework for estimating the density of the radial component, conditioned on a fixed angle. In Section 5, we introduce tools for quantifying uncertainty and assessing goodness of fit when applying the SPAR model in practice. In Section 6 and 7, we apply the proposed framework to simulated and real data sets, respectively, illustrating the proposed framework can accurately capture a wide range of extremal dependence structures for both prescribed and unknown marginal distriutions. We conclude in Section 8 with a discussion and outlook on future work.

Let $(X,Y)$ denote a random vector in $\mathbb{R}^{2}$ with continuous joint density function $f_{X,Y}$ and simply connected support containing the point $(0,0)$. The SPAR model for $f_{X,Y}$ requires a transformation from Cartesian to polar coordinates. Polar coordinates can be defined in various ways; see Mackay and Jonathan, [36] for discussion. In this paper, we restrict attention to two particular angular-radial systems corresponding to the $L1$ and $L2$ norms, defined as $\|(x,y)\|_{p}:=(|x|^{p}+|y|^{p})^{1/p},$ $p=1,2$, for $(x,y)\in\mathbb{R}^{2}$. We define $R_{p}:=\|(X,Y)\|_{p}$, $p=1,2$, and consider these variables as radial components of $(X,Y)$. Such definitions of radial variables are common in multivariate extreme value models [e.g., 10, 65]. When using the $L2$ norm to define the radial variable, the corresponding angular variable is usually defined as $\Theta=\text{atan2}(X,Y)$, where atan2 is the four-quadrant inverse tan function. The map between $(X,Y)$ and $(R_{2},\Theta)$ is bijective on $\mathbb{R}^{2}\setminus\{(0,0)\}$. When using the $L1$ norm to define radii, the angular variable is typically defined as $W:=X/\|(X,Y)\|_{1}$ [e.g., 53, chapter 5]. The random vector $(R_{1},W)$ has a one-to-one correspondence with $(X,Y)$ in the upper half of the plane ($Y\geq 0$), but the use of the vector $(R_{1},W)$ becomes ambiguous if we are interested in the full plane, since $W$ contains no information about the sign of $Y$.

With this in mind, we follow Mackay and Jonathan, [36] and define the bijective angular functions $\mathcal{A}_{p}:\mathcal{U}_{p}\to(-2,2]$, where $\mathcal{U}_{p}:=\{(u,v)\in\mathbb{R}^{2}\mid\|(u,v)\|_{p}=1\}$ is the unit circle for the $L_{p}$ norm. For $p=1,2$, these are defined as

where $\varepsilon(v)=1$ for $v\geq 0$ and $-1$ otherwise, is the generalised signum function. The functions $\mathcal{A}_{p}(u,v)$ give a scaled measure of the distance along the unit circle $\mathcal{U}_{p}$ from the point $(1,0)$ to $(u,v)$, measured counter-clockwise.

With angular functions established, we define the angular variables of $(X,Y)$ to be $Q_{p}:=\mathcal{A}_{p}(X/R_{p},Y/R_{p})$, $p=1,2$. The corresponding radial-angular mapping $t:\mathbb{R}^{2}\setminus\{(0,0)\}\to(0,\infty)\times(-2,2]$ given by

is bijective for $p=1,2$. Consequently, we can recover $(X,Y)$ from its radial and angular components, i.e., $(X,Y)=R_{p}\mathcal{A}^{-1}_{p}(Q_{p})$ for $p=1,2$. We note that $Q_{2}=2\Theta/\pi$. However, we use the variable $Q_{2}$ here, in preference to $\Theta$, so that the angular range is the same for both $Q_{1}$ and $Q_{2}$. The joint density of $(R_{p},Q_{p})$ can be written in terms of the joint density of $(X,Y)$,

where the terms $r_{1}$ and $(\pi r_{2})/2$ are the Jacobians of the respective transformations. For ease of notation, we henceforth drop the subscripts on the radial and angular components and simply let $(R,Q)$ denote one of the coordinate systems, with corresponding joint density function $f_{R,Q}$.

Mackay and Jonathan, [36] showed that the choice of coordinate system does not affect whether the SPAR model assumptions (discussed below) are satisfied. However, the coordinates may affect the inference, so in the examples presented in Sections 6 and 7, we consider both $L1$ and $L2$ polar coordinates.

Applying Bayes theorem, the joint density $f_{R,Q}$ can be written in the conditional form $f_{R,Q}(r,q)=f_{Q}(q)f_{R_{q}}(r\mid q),$ where $f_{Q}(q)$ denotes the marginal density of $Q$, $R_{q}:=R\mid~{}(Q~{}=~{}q),$ $q\in(-2,2]$ and $f_{R_{q}}(r\mid q)$ denotes the density of $R_{q}$, with corresponding distribution function $F_{R_{q}}(r\mid q)$. Viewed in this way, the modelling of joint extremes is reduced to the modelling of the angular density, $f_{Q}$, and the tail of the conditional density, $f_{R_{q}}$.

Given any $\gamma\in(0,1)$, define $u_{\gamma}:(-2,2]\to\mathbb{R}_{+}$ as $u_{\gamma}(q)=\inf\{r\mid F_{R_{q}}(r\mid q)\geq\gamma\}$ for all $q\in(-2,2]$, implying $\Pr(R_{q}\leq u_{\gamma}(q))=\gamma$. We refer to $u_{\gamma}(q),q\in(-2,2]$ as the threshold function henceforth. For the SPAR model, we assume that for all $q\in(-2,2]$, there exists a normalising function $c_{q}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ such that

as $\gamma\to 1^{-}$, with $\xi(q)\in\mathbb{R}$. The right hand side of equation (2.2) denotes the cumulative distribution function of a generalised Pareto (GP) distribution, and we term $\xi(q)$ the shape parameter function. The case $\xi(q)=0$ can be interpreted as the limit of equation (2.2) as $\xi(q)\to 0$. Assumption (2.2) is equivalent to the assumption that $R_{q}$ is in the domain of attraction of an extreme value distribution [1]. Given the wide range of univariate distributions satisfying this assumption, it is reasonable to expect the convergence of (2.2) to hold in many cases for $R_{q}$ also. Mackay and Jonathan, [36] showed that this assumption holds for a wide variety of theoretical examples.

This convergence motivates a model for the upper tail of $R_{q}$. Assuming that equation (2.2) approximately holds for some $\gamma<1$ close to $1$, we have

for some $\tau(q)\in\mathbb{R}_{+}$ which we refer to as the scale parameter function. The inclusion of the scale parameter removes the need to estimate the normalising function $c_{q}$, and this is equivalent to the standard peaks over threshold approximation used in univariate extreme value theory [9].

Given $q\in(-2,2]$ and $r\geq u_{\gamma}(q)$, assumption (2.3) implies that

where $\bar{F}_{-}(\cdot):=1-F_{-}(\cdot)$ denotes the survivor function. The joint density of $(R,Q)$ in the region $\mathcal{U}_{\gamma}:=\{(r,q)\in(0,\infty)\times(-2,2]\mid r\geq u_{\gamma}(q)\}$ is then given by

where $f_{GP}$ is the GP density function. Equation (2.4) implies that the SPAR model is defined within the region $\mathcal{U}_{\gamma}$.

To simplify the inference, we also assume that the functions $f_{Q}(q)$, $u_{\gamma}(q)$, $\tau(q)$ and $\xi(q)$ are finite and continuous over $q\in(-2,2]$ and satisfy the periodicity property $\lim_{q\to-2^{+}}f(q)=f(2)$. Such conditions are not guaranteed in general, and whether they are satisfied depends on the choice of margins, alongside the form of the dependence structure. Mackay and Jonathan, [36] showed that the assumptions are valid for a wide range of copulas on Laplace margins, but using one-sided margins (e.g., exponential) or heavy-tailed margins can result in the assumptions not being satisfied for the same copulas.

The SPAR model does not require specific marginal assumptions, and SPAR representations exist for variables with different marginal domains of attraction; however we consider these characteristics unlikely for phenomena in the Earth’s environment. In applications, we typically assume either (i) a practical environmental setting, in which it is reasonable to assume that all variables are bounded (and then apply the model to standardised variables with zero mean and unit variance), or (ii) make marginal transformation to common scale. As discussed in [36], there are theoretical reasons to prefer transformation to Laplace margins.

We compute local stationary estimates at a fixed grid of values $\mathcal{Q}_{grid}:=\{-2+4i/M\mid i=1,2,\dots,M\}\subset(-2,2]$, where $M$ denotes some large positive integer, selected to ensure $\mathcal{Q}_{grid}$ has sufficient coverage on $(-2,2]$. For each $q\in\mathcal{Q}_{grid}$, we assume there exists a local neighbourhood $\mathcal{Q}_{q}=[q-\delta,q+\delta]$, $\delta>0$, such that the distribution of $R_{q^{*}}$ is stationary for $q^{*}\in\mathcal{Q}_{q}$. This is true in the limit as $\delta\to 0$, and a reasonable approximation for small $\delta$.

In practice, rather than fixing the size of the interval, we select the $N$ nearest observations in terms of the angular distance from $q$, defined as $d(q)_{i}:=\min\{\lvert q_{i}-q\rvert,4-\lvert q_{i}-q\rvert\}$, $i=1,...,n$, for some value $N\ll n$. Define $\mathcal{I}_{q}\subset\{1,2,\dots,n\}$ to be the index set of the $N$ smallest order statistics of $d(q)_{i}$. Local estimates of the threshold and parameter functions can be obtained from the corresponding radial set $\mathcal{R}_{q}:=\{r_{i}\mid i\in\mathcal{I}_{q}\}$. Specifically, we define $\hat{u}^{l}_{\gamma}(q)$ to be the $\gamma$ empirical quantile of $\mathcal{R}_{q}$, and $\hat{\sigma}^{l}(q)$ and $\hat{\xi}^{l}(q)$ to be maximum likelihood estimates of the GP distribution parameters obtained from the set $\{r_{i}-\hat{u}^{l}_{\gamma}(q_{i})\mid i\in\mathcal{I}_{q},r_{i}>\hat{u}^{l}_{\gamma}(q_{i})\}$. Choosing an appropriate value for $N$ involves a bias-variance trade-off; selecting too large (small) a value will increase the bias (variability) of the resulting pointwise threshold and parameter estimates. For our modelling procedure, this selection is not crucial, since local estimates are merely used as a means to inform the smooth estimation procedure presented in Section 4.2.

We now consider smooth estimation of the threshold and parameter functions. For this, we employ the approach of Youngman, [71], in which GAMs are used to capture covariate relationships; software for this approach is given in Youngman, [70]. Our procedure is two-fold; we first estimate the threshold function $u_{\gamma}(q)$ for a given $\gamma$, then estimate the parameter functions $\tau(q)$ and $\xi(q)$ using the resulting threshold exceedances.

This section is structured as follows. First, we provide a high-level overview of GAM-based modelling techniques. We then introduce procedures for estimating the threshold and parameter functions via the GAM framework. Finally, we discuss the selection of the basis dimensions required for the GAM formulations.

When applying the SPAR modelling framework in practice, uncertainty will arise for each of the estimated components; namely, the angular density, threshold and parameter functions. In practice, this uncertainty is a result of sampling variability combined with model misspecification. The former arises due to finite sample sizes only partially representing the entire population, while the latter arises from modelling frameworks not fully capturing the complex features of the data. Quantifying this uncertainty is crucial for interpreting statistical results and making informed decisions based on the inherent modelling limitations.

To quantify uncertainty in SPAR model fits, we must consider each model component in turn. For this, we take a similar approach to Haselsteiner et al., [19] and Murphy-Barltrop et al., [41], and consider a fixed angle $q\in\mathcal{Q}_{grid}$, with $\mathcal{Q}_{grid}$ defined as in Section 4.1. We then quantify the estimation uncertainty for each model component while keeping the angle fixed. Specifically, we propose the following bootstrap procedure: for $b=1,\dots,B$, where $B\in\mathbb{N}$ denotes some large positive integer, do the following

1. 1.

   Resample the original data set (with replacement) to produce a new sample of size $n$.

2. 2.

   Compute the point estimate of the angular density at $q$, denoted $\hat{f}_{Q,b}(q)$, using the methodology described in Section 3.

3. 3.

   Compute the point estimates of the threshold, scale and shape parameters at $q$, denoted $\hat{u}_{\gamma,b}(q)$, $\hat{\tau}_{b}(q)$ and $\hat{\xi}_{b}(q)$ respectively, using the methodology described in Section 4.2.

We remark that the choice of resampling procedure can be adapted to incorporate data sets exhibiting temporal dependence. In this case, rather than using a standard bootstrap, one can apply a block bootstrap [29]. This sampling scheme retains temporal dependence in the resampled data set, ensuring the additional uncertainty that arises due to lower effective sample sizes is accounted for [49]. See Keef et al., [27] and Murphy-Barltrop et al., [41] for applications of block bootstrapping within the extremes literature.

Given a significance level $\alpha\in(0,1)$, we use the outputs from the bootstrapping procedure to construct estimates of the median and $100(1-\alpha)\%$ confidence interval for each model component at $q$. Considering the angular density, for example, these quantities are computed using the set $\{\hat{f}_{Q,b}(q)\mid b\leq k\leq B\}$. Assuming the estimation procedure is unbiased, one would expect $\Pr(\hat{f}^{\alpha/2}_{Q,b}(q)\leq f_{Q}(q)\leq\hat{f}^{1-\alpha/2}_{Q,b}(q))\approx(1-\alpha)$, where $\hat{f}^{\alpha/2}_{Q,b}(q)$ and $\hat{f}^{1-\alpha/2}_{Q,b}(q)$ denote the empirical $\alpha/2$ and $1-\alpha/2$ quantile estimates from $\{\hat{f}_{Q,b}(q)\mid 1\leq b\leq B\}$, respectively.

Repeating this procedure for all $q\in\mathcal{Q}_{grid}$ allows one to evaluate uncertainty over the angular domain for each model component, thus quantifying the SPAR model uncertainty. This in turn allows us to evaluate uncertainty in quantities computed from the SPAR model, such as isodensity contours or return level sets; see Sections 6 and 7.

We present a localised diagnostic to assess the relative performance of the SPAR model fits in different regions, similar to that used in [37]. Consider a partition of the angular domain around $q\in\{-1.5,-1,-0.5,0,0.5,1,1.5,2\}$, corresponding to a variety of regions in $\mathbb{R}^{2}$. For each $q$, take the local radial window $\mathcal{R}_{q}$ as defined in Section 4.1. Treating $\mathcal{R}_{q}$ as a sample from $R_{q}$, we compare the observed quantiles with the fitted SPAR model quantiles, resulting in a localised QQ plot. Similar to before, uncertainty can be quantified via non-parametric bootstrapping. Comparing the resulting QQ plots over different angles, and different values of $\gamma$ and $k$, provides another means to assess model performance.

Finally, we propose comparing the estimated angular density, obtained using the methodology of Section 3, with the corresponding density function computed from the histogram. This comparison allows one to assess whether the choice of bandwidth parameter, $h$, is appropriate for a given data structure.

In this section, we evaluate the performance of the smooth inference framework introduced in Section 4.2 via simulation. We do not consider the local estimation approach of Section 4.1 here, since our proposed local estimates are only meant as a means of assessing smooth SPAR estimates when the true values are unknown.

We consider four copulas on standard Laplace margins; Gaussian, Frank, t and Joe, as defined in the Supplementary Material. These distributions represent a range of dependence structures. Note that analogous dependence coefficients to $\chi$ can be defined to quantify the strength of extremal dependence in other regions of the plane [see 36]. In the following, we denote the four quadrants of $\mathbb{R}^{2}$ as $\mathbb{Q}_{1},...,\mathbb{Q}_{4}$. For $\rho>0$, the Gaussian copula has intermediate dependence in $\mathbb{Q}_{1}$ and $\mathbb{Q}_{3}$ [22], and negative dependence in $\mathbb{Q}_{2}$ and $\mathbb{Q}_{4}$. The Frank copula is AI in all quadrants, whereas the t copula is AD in all quadrants. Finally, the Joe copula is AD in $\mathbb{Q}_{1}$, negatively dependent in $\mathbb{Q}_{2}$ and $\mathbb{Q}_{4}$, and AI in $\mathbb{Q}_{3}$. Samples from each copula are shown in Figure 6.1, together with the corresponding values of the copula parameters used in the simulation studies. In each case, the sample size is $n=10,000$. One can observe the variety in dependence structures, as evidenced by the shape of data clouds. For the distributions considered here, the asymptotic shape parameter function is $\xi(q)=0$, $q\in(-2,2]$, and the asymptotic scale parameter functions can be derived analytically [see 36]. The true values of the threshold functions $u_{\gamma}(q)$ and angular density functions $f_{Q}(q)$ can be calculated using numerical integration. Hence, in all cases, the target values for the SPAR model parameters are known.

To evaluate performance, we simulate $500$ samples from each distribution and for every sample, apply the methods outlined in Section 3 and 4.2 to estimate all SPAR model components. Using these estimates, we compute isodensity contours, defined as $\{(x,y)\in\mathbb{R}^{2}\mid f_{X,Y}(x,y)=p\}$ for some $p$. In particular, we consider $p\in\{10^{-3},10^{-4},10^{-5},10^{-6}\}$; the corresponding true contours for each distribution are given in Figure 6.1. These density values represent regions of low probability mass, corresponding to joint extremal behaviour. Moreover, estimates of the joint density are appropriate for evaluating the overall performance of the SPAR model, since in practice, capturing the joint density is crucial for ensuring one can accurately extrapolate into the joint tail.

Alongside isodensity contours, we also compare the estimated GP scale parameter functions and angular density functions to their corresponding target values. For the former, we remark that for each copula, the conditional radial distribution $f_{R_{q}}(r|q)$ only converges to a GP distribution in limit as $r\to\infty$, implying we are unlikely to accurately estimate the asymptotic GP parameter functions for finite samples; see Mackay and Jonathan, [36]. We remark that even with this caveat, we still obtain high quality estimates for the isodensity contours at extreme levels.

Uncertainty in the estimation procedure is quantified by adapting the procedure of Section 5.1 across the $500$ simulated samples. This allows us to compute median estimates and confidence intervals for isodensity contours, scale functions and angular density functions.

Although the choice of coordinate system does not affect whether the SPAR model assumptions are satisfied, the asymptotic SPAR parameter functions do depend on the coordinate system. Since smooth, continuous splines are used to represent the GP parameter functions, the choice of coordinate system may affect the quality of model fits. The simulation study is therefore conducted using both $L1$ and $L2$ polar coordinates.

We first consider the tuning parameters required for the smooth inference procedure, as outlined in Section 4.2.3. For each copula, the threshold and scale parameter functions appeared to vary in a similar manner over angle. Furthermore, we fix a constant value of the shape parameter with angle, i.e., $\xi(q)=\xi\in\mathbb{R}$ for all $q\in(-2,2]$, since for each copula, the tail behaviour remains constant over angles [36]. As discussed in Section 7, this is not true in the general case, so fixing $\xi(q)$ to be constant is an additional constraint imposed on the model.

As noted previously, it is better to select a basis dimension $k$ that is larger than one would expect to be necessary. We considered a range of dimensions in the interval $[5,50]$, and compared the resulting model fits across each of the four copulas. From this analysis, we found that setting $k=25$ was sufficiently flexible for capturing the angular dependence for both the threshold and scale functions.

We are also required to select a non-exceedance probability $\gamma\in(0,1)$. As observed in Mackay and Jonathan, [36], each of the four copulas exhibits a different rate of convergence to the asymptotic form. Therefore, different values of $\gamma$ may be appropriate for these different dependence structures. However, we instead opt to keep $\gamma$ fixed across all copulas. This is for consistency in the estimation framework, as well as to show that even in the case when the model is misspecified, the corresponding inference framework is still robust enough to approximate the true model. We considered a range of $\gamma$ values, restricting our attention to the interval $[0.5,0.95]$, and compared the resulting model fits. As expected, the performance for each copula varied non-homogeneously across $\gamma$ values. Ultimately, we found that setting $\gamma=0.8$ appeared sufficient for approximating the conditional radial tails for each dependence structure. In particular, this value appeared high enough to give approximate convergence to a GP distribution model, without giving a large degree of variability in model estimates.

Finally, for estimation of the angular density, we fix the bandwidth parameter at $h=1/50$ for each copula. Our results show that for these copula examples, this bandwidth is sufficiently flexible to approximate the true angular density functions across the majority of angles.

Figure 6.2 compares the median estimates of isodensity contours, obtained using the $L1$ coordinate system, to the true contours at a range of low density levels; the corresponding plots for the $L2$ coordinate system are given in the Supplementary Material. For both coordinate systems, one can observe generally good agreement between the sets of contours, suggesting the modelling framework is, on average, capturing the dependence structure of each copula. Furthermore, plots comparing the median estimates from both coordinate systems can also be found in the Supplementary Material. These plots show a similar overall performance for both systems, with perhaps a slight preference for the $L1$ estimates.

Plots comparing the estimated GP scale parameter functions, and associated confidence intervals, to the known asymptotic functions are given in the Supplementary Material. In some angular regions, the estimated isodensity contours and scale functions from the SPAR model do not agree with the known values; for example, in the region around $q=0.5$ for the Joe copula. In this case, there is a sharp cusp in the asymptotic GP scale parameter function. Similarly, there is a sharp cusp in the true GP scale parameter for the t copula at $q=\pm 0.5,\pm 1.5$. As the inference framework assumes that the scale is a smooth function of angle, this behaviour is not properly captured. Despite the GAM representation not being able to capture these cusps, the overall performance of the estimated SPAR model is still reasonable in these regions. Furthermore, there is poor agreement between the estimated and asymptotic scale functions for the Frank copula. This is likely due to the relatively slow convergence of this distribution to its asymptotic form, and hence the poorer approximation by the GP distribution.

Next, we consider the uncertainty of the isodensity contours for $p\in\{10^{-3},10^{-6}\}$ in the radial-angular space. Figure 6.3 shows the median contour estimates, along with estimated 95% confidence intervals, obtained under the $L1$ coordinate system; the corresponding plots for $L2$ coordinates are given in the Supplementary Material. One can observe that the true contours are generally well captured within the estimated uncertainty regions. The exception in the $10^{-6}$ density contour for the Frank copula in $\mathbb{Q}_{2}$ and $\mathbb{Q}_{4}$, owing to the aforementioned slow rate of convergence to the asymptotic form for this copula.

Finally, we compare median estimates of the angular density functions, alongside estimated confidence regions, to the corresponding true density functions for each copula. These results are given in the Supplementary Material. Overall, we observe close agreement between the estimated and true functions at the majority of angles. However, we note that the KD estimation framework appeared unable to fully capture the modal regions for the Frank, t and Joe copulas; see Section 8 for further discussion.

Overall, when compared to the truth, the SPAR model estimates perform well for each copula. This observation suggests that our proposed inference framework, with appropriate tuning parameters, can capture the extremal dependence structure across a range of copulas with differing dependence classes. This illustrates both the flexibility and robustness of the SPAR approach, and its advantages over many alternative multivariate models.

The simulation study in Section 6 considered data on standard Laplace margins. This is because the SPAR model assumptions are satisfied for a wide range of copulas on Laplace margins, and the resulting asymptotic representations are relatively simple [36]. However, the SPAR framework does not pre-suppose any particular choice of margins, and Mackay and Jonathan, [36] also showed that SPAR representations arise for random vectors with bounded and heavy tailed margins.

For the case of the metocean data sets, the margins are unknown. In practice, estimation of the marginal distributions, as is common in many extreme value analyses, introduces a high degree of additional modelling uncertainty, and poor marginal estimates affect the quality of the resulting multivariate inference [60]. Therefore, we opt not to transform the margins of the metocean time series, and to instead fit the SPAR model on the original scale of the data. With suitable selections of tuning parameters, we demonstrate below that our inference framework is flexible enough to capture the observed extremal dependence structures for the metocean data sets without the need for marginal transformation.

When modelling phenomena (such as ocean waves) in the Earth’s environment, physical constraints (e.g. limited capacity for energy transfer from surface wind caused by atmospheric low pressure systems, wave steepness limits) typically support the assumption that tails of marginal distributions of random variables are bounded. When variables are presented on different physical scales (e.g. mm vs. km), we are careful to standardise them to zero mean and unit standard deviation prior to analysis. For variables which are believed to have unbounded tails, we would transform to common margins prior to analysis; in this case, as outlined in [36], we favour transformation to common standard Laplace margins.

An important consideration for the SPAR model is where to place the origin of the polar coordinate system. When using Laplace margins, a natural choice is to locate the origin at $(x,y)=(0,0)$. When working on the original scale of the data, the choice is less clear. One option would be to place the polar origin at $(x,y)=(0,0)$. However, this would restrict the range of angles for which the SPAR model offers a useful representation, since both variables we are considering here take only positive values. To account for this, we normalise the data to have zero mean and unit variance, and select the polar origin at $(x,y)=(0,0)$ in the normalised variable space. Define normalised variables $(\tilde{T}_{z},\tilde{H}_{s}):=((T_{z}-\mu_{T_{z}})/\sigma_{T_{z}},(H_{s}-\mu_{H_{s}})/\sigma_{H_{s}})$, where $(\mu_{T_{z}},\mu_{H_{s}})$ and $(\sigma_{T_{z}},\sigma_{H_{s}})$ denote the estimated means and standard deviations of $(T_{z},H_{s})$, respectively.

We henceforth assume that for each data set, the normalised joint density function, $f_{\tilde{T}_{z},\tilde{H}_{s}}$, satisfies the assumptions of the SPAR model, namely that the conditional radial variable converges to a generalised Pareto distribution, in the sense of Equation 2.2, and that the functions $f_{Q}(q)$, $u_{\gamma}(q)$, $\tau(q)$ and $\xi(q)$ satisfy the finiteness and continuity assumptions of Section 2.2. We then apply the statistical techniques introduced in Sections 3, 4 and 5; these results are presented in Section 7.3. Throughout this section, we present all results for the $L1$ coordinate system; the corresponding results for $L2$ coordinates are given in the Supplementary Material, with both systems resulting in similar model fits.

We remark that each metocean time series exhibits non-negligible temporal dependence. Following Section 5.1, we apply block bootstrapping throughout this section whenever quantifying uncertainty, with the block size set to $4$ days. This block size appeared appropriate to account for the observed dependence in each time series. Note that temporal dependence could alternatively be accounted for by ‘declustering’ the data and only modelling peak values. However, in multivariate applications, what constitutes a ‘peak value’ is ambiguous since the extremes of each variable do not necessarily occur simultaneously; see Mackay et al., [34] and Mackay and de Hauteclocque, [33] for further discussion.

Prior to inference, we must first select each of the relevant tuning parameters for the methodologies discussed in Sections 3 and 4. Since the true dependence features are unknown, we use local estimates of the SPAR model, obtained using the framework of Section 4.1, to inform the choice of basis dimensions for the smooth inference procedure of Section 4.2.

To obtain local estimates, we are required to specify the number of reference angles $M$, the number of order statistics $N$, and the non-exceedance probability $\gamma$. For the first two values, we set $M=200$ and $N=500$; we found these values to be adequate to ensure a high degree of coverage over the angular interval $(-2,2]$, and to give angular windows that appeared approximately stationary. For selecting $\gamma$, we tested a range of probabilities in the interval $[0.5,0.95]$. Through this testing, we set $\gamma=0.7$, since this value appeared to give approximate convergence to a GP tail across the majority of local windows. The same non-exceedance probability is also used for the smooth SPAR model estimates.

With local estimates obtained, we then consider smooth estimation of the SPAR components. Notably, for each of the time series, we observe clear trends in the locally estimated shape parameter function; it would therefore not be appropriate to specify this parameter as constant. This makes sense when one considers the shapes of the data clouds illustrated in Figure 1.4; the radial behaviour varies significantly over angles. Furthermore, both metocean variables are bounded below by $0$; therefore, we would expect shorter tails in angular directions that intersect the axes. A range of basis dimensions were tested, and from this analysis, we fixed $k=35$ for the threshold and scale functions, and $k=12$ for the shape function. These values appeared to offer adequate flexibility for capturing the trends observed over the angular variable.

Finally, for angular density estimation, we follow Section 6.2 and set the bandwidth $h=1/50$. This value offered sufficient flexibility to capture the observed angular distributions.

Figure 7.1 compares the threshold and parameter function estimates for data set B from the local and smooth inference procedures; the corresponding plots for data sets A and C are given in the Supplementary Material. The shaded regions in this figure denote the 95% bootstrapped confidence intervals for the smooth model fits. One can observe generally good agreement for each component of the SPAR model. We remark that the local estimates appear unstable, and hence unreliable, for certain angles; this is not surprising, given the small sample size of the angular window. However, the general overall agreement suggests the smooth SPAR estimates are accurately capturing the observed dependence features for each data set, providing evidence that the chosen tuning parameters are appropriate.

Following Section 5.2, we compare the median angular density functions from the KD estimation technique with empirical histograms. These comparisons are given in Figure 7.2 for each data set, and one can observe good agreement between the estimated quantities.

Estimates of isodensity contours are shown in Figure 7.3. These joint density contours are given on the original scale of the data, rather than on the normalised scale, and we consider the joint density levels $p\in\{10^{-3},10^{-6}\}$, corresponding to regions of low probability mass. The estimated isodensity contours appear to capture the shape and structure of each data cloud well. Furthermore, we note that the SPAR model appears able to capture the observed asymmetric dependence structures, illustrating the flexibility and robustness of this modelling framework.

To further demonstrate the utility of the SPAR framework, we also use the fitted model to obtain return level sets for each of the data sets. Return level sets are commonly used in ocean engineering for the design of offshore structures. They can be defined in various ways, but are generally defined in terms of marginal probabilities under various rotations of the coordinate axes, or in terms of the probability of an observation falling anywhere outside the set (the so-called ‘total exceedance probability’ of a contour); see Mackay and Haselsteiner, [35]. Papastathopoulos et al., [47] noted that SPAR-type models offer a natural way for total exceedance probability contours to be constructed. For exceedance probability $a\in[0,1]$ with $a<1-\gamma$, the radius of the contour at angle $q$ is the $(1-a-\gamma)/(1-\gamma)$ quantile of the GP distribution with parameter vector $(u_{\gamma}(q),\xi(q),\tau(q))$. For any angle $q\in(-2,2]$, the probability of an observation exceeding this radius is equal to $a$; consequently, the probability of observing data outside of the resulting contour set is equal to $a$. When observations are independent and the distribution is stationary, we can define such sets in terms of return periods; given a number of years $K\in\mathbb{N}$, the $K$-year return level set is the set corresponding to the probability $a:=1/n_{y}K$, where $n_{y}$ denotes the number of observations per year. One would expect to observe data points outside of the return level set once, on average, every $K$ years. Given the temporal dependence observed within the metocean data sets, we note that such an interpretation is not possible due to clustering of extreme events, and as such these estimates are conservative [34]. However, the resulting return level sets can still provide a useful summary of joint extreme behaviour.

Plots of estimated median $10$ year return level sets for each data set are given in Figure 7.4, along with 95% bootstrapped confidence intervals. These sets, obtained by computing GP distribution quantiles from the fitted model, appear sensible in shape and structure when compared to the data cloud. Moreover, given the lengths of observation windows of each data set, we would not expect to observe many datapoints outside of the return level set; this is clearly true in every case. Furthermore, a comparison of return level sets from the two coordinate systems is given in the Supplementary Material, where one can observe generally good agreement between the estimated sets

We note that simulation from the SPAR model is straightforward; a simulation scheme is given in the Supplementary Material, alongside examples for each of the metocean data sets.

Inspection of the local (angular) QQ plots introduced in Section Section 5.2 suggests that the performance of the fitted SPAR model varies across different angular regions. Although we observe generally good agreement between quantiles, there is clearly better agreement for certain angles, suggesting rates of convergence to the GP tail model may vary over angle for these data sets.

Not Applicable

The data sets analysed in Section 7 are freely available online at https://github.com/EC-BENCHMARK-ORGANIZERS/EC-BENCHMARK.

The authors have no relevant financial or non-financial interests to disclose.

This work was supported by EPSRC grant numbers EP/L015692/1 and EP/Y016297/1.

All authors contributed to the study conception and design. Model development was performed by Ed Mackay and Phil Jonathan. Material preparation and initial analysis were performed by Callum Murphy-Barltrop. The first draft of the manuscript was written by Callum Murphy-Barltrop and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

This paper is based on work partly completed while Callum Murphy-Barltrop was part of the EPSRC funded STOR-i centre for doctoral training (EP/L015692/1). Ed Mackay was funded by the EPSRC Supergen Offshore Renewable Energy Hub, United Kingdom (EP/Y016297/1). We would like to thank Ben Youngman for his assistance with the evgam package in the R computing language.

Supplementary Material to ‘Inference for bivariate extremes via a semi-parametric angular-radial model’