Modelling non-stationarity in asymptotically independent extremes

by   C. J. R. Murphy-Barltrop, et al.

In many practical applications, evaluating the joint impact of combinations of environmental variables is important for risk management and structural design analysis. When such variables are considered simultaneously, non-stationarity can exist within both the marginal distributions and dependence structure, resulting in complex data structures. In the context of extremes, few methods have been proposed for modelling trends in extremal dependence, even though capturing this feature is important for quantifying joint impact. Motivated by the increasing dependence of data from the UK Climate Projections, we propose a novel semi-parametric modelling framework for bivariate extremal dependence structures. This framework allows us to capture a wide variety of dependence trends for data exhibiting asymptotic independence. When applied to the climate projection dataset, our model is able to capture observed dependence trends and, in combination with models for marginal non-stationarity, can be used to produce estimates of bivariate risk measures at future time points.



page 1

page 2

page 3

page 4


Spatial deformation for non-stationary extremal dependence

Modelling the extremal dependence structure of spatial data is considera...

Mar-Co: a new dependence structure to model match outcomes in football

The approaches commonly used to model the number of goals in a football ...

Detecting and modeling worst-case dependence structures between random inputs of computational reliability models

Uncertain information on input parameters of reliability models is usual...

Modelling Extremes of Spatial Aggregates of Precipitation using Conditional Methods

Inference on the extremal behaviour of spatial aggregates of precipitati...

Non-parametric multimodel Regional Frequency Analysis applied to climate change detection and attribution

A recurrent question in climate risk analysis is determining how climate...

Generalised Joint Regression for Count Data with a Focus on Modelling Football Matches

We propose a versatile joint regression framework for count responses. T...

Explaining predictive models using Shapley values and non-parametric vine copulas

The original development of Shapley values for prediction explanation re...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Modelling joint tail behaviour of multivariate datasets is important in a wide variety of applications, including nuclear regulation (OfficeforNuclearRegulation2018), neuroscience (Guerrero2021) and flood risk analysis (Gouldby2017). When analysing multivariate extremes, it is important to capture the dependence structure at extreme levels appropriately. In certain applications, one would expect the extremes to occur simultaneously – a situation termed asymptotic dependence – whilst in others, joint occurrence of the very largest events cannot happen – a situation termed asymptotic independence. Section 2 explains these concepts in detail. The study of extremal dependence structures is well established, and a wide range of statistical modelling techniques have been proposed (Coles1991; Ledford1997; Heffernan2004).

Extremal dependence between two variables may be summarised by bivariate risk measures. A variety of risk measures have been proposed in the literature (Serinaldi2015), and are selected according to the needs of an analysis. For this paper, we restrict attention to one particular measure known as the return curve due to its practical importance (Murphy-Barltrop2021)

. Given a small probability

, the -probability return curve is given by , with corresponding return period . This curve directly extends the concept of a return level from the univariate framework (Coles2001) to the bivariate setting. These curves are widely used in practice to derive extremal conditions during the design analysis of many ocean and coastal structures, including oil rigs (Jonathan2014), railway lines (Gouldby2017) and wind turbines (Manuel2018).

However, in many real world scenarios, datasets exhibit non-stationarity; this feature can result in extremal dependence structures that are not fixed due to covariate influences on the underlying processes. In this setting, there is no longer a meaningful or fixed definition of a return curve. We therefore expand the definition of this risk measure to be covariate-dependent, resulting in a non-stationary counterpart. Given a process with covariates , we define the -probability return curve at a covariate realisation to be .

In a practical setting, we wish to derive estimates of non-stationary return curves for environmental datasets. As our motivating example, we consider relative humidity and temperature data for summer months (June, July and August) obtained from the UK Climate Projections (UKCP18) under emissions scenario RCP 8.5. This corresponds to the ‘worst-case’ scenario, whereby greenhouse gas emissions continue to rise throughout the 21st century (UKMetOffice2018). As such, data from these projections can be used as a risk management tool to help mitigate against the impacts of climate change in a conservative manner. Data is only considered for summer months since extremal dependence structures vary significantly across meteorological seasons; see the Supplementary Material for further details.

Our particular dataset corresponds to temperature and relative humidity projections at a grid cell containing the UK’s Heysham nuclear power station. Denoting relative humidity by for , we define a ‘dryness’ variable as . The dataset of temperature and dryness at the start and end of the time period, along with the temperature time series, are plotted in the left and centre panels of Figure 1, respectively. Extremely high temperatures combined with extremely low humidity (i.e., high dryness) is relevant in the context of nuclear safety (Knochenhauer2004), and therefore one may wish to use these projections to quantify the joint extremal behaviour over the observation period. This dataset provides the subject of our detailed case study in Section 5. Clear non-stationary trends can be observed within both marginal data sets; these trends are likely a result of seasonal behaviour combined with long term trends due to climate change.

Figure 1: Left: Plot of first and last 10 years of combined temperature projections, given in red and light blue, respectively. Centre: Plot of Heysham temperature time series. Right: Plot of estimates over rolling windows (solid black lines), alongside

pointwise confidence intervals (dotted blue lines).

To assess trends in dependence, we can calculate suitable coefficients using rolling windows of data. The right panel of Figure 1 demonstrates a clear trend in an extremal dependence coefficient labelled (Ledford1996); this measure summarises the dependence between the most extreme observations, with larger values corresponding to a higher degree of positive dependence. Further discussion can be found in Section 2.2. The illustrated trend suggests the probability of extreme observations occurring simultaneously is increasing over time, motivating the need for modelling techniques that can capture trends of this nature.

The majority of existing techniques for modelling multivariate extremes assume stationarity in the joint tail structure. Furthermore, of the approaches that can accommodate non-stationarity, most are suitable only for datasets exhibiting asymptotic dependence, as we discuss in Section 2.3. This is restrictive since in practice, asymptotic independence is often observed (Ledford1996; Ledford1997); this is further evidenced by estimated values for the UKCP18 dataset, which indicate the presence of asymptotic independence, at least throughout most of the observation period.

We propose a new method for capturing non-stationary extremal dependence structures when asymptotic independence is present, based on a non-stationary extension to the Wadsworth2013 modelling framework. In doing so, we are able to evaluate and visualise trends across the entire extremal dependence structure. This is in contrast to other approaches, where implementation may be limited to trends in one-dimensional summary measures, such as the coefficient of tail dependence (Ledford1996) or the extremal coefficient (Frahm2006).

This paper is structured as follows: Section 2

recalls existing methodology for capturing tail behaviour in the stationary and non-stationary settings for both univariate and multivariate random vectors. Section


introduces our modelling approach, which relies on quantile regression techniques to derive two estimators of a quantity describing extremal dependence in a non-stationary setting. We also propose methodology for non-stationary return curve estimation using these estimators. Section

4 details a simulation study, showing our estimators to be approximately unbiased over a range of simulated examples. In Section 5, we apply our model to the UKCP18 dataset. Our approach is able to reveal clear trends in the extremal dependence of this process, and estimates of return curves are obtained. We conclude in Section 6 with a discussion and outlook on future work.

2 Background

2.1 Univariate extreme value theory

In the univariate setting, one of the most popular techniques for capturing tail behaviour is known as the peaks-over-threshold approach, whereby a generalised Pareto distribution (GPD) is fitted to all exceedances of some high threshold. This is justified by the Pickands-Balkema-de Haan theorem (Balkema1974; Pickands1975)

, which states that for a random variable

satisfying certain regularity conditions, there exists a normalising function such that


as ; see also Coles2001. Here,

is the cumulative distribution function of a GPD, with scale and shape parameters,

and , respectively, and . The shape parameter dictates the behaviour of the tail, with , and corresponding to bounded, exponential and heavy tails, respectively. In practice, for an observed random variable with a finite sample size, a high threshold is selected and a GPD is fitted to the positive exceedances: we write .

In many contexts, such as financial and environmental modelling, datasets exhibit non-stationarity, whereby the underlying distribution changes with time or other covariates. In most such cases, we can no longer expect a stationary GPD model to capture the tail adequately. This feature can be present in a range of different forms, as exhibited by the seasonal and long term trends present in the UKCP18 dataset introduced in Section 1. Davison1990 addressed this issue by using covariates to capture trends in the GPD parameters. Given a non-stationary process , with covariates , a non-stationary GPD model is given by


for a sufficiently large threshold . More recent extensions to this model also allow the threshold to be covariate dependent. For example, Kysely2010 and Northrop2011 use quantile regression to estimate a threshold with a constant exceedance probability, whereas Sigauke2017 use a cubic smoothing spline. More flexible approaches have been proposed using generalised additive models (GAMs) to capture non-stationary behaviour in univariate extremes (Chavez-Demoulin2005; Youngman2019). GAMs use smooth functions to capture trends due to covariates, and are less rigid than standard regression models. A wide range of statistical software is available for fitting these models (Youngman2020; Wood2021).

All the approaches discussed thus far can only be used to model non-stationarity in the extremes. For many statistics corresponding to joint tail behaviour, such as return curves, one must also be able to capture non-stationary within the body of the data simultaneously. This is because extremes of one variable may occur with average values of another variable. To address such challenges, a range of pre-processing techniques have been proposed that allow marginal non-stationarity to be captured in the body and tail of a dataset simultaneously (Nogaj2007; Eastoe2009; Mentaschi2016). For these approaches, covariate functions are used to capture and effectively ‘remove’ non-stationarity from the body of the data. Once removed, any remaining trends in the tail can be captured using any of the methods introduced above. For a non-stationary process , with covariates , the general set-up of these models is to assume


with and as linear functions of covariates. Here, the residual process is assumed to be approximately stationary, and assigning a distribution to this yields a likelihood for all parameters; Eastoe2009

, for example, adopt a standard normal distribution, with the option to also include a shape transformation. Covariate functions are selected through an analysis of non-stationary trends within the body.

An alternative is given by Krock2021, who propose a single distribution for capturing non-stationary behaviour in the body and tail simultaneously. This distribution, which extends the model for stationary data proposed in Stein2021a, accounts for both seasonal and long term trends. The idea behind this approach is to provide a ‘one size fits all’ model, hence the same trends are assumed to be present within the extremes as in the rest of the data. This is unlikely to be the case in all practical scenarios; as a result, we prefer to adopt pre-processing techniques.

2.2 Bivariate extreme value theory

We briefly recall approaches to modelling extremes in the bivariate setting. To begin, consider a random vector with respective marginal distribution functions . Consider the conditional probability and define the coefficient . The cases and correspond to the aforementioned asymptotic independence and asymptotic dependence schemes, respectively. This distinction is important since many models are suitable for data exhibiting one scheme only.

For mathematical simplicity in the description of extremal dependence, it is common to consider random vectors with standardised marginal distributions. This is achieved in practice through marginal estimation and application of the probability integral transform.

Classical modelling approaches are based on the framework of multivariate regular variation, and are applicable only to asymptotically dependent data. Given a random vector with standard Fréchet margins, we define the radial and angular components to be and , respectively. We say that is multivariate regularly varying if, for all Borel subsets , we have

for any , where is termed the spectral measure (Resnick1987). This assumption implies that, for large radial values, and are independent. The spectral measure captures the extremal dependence structure of

. It must satisfy the moment constraint

, but has no closed parametric form. All asymptotically independent distributions have a spectral measure placing mass at the endpoints and of the unit interval, which is why this modelling framework is unable to capture tail properties under this scheme (Coles1999). Moreover, it has been shown that assuming the incorrect form of extremal dependence will lead to unsatisfactory extrapolation in the joint tail (Ledford1997; Heffernan2004). This has consequently led to the development of flexible modelling approaches that are able to theoretically capture both extremal dependence regimes.

The first such idea was proposed in Ledford1996; Ledford1997. It is assumed that the joint tail of a random vector with standard exponential margins is given by


where is a slowly varying function at infinity, i.e., for , and . The parameter is termed the coefficient of tail dependence, with and corresponding to asymptotic dependence and , or and , corresponding to asymptotic independence. In Figure 1, our estimates of suggest asymptotic independence is exhibited by the UKCP18 data throughout most of the observation period. In practice, this framework is limited by the fact it only characterises the joint tail where both variables are large, and hence is not applicable in regions where only one variable is extreme.

Alternative characterisations of the joint tail have been proposed to circumvent this issue. Heffernan2004 introduce a general, regression-based modelling tool for conditional probabilities. Given a random vector with standard Laplace margins (Keef2013), it is assumed that normalising functions and exist such that the following convergence holds:

for a non-degenerate distribution function . Both regimes can be captured via the functions and , with asymptotic dependence arising when and . Note that one could instead condition on the event . The functions and are typically estimated parametrically, while the distribution function is estimated non-parametrically. This model has been widely used in practice, with applications ranging from ocean engineering (Ross2020) to coastal flood mitigation (Gouldby2017).

Wadsworth2013 provide an alternative representation for the joint tail using a general extension of the framework described in equation (2.4). Given with standard exponential margins, they assume that for each ,


as , where is slowly varying for each ray and is the termed the angular dependence function (ADF). This function, which describes the dependence structure of the joint tail along the ray , generalises the coefficient , with . Both extremal dependence regimes can be captured by , with asymptotic dependence implying the lower bound for all . Pointwise estimates of the ADF can be obtained in practice via the Hill estimator (Hill1975).

Alongside these approaches, we note that there exist several copula-based models that can theoretically capture both extremal dependence regimes, such as those given in Pauli2002, Wadsworth2017 and Huser2019. However, due to the stronger assumptions about the form of parametric family for the bivariate distribution, we prefer instead to use more flexible modelling techniques.

2.3 Non-stationary extremal dependence

Although many extreme value analyses seek to capture marginal non-stationarity, common practice is to assume stationarity in dependence, often without even assessing this feature. Relatively little consideration has been given to this problem in the literature, and most of the approaches that do exist rely on the multivariate regular variation framework, thereby being restricted to asymptotically dependent data. For example, Mhalla2019a and Mhalla2019 propose semi-parametric models to capture trends in parameters of quantities related to the spectral measure, while Carvalho2014 and Castro-Camilo2018 propose flexible modelling techniques for capturing non-stationary trends in the spectral measure under covariate influence.

Mhalla2019 also propose a technique for data exhibiting asymptotic independence, using GAMs to capture trends in the non-stationary extension to the ADF defined in equation (3.1). Given a non-stationary process with standard exponential margins and an external -dimensional covariate , the extended ADF is assumed to take the semi-parametric form


where is a link function, and are subvectors of , or products of covariates if interactions are considered. The vector gathers linear coefficients whereas denote smooth functions for each ; the parameter vector represents all parameters to be estimated, i.e., and the coefficients from each basis function . Estimation of is via a penalised log-likelihood approach, where roughness penalties are imposed to ensure a smooth model fit. The link function is used, resulting in fitted values contained in the interval . However, this range is restrictive since implies positive extremal association. In practice, Mhalla2019 only applied model (2.6) along the ray , corresponding to modelling non-stationarity in only. In Section 3, we propose a novel modelling framework for capturing extremal dependence trends under asymptotic independence via a non-stationary extension to the ADF. A simulation study in Section 4 demonstrates that our approach often results in less bias and variability compared to the technique given in Mhalla2019.

Non-stationary extensions to the Heffernan2004 model also exist: Jonathan2014b propose smooth covariate functions for and , while Guerrero2021 allow theses parameters to vary smoothly over time for blocks of observations via a penalised log-likelihood. However, we note that conditional extremes techniques have been shown to create additional complexities during implementation, requiring more steps compared to alternative approaches because of the need to condition on each variable being extreme separately; see Murphy-Barltrop2021. Our proposed method is simpler to implement in practice compared to the those derived under this framework.

3 Non-stationary angular dependence function

3.1 Introduction

We describe a non-stationary extension to the ADF of Wadsworth2013, which is the key building block for estimating non-stationary return curves. We assume stationary marginal distributions throughout this section, allowing us to separate out the two forms of trends; further discussion on the separate treatment of these trends can be found in Section 5.

Given a non-stationary process with standard exponential margins and an external -dimensional covariate , we assume that for all and ,


where denotes a slowly varying function and denotes the non-stationary counterpart of the ADF at time . Define : we refer to this variable as the min-projection. Equation (3.1) implies that, for each , and ,

However, unlike its stationary counterpart, the non-stationary ADF cannot be estimated via the Hill estimator; this is because we typically do not have repeated observations for a covariate realisation. Moreover, even with repeated observations, the resulting sample sizes would typically be too small for reliable estimation. As such, a new estimation procedure is required for this function. Given and two quantiles close to one with , we assume the existence of positive sequences and such that


for all . Assuming strict monotonicity of the cumulative distribution function for , we deduce that for all . Furthermore, the quantile being close to one implies values of the sequence are large in magnitude. Under the model assumptions, we can therefore deduce that

which is rearranged to give


for all . Hence, estimates of the sequence lead to a point-wise estimator for the non-stationary ADF at a given angle . We denote this estimator by , and describe improvements to its stability in Section 3.3.

3.2 Estimating quantiles of the min-projection

The sequence corresponds to differences in covariate-varying quantiles of the univariate min-projection for each . Quantile regression methods therefore provide a natural solution to the problem of its estimation. Such techniques have successfully been applied in a variety of contexts, ranging from ecology (Cade2003) to growth charts (Wei2006). Here, we describe the most commonly-used approach in terms of the min-projection variable. Given a value , the -th quantile of is

We assume that the conditional quantile function is linear in , implying , where denotes a vector of coefficients. The vector

is estimated through a minimisation of a suitable loss function; see

Koenker2017 for further details.

There is also a range of literature available on the topic of extremal quantile regression (Chernozhukov2004; Youngman2019; Velthoen2021), applicable when is very close to one. However, many of the derived extreme value laws involve unknown constants which cannot be estimated without strong modelling assumptions (Chernozhukov2017). Furthermore, applying the simple ‘rule-of-thumb’ proposed in Chernozhukov2011 for deciding between extreme value and normal quantile approximations, the quantile levels considered in Sections 4 and 5 would not be deemed extreme enough to warrant the former. Moreover, we have found standard quantile regression techniques are able to produce accurate estimates of non-stationary ADFs across a range of simulated examples.

3.3 Averaging over quantiles

Prior to applying the proposed model, one must first select and for estimating quantiles of

. This selection represents a bias-variance trade off, as is often observed in applications of extreme value theory: quantiles that are not sufficiently extreme (close to one) will induce bias in results, while quantiles that are too large will result in highly variable estimates. Moreover, considering only a single pair of quantiles will lead to higher variability in ADF estimates. To address these issues, we consider a range of quantile pairs simultaneously and compute an average estimator over these values. Specifically, let

be quantiles near one, with for . For each , the pair is used to derive an estimator , as in equation (3.3). Our final estimator is derived to be the average of these:

for all and . In unreported simulations, we have found this aggregated estimator to outperform estimators obtained from any individual pair of quantiles considered. Furthermore, a range of quantile sets were compared for the examples discussed in Section 4, with the resulting ADF estimates showing very little difference in variability or accuracy. Our choices for and are detailed in Section 4.

3.4 Bernstein-Bézier polynomial smooth estimator

One drawback of the average estimator proposed in Section 3.3 is that it is pointwise for each ray . This typically leads to non-smooth estimates of the ADF that one would not expect to observe in practice.

In this section, we extend this estimator to give smooth functional estimates using a parametric family derived from the set of Bernstein-Bézier polynomials. These polynomials have been applied in many approaches to estimate Pickands dependence function (Guillotte2016; Marcon2016; Marcon2017a), a quantity related to the spectral measure which bears many similarities to the ADF (Wadsworth2013). For dimension , the family of Bernstein-Bézier polynomials of degree is defined to be

We use this family of functions to approximate . However, for any , we have for all . As such, this family of polynomials can only approximate ADFs representing non-negative dependence in the extremes. Furthermore, we wish to allow for covariate influence in the dependence structure; this corresponds to covariate influence in the coefficient vector . We therefore propose extending this family of polynomials to the following set:

where denote functions of the covariates. For any , let represent a form of the non-stationary ADF given by this family of functions. Our objective is to find an estimator that minimises the equation


over all rays and for ; this is achieved through estimation of the coefficient functions . Since is unobserved in practice, we consider the objective function


with defining a finite set spanning the interval and denoting the parameter vector corresponding to the coefficient functions . The intuition here is that gives an approximation of the absolute value in (3.4) integrated over and : it is therefore desirable to find a value of which minimises .

To estimate , we must specify the form of coefficient functions. To start, we impose that for all ; any function satisfying these conditions has the property that , corresponding to the theoretical end-points of the ADF: . For , we assume that , where denotes a link-function and denotes a vector of coefficients for each . The entire parameter vector is therefore , with an estimator defined as

Finally, a smooth estimator of is given by .

3.5 Incorporating theoretical properties

Wadsworth2013 show that the ADF is bounded from below, i.e, for all . This bound corresponds to asymptotic dependence, with larger values corresponding to a form of asymptotic independence. Neither of our estimators or are required to satisfy this constraint, and so we impose this bound retrospectively. In each case, we define or as

for all , and similarly for . In unreported results, we found that incorporating this lower bound improves estimation quality within the resulting ADF estimates.

For , we also impose end-point conditions retrospectively: for all . Recall that this requirement is already satisfied by the smooth estimator by setting for all .

3.6 Estimating non-stationary return curves

We now consider the problem of estimating at some fixed using an estimator of the non-stationary ADF. Let denote either the estimator or . Given the set of rays defined in Section 3.4 and a quantile close to one, we let be defined as in equation (3.2). Then, for all , define as

implying . Define . We have

meaning that the set provides an approximation of . Similarly to the estimation of and , as described in Sections 3.3 and 3.4, we denote for each quantile close to one, and take our final estimator of the return curve to be .

4 Simulation study

4.1 Overview

We use simulation to evaluate the properties of the estimators and derived in Section 3. Section 4.2 introduces a range of examples exhibiting non-stationary extremal dependence. The variation in dependence structures allows us to assess the relative strengths and weaknesses of each estimator. In Section 4.3, we consider parametric forms for the covariate and link functions introduced in Section 3. In Section 4.4, we evaluate the bias and variability that arises from and , which we find to outperform the corresponding estimator from Mhalla2019. To the best of our knowledge, these are the only approaches within the literature for capturing trends in extremal dependence under asymptotic independence. Finally, in Section 4.5, we illustrate that our proposed estimators can be used to derive accurate estimates of non-stationary return curves across each of the considered examples.

4.2 Simulated examples of non-stationary dependence structures

We introduce several examples exhibiting non-stationary extremal dependence under asymptotic independence. In each case, the non-stationarity is over the time covariate , with denoting the sample size. The first two examples are obtained using the bivariate normal copula, for which the dependence is characterised by the coefficient . For the first example, we take , so that and , i.e., moving from independence to perfect positive dependence. For the second example, we take , giving and , i.e., moving from a strong negative correlation to independence.

For the third, fourth and fifth examples, we use the inverted extreme value copula (Ledford1997) with logistic, asymmetric logistic and Hüsler-Reiss families, respectively. For the logistic and asymmetric logistic, the dependence is characterised by the parameter , with the degree of positive dependence increasing at approaches . We take , hence moving from strong positive dependence at to close to independence at . The asymmetric logistic distribution also requires two asymmetry parameters (Tawn1988): we fix , noting this does not change the overall trend in dependence. The Hüsler-Reiss family is characterised by the dependence parameter , with independence and complete dependence obtained as approaches and , respectively. We take , resulting in increasing dependence over time.

For the final example, we start with a specified ADF and use a method given in Nolde2021 to construct a copula with this ADF. Given the dependence parameter , we take , which is the ADF of the density proportional to


and simulate from this density using MCMC. Such a distribution does not have exactly exponential margins, however in this case the transformation to exponential margins, via the probability integral transform, yields a density with the same ADF. We refer to this example as the copula of model (4.1) henceforth and set ; this results in a similar dependence trend to the inverted logistic example. Illustrations of the resulting ADFs over time for each example are given in Figure 2.

Figure 2: Illustration of true ADFs over time for each copula example. Colour change is used to illustrate trends in extremal dependence structure over the time frame, with red and blue corresponding to the start and end of time frame, respectively.

4.3 Covariate and link functions

Both estimators derived in Section 3 require estimates of sequences obtained via a quantile regression procedure. For this, one must first specify the functional form for the relationship between quantiles and covariates . Since the data we have simulated has a dependence structure directly related to the covariate , we propose the covariate set . For any given quantile and ray , it is assumed that the quantile function is given by , where . Additional polynomial terms were considered, but we found that a cubic expression was flexible enough to accurately capture quantiles trends for each of the studied examples. Applying the methodology described in Koenker2017, estimates of the sequences described in equation (3.2) can be obtained for any pair of quantiles .

As noted in Section 3.3, the choice of quantile sets did not appear to significantly alter the resulting ADF estimates, so long as a range of different quantile pairs were considered. Therefore, for , we take equally spaced points in the interval , and set for .

For , estimates of non-stationary ADFs can be derived directly using estimated sequences, while specification of a coefficients functions is also required for . Defining , we set , with for each , thereby ensuring positive coefficient functions. We found this form to be flexible enough to capture the range of dependence trends considered in the study, though additional polynomial terms and/or coefficient functions could be added for more complex data structures.

4.4 Simulation study

We evaluate the properties of estimators proposed in Section 3 and compare these to the estimator derived in Mhalla2019, which we denote . For the latter approach, the thresholds used to define exceedances of the min-projection are set to be constant over the covariate space: this assumption is not required for or .

To obtain ADF estimates using , the same covariate set is used as for the quantile regression procedure, i.e., . This link function proposed in Mhalla2019 imposes the lower bound of , but also imposes an upper bound of , meaning negative dependence cannot be captured. To relax this constraint, we instead consider the link function , , when estimating the ADFs for the negatively correlated Gaussian example. This ensures the lower bound is still satisfied for each example, while allowing all forms of extremal dependence to be captured. We also impose the boundary conditions discussed in Section 3.5 for and . Furthermore, to ensure comparability, the final estimator is obtained as an average over equally spaced quantiles in the interval . Such quantiles are used to define threshold exceedances of the min-projection variable prior to GAM fitting. We observe that imposing the upper bound for in the majority of considered examples results in this estimator having a practical advantage over and , since the latter estimators can exceed this upper bound.

To evaluate the performance of the estimators, estimates of the mean integrated squared error (MISE) were obtained using samples from each copula. Given an estimator , the MISE at time is given by

with smaller MISE values corresponding to estimators with lower bias and variance. Three different time points, and , were considered, corresponding to the start, middle and end of the simulated time frame, respectively.

Table 1 gives MISE values for each estimator and copula example at each time point. One can observe the lowest MISE values are always for or ; this is likely due to the reduced variance of these estimators compared to , owing to their semi-parametric forms. We also compared the integrated squared error (ISE) of median estimators, obtained pointwise over for each method; in this case, the median estimator for often gives the lowest ISE values. These additional results can be found in the Supplementary Material.

Copula Times MISE - MISE - MISE -
Gaussian (Positive Correlation) Start 0.03375 0.01012 0.00497
Gaussian (Positive Correlation) Middle 0.00356 0.0021 0.00229
Gaussian (Positive Correlation) End 0.01978 0.00407 0.00143
Gaussian (Negative Correlation) Start 32.24126 35.76231 30.10062
Gaussian (Negative Correlation) Middle 0.08607 0.08646 0.06309
Gaussian (Negative Correlation) End 0.02546 0.00849 0.00611
Inverted Logistic Start 0.02075 0.00541 0.00315
Inverted Logistic Middle 0.00209 0.00107 0.00627
Inverted Logistic End 0.03539 0.00893 0.02744
Inverted Hüsler-Reiss Start 0.05748 0.02017 0.03849
Inverted Hüsler-Reiss Middle 0.00177 0.00065 0.00089
Inverted Hüsler-Reiss End 0.01741 0.0057 0.00367
Inverted Asymmetric Logistic Start 0.02771 0.00602 0.01863
Inverted Asymmetric Logistic Middle 0.00279 0.00142 0.01995
Inverted Asymmetric Logistic End 0.02956 0.00641 0.01406
Copula of model (4.1) Start 0.0156 0.00621 0.00683
Copula of model (4.1) Middle 0.00208 0.00088 0.01775
Copula of model (4.1) End 0.04873 0.01486 0.00392
Table 1: MISE values for each estimator at start, middle and end of simulated time frame. Smallest MISE values in each row are highlighted in bold.

We observe that the MISE values for the negatively correlated Gaussian copula at the start of the interval () are significantly larger than other MISE values. This is due to the fact that for strongly negatively dependent data structures, the value of the true ADF tends towards infinity as approaches for any ray , which is difficult to capture in practice. However, we note that while significant bias appears to exist in ADF estimates, we are still able to obtain accurate return curve estimates, as discussed in Section 4.5.

Focusing on and , the MISE results for and appear comparable; however, this is not the case if we consider the time frame as a whole. For this, we fix rays and compute estimates of the ADF across 250 simulated examples. The median of the ADF estimates, alongside and quantiles, are calculated for each ray and time point and these estimates are then plotted over time for fixed rays. The resulting plots for the inverted logistic copula are illustrated in Figure 3. As can be observed, the median estimates appear very close to the true ADF values for both our estimators proposed in Section 3. The uncertainty is lower for compared to , owing to the former’s semi-parametric form. Furthermore, both of our estimators appear to show lower bias on average compared to , particularly for and . Similar results are obtained for each of the copula examples, with the exception of the negatively correlated Gaussian copula: in this case, significant bias arises nearer the start of the time interval for each estimator, owing to the infinite upper bound for as . The resulting plots can be found in the Supplementary Material. These results suggest higher bias for compared to the estimators proposed in Section 3.

Figure 3: Plots of median and 95% confidence interval estimates over time at rays for the inverted logistic copula. Red lines correspond to true ADF values, while the black, green and blue lines correspond to the median estimates for , and , respectively, with the coloured regions representing the estimated confidence intervals.

Finally, we fix the time points and and evaluate the variability in ADF estimates. Median curve estimates, alongside pointwise confidence intervals, are obtained for each of the copula examples. Similar findings are reported for each of the estimators as have been discussed above. Plots of the estimated median curves and confidence regions for each estimator can be found in the Supplementary Material.

From these results, we conclude that the proposed estimators from Section 3 outperform the GAM-based estimator in a wide range of scenarios. This can be partly explained by the fact these estimators do not make the assumption of a constant threshold over time for the min-projection variable, leading to more realistic estimates. Moreover, we note that these results hold even though has the additional advantage of a link function which bounds the resulting ADF estimates for examples exhibiting positive extremal dependence. We therefore choose not to consider this GAM-based estimator further.

4.5 Return curve estimates

We now consider our ultimate goal of estimating return curves, , for extreme survival probabilities. Curve estimates at under and were obtained for 250 simulated examples from each copula. To give an overall impression of the bias from each estimator, we fix a time point and plot the median of the 250 estimates for . Specifically, since each coordinate of is associated with a ray , we take the median along each ray. Median curve estimates for the negatively correlated Gaussian example are given alongside the true return curves in Figure 4. As can be observed, the estimated curves in each case closely resemble the true return curves. This is even true at the start of the observation period, for which significant bias in ADF estimators was observed. Return curve estimates for each of the other simulated examples were similar in accuracy; these plots can be found in the Supplementary Material. Note than in each case, none of the theoretical properties for return curves derived in Murphy-Barltrop2021 have been applied.

Figure 4: Plots of median return curve estimates over time with for the Gaussian copula with negative correlation. Colour change illustrates extremal dependence trends over time, with red and blue correspond to start and end of time frame, respectively. True curves given in left panel, with estimated curves from the median estimators of and given in centre and right panels, respectively.

5 UKCP18 case study

5.1 Properties of data

We denote the dataset introduced in Section 1 as for , with and corresponding to the temperature and dryness variables, respectively. In this case, we have , corresponding to 100 years of summer projections from June 1981 to August 2080. We treat the time index as a covariate for this data; while this does not correspond to any physical process, it can be used to capture the non-stationarity present in the data. For the marginal time series, empirical evidence indicates the presence of seasonal and long term trends within the main bodies of both variables. Further exploratory analysis suggests the presence of non-stationary behaviour within dryness extremes and that non-stationarity is present within the extremal dependence structure, as evidenced by the trend in in Figure 1. In this section, we attempt to account for all three forms of non-stationary trends and produce return curve estimates up to the end of the observation period.

5.2 Capturing marginal non-stationarity

To capture marginal non-stationarity, we extend the pre-processing technique described in equation (2.3), with the goal of removing any marginal trends from the data. Rather than specifying linear parametric forms for the covariate functions functions, as in Eastoe2009, we instead assume the residual process is a sequence of standard normal variables and allow and to be general smooth functions of covariate vectors. These functions can be estimated using a GAM framework, allowing for flexible modelling of the marginal trends.

The time series in both cases appear to exhibit long term trends in location and scale, along with seasonal behaviour within the former. Therefore, for the location function , and scale function , we take and , respectively, where denotes the day index of the process at time . For example, equals 1, 45, and 90 for June 1st, July 15th and August 30th, respectively. A thin plate regression spline is used for the covariate , while for , a cubic regression spline of dimension (corresponding to the number of data points in each year) is used.

The fitting of the location and scale covariate functions is carried out using the R package mgcv (Wood2021), with a gaulss family. Model optimisation is achieved via restricted maximum likelihood estimation, with smoothness penalties and knot quantities selected automatically using generalised cross validation. Further details for these modelling procedures can be found in Wood2017. The resulting fitted trends are in good agreement with empirical estimates from the marginal time series: see the Supplementary Material for the corresponding plots. Moreover, the transformed series appear stationary, with no obvious long term or seasonal trends in either the location or scale.

With non-stationary trends in both marginal bodies accounted for, residual processes can be obtained through the transformations and , where and denote the estimated covariate functions for and , respectively. Assuming an accurate model fit, these processes should be approximately stationary within the body of the data. However, non-stationary trends are likely to remain in the tails since GAM fitting is driven by the body. Following Eastoe2009, we fit the non-stationary GPD described in equation (2.2) to capture any remaining trends. Significant linear and harmonic trends are shown to exist within the scale parameter of the residual process for dryness, while the shape parameter is assumed to be fixed over time. This assumption is common within the analysis of non-stationary univariate extremes (e.g. Eastoe2009; Chavez-Demoulin2012), since the shape parameter is often seen seen as an intrinsic property of a physical process. No significant trends were found for the scale parameter corresponding to the temperature variable. Let denote a data sample corresponding to the residuals vector . An estimate of the marginal distribution function is given by


where denotes an indicator function, is the empirical quantile of , and are the MLEs of the GPD scale and shape parameters. The stationary GPD, denoted in equation (2.1), is used to estimate the upper tail of .

Finally, the data is transformed to standard exponential margins via the probability integral transform. To assess the outcome of the pre-processing procedure, exponential rate parameters were estimated over time for both marginal processes. The resulting estimates remain approximately constant at one (the target value) throughout the observation period; an illustrative plot can be found in the Supplementary Material.

5.3 Model fitting

With the data transformed to standard exponential margins, we apply the methodology proposed in Section 3 for return curve estimation. The extremal dependence trend observed in Figure 1 appears to suggest the extremes of the process are becoming more positively dependent over the time frame; therefore, one may expect the ADF estimates to tend towards the lower bound as .

We set for both the quantile regression procedure end estimation of Bernstein polynomial coefficients. This reduced covariate space was flexible enough to capture the observed extremal dependence trend within the data. Moreover, the same set of quantile pairs was considered as defined in Section 4.3. Two sets of ADF estimates are obtained over the observation period corresponding to and ; these are plotted in the left and centre panels of Figure 5, respectively. The selected values of for the plotted curves correspond to July 15th for an increasing subset of years between 1981-2080.

Figure 5: Left, centre and right: ADF estimates over time for , and constrained , respectively. Curves change from red to blue over the observation period.

However, we note that for both estimators, we obtain ADFs which exceed one at several rays; this is unexpected, since exploratory analysis of the data appears to suggest the presence of positive dependence over the entire time frame. While or for all indicate positive and negative extremal dependence, respectively, intermediate scenarios where for some rays and for others are theoretically possible. Having investigated this aspect, the evidence that for small is not strong, and so we seek to implement an estimator which is bounded above by one. We therefore re-apply the Bernstein polynomial estimator, this time constraining the coefficient functions to be in the interval

via a logit link function, i.e.,

, , . The resulting ADF estimates are illustrated in the right panel of Figure 5. The estimated trends are again in good agreement with what is observed in the data.

To evaluate uncertainty in estimates, we propose a block bootstrapping procedure. First, the data on exponential margins is split into segments of size , corresponding to five years of observations. The extremal dependence structure within such segments is then assumed to be approximately stationary. Within each segment, data is then resampled in blocks of size with replacement to account for temporal dependence. These blocks are combined to form a resampled segment, then each of the segments are merged in order to obtain a new dataset. This process is repeated times to generate sets of ADF estimates while accounting for complex structures in the data. The median of the estimated ADFs for , alongside pointwise confidence intervals, are illustrated in the left panel of Figure 6 for the Bernstein polynomial estimator with constrained coefficients.

Finally, to assess the quality of the ADF estimates, we have compared the model estimates of over time, obtained via the constrained version of , to the empirical estimates introduced in Section 1. For each rolling window, we have taken the average estimate from the fitted model. As illustrated in the centre panel of Figure 6, the model estimates appear similar over time, suggesting we have accurately captured the extremal dependence trend at this ray.

5.4 Return curve estimates

We use our estimated ADFs to estimate return curves up to the year 2080. This is done in two steps: first, we apply the method introduced in Section 3.6 to obtain return curve estimates on standard exponential margins. These curves are then transformed back to the original scale by applying the the inverse of the empirical distribution given in equation (2.2), followed by the inverses of the transformations used to obtain the variables and . The right panel of Figure 6 illustrates the return curve estimates for the vector , obtained using the constrained version of , with ; this corresponds to an extreme probability, relative to the size of dataset. The selected values of for the plotted curves again correspond to July 15th for an increasing subset of years between 1981-2080. All of the theoretical properties for return curves introduced in Murphy-Barltrop2021 have been imposed to ensure the resulting estimates are both theoretically possible and realistic.

From these curve estimates, two conclusions are evident. Firstly, clear marginal trends can be observed for both time series. For example, the temperature values at the point where the curves intersect the -axis, which equate to -th non-stationary quantile estimates, increase by over C through the observation period. Moreover, a dependence trend is evidenced by the changing shape in return curve estimates over time, with the curves becoming increasingly ‘square’ shaped from 1981-2080. This suggests that joint extremes of temperature and dryness are becoming more likely to occur over the observation period.

Figure 6: Left: median ADF estimate for obtained through bootstrapping procedure, with coloured region representing uncertainty bounds. Centre: Comparison of averaged model estimates for rolling windows to empirical estimates, with black, green and dotted blue lines corresponding to empirical, model, and 95% confidence interval estimates respectively. Right: Return curve estimates on original margins with at July 15th over the observation period. Time is illustrated using a colour transition, with red and blue curves correspond to start and end of time frame, respectively.

6 Discussion

We have proposed a novel method for capturing non-stationary extremal dependence structures under asymptotic independence. Our method has been successfully applied to heavily non-stationary data from the UKCP18, allowing us to obtain return curve estimates up to the year 2080 and thereby illustrating how this framework could help improve practical risk management under future climate scenarios.

We have investigated the properties of our estimators via simulation and have observed them generally perform well in terms of bias. While it would be desirable also to have theoretical results on bias of the estimators, this is very challenging in practice and is likely to require assumptions that are too strict to make the results worthwhile. Potentially a more promising line of future work would be the development of diagnostic plots for non-stationary ADF and return curve estimates. Our diagnostic procedures were limited to comparison of rolling window estimates against those derived from the ADF. In general, diagnostic plots for return curves, such as those introduced in Murphy-Barltrop2021, require stationarity assumptions.

For the sake of simplicity, we have restricted attention to the bivariate setting. However, it is worth noting that all methods could in principle be extended to the general multivariate setting. This scenario results in additional modelling complexities, since different extremal dependence regimes can exist within subvectors of a multivariate random vector. Moreover, applications of multivariate extremal risk measures are limited, owing in part to the fact visualisation and interpretation becomes more challenging in this setting.

As noted in Section 3.5, imposing the lower bound on the estimates of the ADF was shown to improve the estimation procedure. This was an ad-hoc post-processing step, with many of the obtained ADF estimates below the bound. Future work could explore how this lower bound could be incorporated into the theoretical quantile regression framework, and whether this would further improve the quality of ADF estimates. We note that when asymptotic dependence is present, will cease to depend on . However, return curve estimates are also affected by the sequences as described in Section