Understanding the relationship between the extremes of multiple variables, termed extremal dependence, is an important area of research, with applications in environmental science (Brunner et al. 2016), engineering (Tendijck et al. 2021), actuarial science (Quinn et al. 2019) and finance (Nolde and Zhou 2021). Typically, modelling procedures aim to estimate fixed coefficients or parameters related to the extremal dependence structure, which are subsequently used to infer extremal statistics. This field of study, known as multivariate extreme value theory, helps practitioners to carry out robust joint risk assessments across multiple variables, or locations, relevant to a given application.

In the bivariate setting, there exist two distinct subclasses of extremal dependence. Given variables \(X_i \sim F_i, i=1,2\), define the conditional probability

and set \(\chi = \lim _{u \rightarrow 1^-} \chi (u) \in [0,1]\), where this limit exists (e.g., Coles et al. 1999). This dependence measure provides a summary of joint tail behaviour, and when \(\chi > 0\) and \(\chi = 0\), we say that the random vector \(\varvec{X}:= (X_1,X_2)\) exhibits asymptotic dependence (AD) and asymptotic independence (AI), respectively. Notably, classical extremal dependence modelling procedures, which are developed under the framework of regular variation, have been shown to perform poorly for data exhibiting AI (Ledford and Tawn 1997; Heffernan and Tawn 2004). Since AI is commonly observed in many relevant data sets, this is significant drawback of classical approaches, and many recent works have argued against their use when AI is present (e.g., Opitz 2016, Murphy-Barltrop et al. 2024b, Huser et al. 2025). In this paper, we pioneer the development of non-stationary extremal dependence modelling via a geometric approach based on limit sets (Nolde 2014; Nolde and Wadsworth 2022). The framework can accommodate both AD and AI structures, but the modelling choices we subsequently make are better suited to AI. For simplicity, we also restrict attention to the bivariate setting throughout, noting that many of the introduced techniques can, in principle at least, be extended to higher dimensions; see Section 6 for further discussion.

Most modelling frameworks for extremal dependence make the simplifying assumption that data are independent and identically distributed (IID). However, many real-world data sets exhibit non-stationarity, whereby the distribution of the data is not fixed in time. Non-stationarity can be present in two distinct forms: within the marginal distributions \(F_1,F_2\), and within the (extremal) dependence structure of \(\varvec{X}\). We restrict attention to the latter form, noting that stationary margins can be obtained by fitting non-stationary marginal models and transforming the data via the probability integral transform to a standard marginal scale. Such transformations do not alter the dependence structure of \(\varvec{X}\), owing to Sklar’s theorem (Sklar 1959). However, as recently demonstrated in Kakampakou et al. (2024), selecting the appropriate marginal model for non-stationary data is not straightforward, and poor estimates of the marginal distribution can strongly affect the representation of extremal dependence; see also Murphy-Barltrop and Wadsworth (2024).

In the context of extremal dependence, non-stationarity can result in dependence coefficients or parameters that are not fixed in time, requiring modelling procedures that can capture trends in these quantities. Let \(\varvec{X}_t:= (X_{1,t},X_{2,t})\) denote a continuous, non-stationary bivariate process with standardised IID marginal distributions indexed by \(t \in \mathbb {Z}\), where t represents discrete-time, and \(\varvec{Z}_t\) is a continuous, p-dimensional covariate process (the predictor variables) that directly influences \(\varvec{X}_t\) (the response vector). Our interest lies in understanding, and modelling, the relationship between the extremal dependence structure of \(\varvec{X}_t\) and the covariate process \(\varvec{Z}_t\). Relatively few approaches have been proposed in this context, even though non-stationary dependence is a common feature across many relevant data sets (e.g., Jonathan et al. 2014b, Castro-Camilo et al. 2018, Murphy-Barltrop and Wadsworth 2024). We note here that typically we only observe any given non-stationary process for a discrete number of time points; we therefore let \(\{1,\hdots ,T \}\) denote a finite set of arbitrary, fixed indices.

The majority of the available models for non-stationary dependence are targeted at data structures exhibiting AD (see, e.g., de Carvalho and Davison 2014, Mhalla et al. 2017, Castro-Camilo et al. 2018, Mhalla et al. 2019), and will therefore perform poorly for data exhibiting AI. In the AI setting, Jonathan et al. (2014b); Guerrero et al. (2023) and Talento et al. (2025) consider non-stationary extensions of the conditional extremes model proposed by Heffernan and Tawn (2004), while Mhalla et al. (2019); Murphy-Barltrop and Wadsworth (2024) and André et al. (2025) propose non-stationary extensions, and inference techniques, for a quantity termed the angular dependence function (Wadsworth and Tawn 2013). Alongside these approaches, Mhalla et al. (2019) and Lee et al. (2024) model coefficients of extremal dependence in the non-stationary setting; in particular, the limit \(\chi \) and the coefficient of tail dependence \(\eta \in (0,1]\), first proposed by Ledford and Tawn (1996), which quantifies the dependence amongst asymptotically independent data structures; see Appendix B for a formal definition.

The available approaches for modelling non-stationary extremal dependence offer limited utility. For example, coefficients such as \(\eta \) and \(\chi \) merely summarise the extremal dependence, and provide no information outside of the region where all variables are jointly large (i.e., the top right corner of the positive quadrant \(\mathbb {R}^2_+\)). Furthermore, the Heffernan and Tawn (2004) and Wadsworth and Tawn (2013) modelling frameworks suffer from various drawbacks; the former requires one to fit several separate models, which can result in contradictory conclusions (Liu and Tawn 2014; Simpson and Tawn 2024a), while the latter can only be used for specific probability calculations (e.g., joint survivor sets). These limitations, combined with the sparsity of existing approaches, motivate novel developments for modelling non-stationary extremal dependence in the AI setting.

An increasingly popular modelling tool for multivariate extremes involves the use of limit sets, which arise from the asymptotic convergence of scaled sample clouds for suitably chosen scaling sequences and marginal distributions, to infer extremal dependence properties and perform inference. We briefly summarise this convergence property below. Let \(\varvec{X}= (X_1,X_2)\) be a continuous random vector with a stationary copula and marginal distribution functions asymptotically equal to a von Mises function, i.e., light-tailed margins (Nolde and Wadsworth 2022). Consider the set \(C_n:= {\{ \varvec{X}_i/r_n\}^n_{i=1}}\) denoting n independent copies of \(\varvec{X}\) scaled by a suitably chosen sequence \(r_n>0\) satisfying \(r_n \rightarrow \infty \) as \(n \rightarrow \infty \). Under mild conditions, \(C_n\) converges in probability, as \(n \rightarrow \infty \), onto a deterministic set \(\mathcal {G}:= \{\varvec{x}: g(\varvec{x}) \le 1 \}\) (Davis et al. 1988; Kinoshita and Resnick 1991; Balkema and Nolde 2010). The set \(\mathcal {G}\) is known as the limit set, and \(g:\mathbb {R}^2 \mapsto \mathbb {R}_+\) is termed the gauge function of \(\mathcal {G}\).

The shape of \(\mathcal {G}\) is affected by the margins and dependence structure of \(\varvec{X}\). To focus on dependence only, a common approach is to standardize margins. We adopt standard Laplace margins in this article. The theoretical limit sets for three copulas on Laplace margins are illustrated in Fig. 1, alongside the scaled sample clouds \(C_n\) for a large, but finite, sample size n. One can observe that even at the finite level, the scaled observations lie approximately within the interior \(\mathcal {G}\). Further background on limit sets and conditions for convergence of \(C_n\) are discussed in Section 2.2.

Many recent works have introduced techniques for estimating \(\mathcal {G}\). These approaches can be divided into two categories: those that only target limit set estimation, and those that use the geometric framework as a vehicle for full joint tail inference. For the former, Simpson and Tawn (2024a) and Majumder et al. (2025) introduce semi-parametric techniques using generalised additive models (GAMs) and Bézier splines, respectively, to approximate \(\mathcal {G}\). Note that we also employ GAMs and splines for our proposed modelling approach, although in an entirely different context (i.e., non-stationary extremal dependence) requiring more complex formulations. For the latter category, Wadsworth and Campbell (2024) propose a parametric approach, Campbell and Wadsworth (2024) present a semi-parametric piecewise linear model, Papastathopoulos et al. (2023) introduce a latent Gaussian model, and Murphy-Barltrop et al. (2024a) put forward a deep learning approach. In each of these works, the authors perform rigorous simulation studies, and demonstrate that the proposed techniques are able to accurately capture a wide variety of extremal dependence structures encompassing both dependence classes.

It is clear that modelling multivariate extremes through limit sets represents an important line of research. Furthermore, existing estimation approaches for \(\mathcal {G}\) have clearly demonstrated many advantages of the geometric framework; namely, robustness across difference dependence structures (Majumder et al. 2025), self-consistent conclusions (Simpson and Tawn 2024a), and greater applicability in higher dimensions (up to \(d = 10\)) than many competing approaches for modelling extremal dependence (De Monte et al. 2025). However, to date, the literature on limit set estimation has exclusively focused on the case when data are IID, and thus cannot be applied for data exhibiting non-stationary dependence.

In this article, we introduce a novel non-stationary framework for modelling extremal dependence via limit sets. Our approach can be applied to data exhibiting both AD and AI, with more of a focus placed on the latter class. Moreover, inference from our approach is not limited to summary dependence measures or specific joint tail probabilities: for any index t, our model can be used for simulation from the joint tail of \(\varvec{X}_t\), for probability estimation in general joint tail regions, and for approximating bivariate risk measures.

In Section 2, we outline a natural extension of the geometric framework to the non-stationary setting, and detail our asymptotic assumptions. In Section 3, we explain our GAM framework for estimating the non-stationary counterpart of the gauge function (equivalently, the limit set). In Section 4 we apply our framework to a range of simulated data sets, demonstrating the proposed methodology can approximate non-stationary boundary sets across many dependence structures. In Section 5, we apply our model to financial data, with diagnostics indicating good performance. We conclude in Section 6 with a discussion and outlook on future work.

The limit set \(\mathcal {G}\) of the scaled sample cloud \(C_n\) was introduced in Section 1.3. When \(\varvec{X}\) has standard Laplace margins, Balkema and Nolde (2010) and Nolde and Wadsworth (2022) provide a sufficient condition for \(C_n\) to converge onto \(\mathcal {G}\).

(Corollary of Proposition 2.2 of Nolde and Wadsworth 2022). Suppose the joint density, \(f_{\varvec{X}}(\cdot )\), of \(\varvec{X}\) exists, and that

for a continuous gauge function \(g(\cdot )\). Then, \(C_n\) converges in probability onto \(\mathcal {G}\) as \(n\rightarrow \infty \).

The limit set \(\mathcal {G}\) is star-shaped at \(\varvec{0}:=(0,0)\) and compact (Kinoshita and Resnick 1991). It is clear from Eq. 1 that the gauge function \(g(\cdot )\) is 1-homogenous, i.e., \(g(c\varvec{x}) = cg(\varvec{x})\) for any \(c > 0\), \(\varvec{x}\in \mathbb {R}^2\), linking to the star-shaped property. In standard Laplace margins, \(\mathcal {G} \subseteq [-1,1]^2\), with the corresponding limit set possessing componentwise maxima and minima \(\varvec{1}\) and \(-\varvec{1}\), respectively. For any random vector exhibiting AD, such as the student-t copula, we have \(\varvec{1} \in \mathcal {G}\) (see Fig. 1). Here and in what follows, all convergence results involving sets are understood in the Hausdorff sense; we refer to Nolde and Wadsworth (2022) for a detailed discussion and theoretical treatment of limit sets.

Following Nolde (2014), we also define the boundary set \(\partial \mathcal {G}:= \{\varvec{x}: g(\varvec{x}) = 1 \} \subset \mathcal {G}\), noting this set possesses the same componentwise maxima and minima as \(\mathcal {G}\). The boundary sets are also illustrated for the three copula examples in Fig. 1. We remark that one only needs to know \(\partial \mathcal {G}\) in order to specify \(\mathcal {G}\); hence no knowledge is lost through just considering the boundary set. We therefore use the term ‘limit set’ to refer to both \(\partial \mathcal {G}\) and \(\mathcal {G}\) throughout this work.

Recent works have demonstrated that the limit set is directly related to extremal dependence structure of \(\varvec{X}\) (Nolde 2014; Nolde and Wadsworth 2022). Specifically, \(\partial \mathcal {G}\) links several representations for multivariate extremes, immediately providing coefficients and parameters associated with the models proposed in Ledford and Tawn (1996); Heffernan and Tawn (2004); Wadsworth and Tawn (2013) and Simpson et al. (2020). For example, taking \(\eta \in (0,1]\), we have that \(\eta := \max \{s \in [0,1]: [s,\infty )^2 \cap \partial \mathcal {G} \ne \emptyset \}\). In this sense, the study of deterministic geometric shapes is closely linked to the study of extremal dependence, and we refer to this topic as geometric extremes herein. Furthermore, in a practical setting, estimates of \(\partial \mathcal {G}\) can be used to obtain estimates of risk measures, including, but not limit to, joint tail probabilities (Wadsworth and Campbell 2024), return curves (Murphy-Barltrop et al. 2024b), and return level sets (Papastathopoulos et al. 2023). Thus, knowledge of \(\partial \mathcal {G}\) provides not only a great deal of information about the extremal dependence structure, but also supplies practitioners with a framework for performing estimation on multivariate extremes.

A brief review of statistical estimation techniques for \(\mathcal {G}\) was given in Section 1.3. When performing statistical estimation in the geometric extremes framework, it is helpful to decompose a random vector into angular and radial components. However, the method by which one performs this decomposition is ambiguous, owing to the theoretically infinite number of angular-radial systems that are available; see Mackay and Jonathan (2023) for further discussion. As such, we define angular radial systems in a somewhat general manner, and leave the details of our chosen system(s) to Section 3.

Let \(\Vert \cdot \Vert :\mathbb {R}^2 \rightarrow \mathbb {R}_+\) denote any valid norm. Given \(\varvec{X} \in \mathbb {R}^2\setminus \varvec{0}\), define radial R and angular \(\varvec{W}\) variables by \(\varvec{X}\mapsto (R,\varvec{W}):=(\Vert \varvec{X}\Vert ,\varvec{X}/\Vert \varvec{X}\Vert )\), such that \(R >0\) and \(\varvec{W}\in \mathcal {S}^{1}\), where \(\mathcal {S}^{1}:= \{ \varvec{x} \in \mathbb {R}^{2}: \Vert \varvec{x}\Vert = 1\}\) denotes the (closed) unit ball of \(\Vert \cdot \Vert \). Clearly \(\varvec{X} = R\varvec{W}\), implying the behaviour of \(\varvec{X}\) is directly related to the joint distribution of \((R,\varvec{W})\). It is trivial to show that the mapping \(t:\mathbb {R}^2\setminus \varvec{0} \mapsto \mathbb {R}_+ \times \mathcal {S}^{1}\), where \(t(\varvec{x}):= (\Vert \varvec{x} \Vert , \varvec{x}/\Vert \varvec{x} \Vert )\), is bijective; thus, no information is lost through considering \((R,\varvec{W})\). Owing to the star-shaped property of \(\mathcal {G}\) and the homogeneity of \(g(\cdot )\), one can then reformulate the boundary set in terms of the angular component; specifically, we have

This implies that \(\partial \mathcal {G}\) can be specified by only evaluating the gauge function on \(\mathcal {S}^{1}\).

All of the existing limit set estimation approaches make model assumptions about the tail of \((R \mid \varvec{W} = \varvec{w})\), \(\varvec{w} \in \mathcal {S}^{1}\), from which the limit set can be approximated. Treating the angle as a covariate, or predictor variable, this modelling strategy is similar to a non-stationary univariate extreme value analysis, which is also the foundation of the recently-developed SPAR model for multivariate extremes (Mackay and Jonathan 2023), and for which a wide range of modelling approaches are available (see, e.g., Murphy-Barltrop et al. 2025, Mackay et al. 2025).

Let \(\varvec{X}_t\) denote a non-stationary, bivariate process with common margins (e.g., Laplace). To extend the geometric extremes framework to the non-stationary setting, we first require additional notation and assumptions. Letting \(\varvec{z}\in \mathbb {R}^p\) denote any fixed covariate realisation, we make the following assumptions.

The process \(\varvec{X}_t\) is conditionally stationary w.r.t \(\varvec{Z}_t\), i.e., \((\varvec{X}_t \mid \varvec{Z}_t = \varvec{z}) \overset{d }{=} (\varvec{X}_{t'} \mid \varvec{Z}_{t'} = \varvec{z})\) for all \(t,t' \in \{1,\hdots ,T \}\).

The process \(\varvec{X}_t\) is conditionally independent, i.e, \((\varvec{X}_t \mid \varvec{Z}_t) \perp \!\!\!\perp (\varvec{X}_{t'} \mid \varvec{Z}_{t'})\) for all \(t \ne t'\).

The joint density of \((\varvec{X}_t \mid \varvec{Z}_t = \varvec{z})\), denoted \(f_{\varvec{X}_t\mid \varvec{Z}_t}(\cdot \mid \varvec{z})\), exists and satisfies

for all \(t \in \{1,\hdots ,T \}\), where \(g_{\varvec{z}}(\cdot )\) is termed the non-stationary gauge function.

Let \(C^{\varvec{z}}_n:= {\{ \varvec{X}'_t/r_n\}^{i_n}_{t=i_1}}\), where \(\varvec{X}'_t:= \varvec{X}'_t(\varvec{z}) = (\varvec{X}_t \mid \varvec{Z}_t = \varvec{z})\) for all \(t \in \{i_1,\hdots ,i_n\} \subseteq \{1,\hdots ,T \}\) such that \(n \le T\) and \(T \rightarrow \infty \) as \(n \rightarrow \infty \). Assumptions 2.1–2.3 yield the following result.

Suppose that \(\varvec{X}_t\mid \varvec{Z}_t\) obeys Assumptions 2.1–2.3. Then, the scaled sample cloud \(C^{\varvec{z}}_n\) converges in probability onto \(\mathcal {G}_{\varvec{z}}:= \{\varvec{x}: g_{\varvec{z}}(\varvec{x}) \le 1 \}\) as \(n \rightarrow \infty \), for all \(\varvec{z}\in \mathbb {R}^p\).

Under Assumptions 2.1 and 2.2, \(\{\varvec{X}'_t\}^{i_n}_{t=i_1}\) is an IID sequence of n random vectors for all \(\varvec{z} \in \mathbb {R}^p\). Taking Assumption 2.3 and applying Proposition 2.1 pointwisely for all \(\varvec{z} \in \mathbb {R}^p\), the result follows. \(\square \)

As in Section 2.1, we define the boundary set \(\partial \mathcal {G}_{\varvec{z}}:= \{\varvec{x}: g_{\varvec{z}}(\varvec{x}) = 1 \} \subset \mathcal {G}_{\varvec{z}}\), which we term the non-stationary limit set henceforth. Two comments on the connection between theory and practice are in order. First, in practical applications with continuous covariates, one will never observe repeated observations at a fixed realisation \(\varvec{z}\); consequently, our estimation strategy, outlined in Section 3, uses an established regression approach to smoothly approximate \(\partial \mathcal {G}_{\varvec{z}}\) over the covariate domain. Second, standard Laplace marginals for \(\varvec{X}_t\) are a natural choice in practice as they offer a more complete picture of extremal dependence compared to single-tailed marginal distributions, especially in the case of negative dependence (Murphy-Barltrop et al. 2025).

We emphasise that the framework outlined in this section is more general than most existing approaches for modelling non-stationary extremal dependence detailed in Section 1.2. In particular, knowledge of \(\partial \mathcal {G}_{\varvec{z}}\) allows us to estimate a wide range of dependence coefficients, parameters and risk measures. For instance, taking the coefficient of tail dependence \(\eta \in (0,1]\), a non-stationary counterpart is given by \(\eta _{\varvec{z}}:= \max \{s \in [0,1]: [s,\infty )^2 \cap \partial \mathcal {G}_{\varvec{z}} \ne \emptyset \}\). Such quantities are useful for summarising trends in extremal dependence and performing probability estimation. As another example use case, consider the following: given any \(p \in (0,1)\) close to 1, \(\partial \mathcal {G}_{\varvec{z}}\) allows us to specify a set \(\mathcal {A}_{\varvec{z}}^p\) satisfying \(\Pr [\varvec{X}_t \in \mathcal {A}_{\varvec{z}}^p \mid \varvec{Z}_t = \varvec{z} ] = p\). Such sets are known as environmental contours (Haver and Kleiven 2004) or return level sets (Papastathopoulos et al. 2023). In the geometric extremes setting, this set is centred at \(\textbf{0}\) and is computed using high quantiles at a fixed probability level from the conditional variables \((R \mid \varvec{W} = \varvec{w}), \varvec{w} \in \mathcal {S}^{1}\), thus ensuring an equal probability of exceedance in any angular direction. We term \(\mathcal {A}_{\varvec{z}}^p\) a non-stationary return level set, and remark that such sets provide interpretable and intuitive summaries of joint extremal risks that are used to aid with design analysis in practice (Mackay and Haselsteiner 2021); see Section 5 for further discussion.

Finally, the framework outlined in this section can be used for full tail inference on \((\varvec{X}_t \mid \varvec{Z}_t = \varvec{z})\), including, but not limited to, simulation. This represents a significant practical advantage since many applications of multivariate extreme value theory require a means of efficiently simulating synthetic data in the joint tail, e.g., flood risk mitigation (Keef et al. 2013) and catastrophe modelling (Quinn et al. 2019). Moreover, simulated data allows one to estimate non-stationary tail probabilities for \((\varvec{X}_t \mid \varvec{Z}_t = \varvec{z})\); such probabilities are often used to aid with structural engineering (Jonathan et al. 2014a), portfolio management (Nolde and Zhou 2021) and environmental planning (Gouldby et al. 2017).

Having outlined a non-stationary framework for limit sets in Section 2, we now turn our attention to estimation of the non-stationary gauge function. As illustrated in Section 2.1, we are only required to estimate \(g_{\varvec{z}}(\cdot )\) on any unit ball to obtain the boundary set, and we therefore restrict our estimation domain to \(\mathcal {S}^{1}\). Following Section 2.1, we define \(R_t:= \Vert \varvec{X}_t\Vert \) and \(\varvec{W}_t:=\varvec{X}_t/\Vert \varvec{X}_t\Vert \), and turn our interest to understanding the tail of \((R_t \mid \varvec{W}_t = \varvec{w}, \varvec{Z}_t = \varvec{z})\) for any given \(\varvec{w} \in \mathcal {S}^{1}\) and \(\varvec{z}\in \mathbb {R}^p\). Naturally we wish to avoid making strong assumptions about the tail of \((R_t \mid \varvec{W}_t = \varvec{w}, \varvec{Z}_t = \varvec{z})\) while simultaneously avoiding over-fitting. Furthermore, we desire an estimator that allows us to simultaneously regress over the angular and covariate domains.

The framework introduced Section 2.2 is very general, and could represent a large variety of covariate interactions (e.g., continuous and discrete effects, abrupt regime change), making it challenging to specify an framework that can capture all such possible interactions. As such, we opt to simplify the representation for our estimation procedure, and we take two additional steps. First, we define the polar angular variable \(\Phi _t:= {{\,\textrm{atan2}\,}}(X_{2,t},X_{1,t}) \in [0,2\pi )\), where \({{\,\textrm{atan2}\,}}\) denotes the 2-argument arctangent function. Given any norm \(\Vert \cdot \Vert \), one can identify any point \(\varvec{w} \in \mathcal {S}^{1}\) with an angle \(\phi \in [0,2\pi )\): in particular, the transformation \(\varvec{v}:[0,2\pi ) \mapsto \mathcal {S}^{1}\) given by \(\varvec{v}(\phi ) = (\cos (\phi ),\sin (\phi ))/\Vert (\cos (\phi ),\sin (\phi ))\Vert \) is one-to-one. Consequently, we can represent the unit ball for any norm via the polar angular interval \([0,2\pi )\), and no information is lost by considering \(\Phi _t\) in place of \(\varvec{W}_t\). Furthermore, it is more straightforward to define functions on the interval \([0,2\pi )\) compared to \(\mathcal {S}^{1}\), owing to the fact the former requires only a univariate periodic formulation to be valid. We therefore reformulate the non-stationary gauge function \(g_{\varvec{z}}(\cdot )\) via \(m(\phi ,\varvec{z}):=g_{\varvec{z}}(\varvec{v}(\phi ))\), \(\phi \in [0,2\pi )\); we opt for the notation \(m(\phi ,\varvec{z})\) over \(m_{\varvec{z}}(\phi )\) to emphasise the fact that we treat both \(\phi \) and \(\varvec{z}\) as predictor variables in our framework; see Section 3.3. Under this reformulation, we have \(\partial \mathcal {G}_{\varvec{z}} = \{ \varvec{v}(\phi )/m(\phi ,\varvec{z}): \phi \in [0,2\pi ) \}\); see Section 2.1.

To further simplify the formulation, we treat the discrete index t as a continuous covariate (i.e., we set \(\varvec{Z}_t = t\) for all t), and explicitly assume that any relevant or unobserved covariates vary smoothly as a function of t. Combined, these steps allow us to simplify the formulation of the framework introduced in Section 2 and henceforth, we use t and \(\Phi _t\) in place of \(\varvec{Z}_t\) and \(\varvec{W}_t\), respectively. We stress that our proposed methodology for this simplified representation is specifically adapted for modelling extremal dependence that smoothly varies in time; see also Section 6.

To approximate limit sets in the stationary setting, Wadsworth and Campbell (2024) study the tail behaviour of the variable \((R\mid \Phi = \phi )\). While the authors did not consider the non-stationary setting, their derivations extend to this case. In particular, if Assumption 2.3 holds, we have that \(f_{R_t\mid \Phi _t}(r\mid \phi ) \propto r\exp \{-rm(\phi ,t)[1+o(1)]\}\), \(r \rightarrow \infty \), for any \(t \in \{1,\hdots ,T \}\), where \(f_{R_t\mid \Phi _t}(\cdot \mid \phi )\) denotes the density function of \((R_t \mid \Phi _t = \phi )\). In addition, their arguments imply that for many dependence structures

giving asymptotic equivalence of Eq. 2 with a gamma kernel. This motivates the following modelling assumption:

where ‘truncGamma’ is shorthand for the truncated gamma distribution with shape and rate parameters 2 and \(m(\phi ,t)\), respectively, and \(r^\tau (\phi ,t)\) is the \(\tau \)-quantile of \((R_t \mid \Phi _t = \phi )\) for all \(t \in \{1,\hdots ,T \}\) and some \(\tau \in (0,1)\) close to 1, i.e., \(\Pr (R_t \le r^\tau (\phi ,t) \mid \Phi _t = \phi ) = \tau \). Selection of \(\tau \) is discussed in Section 4. In this context, we refer to observations satisfying \(\Vert \varvec{X}_t\Vert > r^\tau ({{\,\textrm{atan2}\,}}(X_{2,t},X_{1,t}),t)\) (i.e., exceeding the radial quantile) as joint tail observations.

For our framework, we assume that the assumption outlined in Eq. 3 holds exactly at some fixed \(\tau < 1\) for all t and \(\phi \in [0,2\pi )\). Empirical evidence from many approaches indicates that this assumption is flexible enough to capture a wide range of dependence structures (Wadsworth and Campbell 2024; Majumder et al. 2025; Murphy-Barltrop et al. 2024a). This resulting inference procedure involves two-steps: estimation of the non-stationary quantile function \(r^\tau (\phi ,t)\), and estimation of the non-stationary gauge function \(m(\phi ,t)\). We represent both of these functions using GAMs, an overview of which is provided in Section 3.3.

We remark that one can also leave the gamma shape parameter of Eq. 3 as a free parameter to estimate, as has been done in several existing limit set estimation approaches (e.g., Wadsworth and Campbell 2024, Murphy-Barltrop et al. 2024a). Fixing this quantity to 2 corresponds to assuming that the asymptotic form of Eq. 2 holds exactly at some finite level. In unreported results, we found that including the shape parameter in the non-stationary framework did not improve the quality of limit set estimates, and in some cases even reduced their quality. We therefore chose to fix the shape at 2, which has the added benefit of reducing parameter variability.

GAMs (Wood 2017) provide a flexible, semi-parametric framework for capturing complex functional forms without requiring rigid modelling assumptions. Typically, GAMs are used in practice to capture the relationship between a response variable and a set of predictor variables; the resulting fitted model can subsequently be used for prediction and interpolation across the joint domain of the predictor variables. Unlike classical regression techniques, which tend to assume simple linear relationships between variables, GAMs exploit smooth functional forms to capture non-linear patterns, alongside complex interactions that may arise between dependent predictor variables. Moreover, in the extreme value context, theoretical properties regarding convergence and asymptotic normality of GAMs have also been established (Yoshida 2026).

When fitting GAMs in practice, the optimal choice of smoothing functions and tuning parameters is context dependent, and fine tuning is typically required to ensure the model accurately represents the underlying data structure. Consequently, we only describe GAMs in the context of our proposed framework. As noted in Section 3.2, we wish to estimate the functions \(r^\tau (\phi ,t)\) and \(m(\phi ,t)\) for all t and \(\phi \in [0,2\pi )\). Thus, under the GAM framework, we view t and \(\phi \) as predictor variables. The interaction between these variables arises from the fact that for a fixed angle \(\phi \), the conditional radial tail behaviour at \(\phi \) will vary over time, while for a fixed time point t, the extremal dependence structure at t will vary across angles.

We propose the following GAM formulations

where \(\beta _0,\beta ^*_0 \in \mathbb {R}\) are intercept terms and \(s_{tp}(\cdot , \cdot ), s_{tp}^*(\cdot , \cdot )\) are smooth tensor product splines capturing the interaction between t and \(\phi \). Tensor splines are formed by taking products of univariate smooth spline functions defined on the predictor variables. Each univariate spline is formed of a piecewise combination of flexible polynomial functions, and the points at which these polynomials are connected are called knots. In our case, we use cubic and cyclic cubic splines for t and \(\phi \), respectively, where latter spline ensures the periodicity over the polar angular variable. Cubic splines posses many desirable properties, such as optimality in terms of smoothness, high flexibility and computational efficiency (Wood 2017). More complex modelling frameworks involving adding univariate smooth splines to the tensor splines of Eq. 4 did not lead to noticeable improvements within our modelling framework. An illustration of tensor product spline basis functions are provided in Appendix A.

Identical GAM formulations are taken for \(r^\tau (\phi ,t)\) and \(m(\phi ,t)\); this is due to the fact, as \(\tau \rightarrow 1\), these functions are approximately inversely proportional (Wadsworth and Campbell 2024) and thus will vary in a similar manner over the predictor variable domain(s). Let \(\varvec{\beta }\) and \(\varvec{\beta }^*\) denote the spline coefficients associated with the tensor product functions \(s_{tp}(\cdot , \cdot )\) and \(s_{tp}^*(\cdot , \cdot )\), respectively. Denote the number of knots (i.e., the basis dimensions of the univariate splines) as \(\kappa _t\) and \(\kappa _{\phi }\) for t and \(\phi \), respectively. Note that the basis dimensions correspond to the flexibility of the model, and their selection typically represents a trade-off; higher dimensions result in greater flexibility but increased parameter variability. In practice, the tensor product functions in Eq. 4 are also penalised to avoid over-fitting, which commonly occurs for models with large parameter sets. We return to this discussion in Section 4.4. The resulting parameter sets are given by \(\varvec{\theta } \!:=\! (\beta _0,\varvec{\beta })\) and \(\varvec{\theta }^* \!:=\! (\beta ^*_0,\varvec{\beta }^*)\), with tuning parameters \(\kappa _t\) and \(\kappa _{\phi }\) and joint predictor domain \(\mathcal {D} \!:=\! \{1,\hdots ,T \} \!\times \! [0,2\pi )\). We now consider estimation of \(\varvec{\theta }\) and \(\varvec{\theta }^*\).

Estimation of \(r^\tau (\phi ,t)\) represents a quantile regression problem, for which a wide range of approaches are available (Koenker et al. 2017). These techniques are typically based on the pinball loss function, which, when minimised with respect to a parameter set, returns quantile estimates for any probability level \(\tau \). Many recent works have incorporated GAM formulations in the context of quantile regression, allowing for flexible semi-parametric inference (e.g., Yu and Moyeed 2001, Oh et al. 2011, Koenker 2011, Youngman 2019). However, there are several issues with these approaches. First, the standard pinball function is piecewise linear and is consequently difficult to optimise, leading to difficulty when incorporating complex model formulations with large parameter vectors, such as Eq. 4. Furthermore, some existing approaches require one to manually select ‘smoothing’ parameters (i.e., parameters to mitigate overfitting for highly parametrised models), even though this selection is not straightforward and is crucial for balancing the trade-off between model fit and complexity. Moreover, to estimate GAM parameters, many approaches adopt an asymmetric Laplace model as a working assumption, even though such distributional assumptions are problematic and can result in poorly calibrated, inaccurate quantile estimates (Fasiolo et al. 2021a).

In this work, we employ the approach introduced in Fasiolo et al. (2021a), which overcomes the drawbacks of existing techniques. They introduce a novel technique for fitting GAM-based quantile regression models and demonstrate that the resulting estimator is well calibrated, stable and efficient. They also introduce a novel loss function which does not suffer from the optimisation issues faced by the pinball function. This novel loss function is penalised to impose smoothness on the resulting splines (i.e., avoid overfitting), resulting in penalty parameters. Penalty parameter estimates are obtained via a Bayesian approach, with a so-called ‘learning rate’ which determines the relative weights of the loss and penalty functions. Selection of basis dimensions is discussed in Section 4.4, and for every t and \(\phi \in [0,2\pi )\), we denote the corresponding quantile estimate by \(\hat{r}^\tau (\phi ,t)\).

With quantile estimates obtained over \(\mathcal {D}\), what remains is to estimate the non-stationary gauge function via the model in Eq. 3. For this, we employ the framework of Youngman (2019) to estimate the parameter vector \(\varvec{\theta }^*\). Taking the log-likelihood function associated with Eq. 3, this approach employs a technique known as restricted maximum likelihood (REML), whereby the log-likelihood function is penalised to avoid over-fitting, thus allowing for flexible GAM formulations in the likelihood function (Wood 2011; Wood et al. 2016). In a similar manner to Section 3.4, this results in penalty parameters that also need to be estimated as part of the inference procedure. Under the restricted maximum likelihood framework, this is achieved by making distributional assumptions on the smoothing parameters, then maximising the resulting likelihoods. For rigorous mathematical details and guidance on model fitting, see Wood (2017) and Youngman (2019).

Given t and \(\phi \in [0,2\pi )\), we denote the corresponding non-stationary gauge function estimate by \(\hat{m}(\phi ,t)\), with the corresponding non-stationary boundary set estimate given by \(\widehat{\partial \mathcal {G}}_t = \{ \varvec{v}(\phi )/\hat{m}(\phi ,t): \phi \in [0,2\pi ) \}\). Owing to the continuity of \(s_{tp}^*(\cdot , \cdot )\) (Eq. 4), the limit set associated with this estimate is always compact (Murphy-Barltrop et al. 2024a). However, there is no guarantee that \(\widehat{\partial \mathcal {G}}_t \in [-1,1]^2\), or that the boundary set estimate possesses componentwise maxima and minima \(\varvec{1}\) and \(-\varvec{1}\), respectively. See Sections 5 and 6 for further discussion.

Thus far, we have provided a general framework for inference on non-stationary boundary sets, but we have been intentionally ambiguous regarding specific features of the model. In particular, we have not yet specified a norm \(\Vert \cdot \Vert \), the quantile level \(\tau \), or the basis dimensions \(\kappa _t\) and \(\kappa _{\phi }\); these are important quantities for model inference, and as it turns out, the appropriate values vary depending on the form(s) of extremal dependence. Furthermore, one must also consider the effect of the observation window T and account for variability for any given sample size. Selection of these features is discussed in Sections 4 and 5.

To assess model fits, we introduce two diagnostic tools. Example applications of each diagnostic discussed in this section can be found in Section 5. Since we do not observe repeated observations for any time point t, one cannot compute scaled sample clouds and compare these to fixed boundary set estimates (see, e.g., Fig. 1), as is common for many stationary approaches. Therefore, a more nuanced approach is required for assessing model fits in the non-stationary setting.

Firstly, we adopt the goodness-of-fit metrics proposed by Wadsworth and Campbell (2024) and Campbell and Wadsworth (2024). In particular, we compute probability-probability (PP) and quantile-quantile (QQ) plots associated with the truncated gamma model fit. For the former, we produce two plots, one comparing observed and theoretical probabilities, and the other illustrating the differences between these probabilities over the observation rank. For the latter, all threshold exceeding observations are transformed to a standard exponential scale using Eq. 3. These transformed observations can be compared to theoretical quantiles. Combined, these metrics provide global diagnostics for model fits over both the polar angle and time variables. Confidence bounds can be added to these plots by exploiting the fact the k-th order statistic of a uniform distribution theoretically follows a \(\text {Beta}(k, n + 1 - k)\) distribution (David and Nagaraja 2004).

Next, we consider the return level set diagnostic proposed in Murphy-Barltrop et al. (2024a). Introduced in the stationary setting, this model computes the empirical probability of points inside an estimated return level set \(\hat{\mathcal {A}}^p_t\) for some p close to one, and compares this value to the nominal level. Repeating this procedure over a range of probability levels, one obtains a QQ plot representing the nominal and estimated probabilities. This diagnostic is easily extended to the non-stationary setting by averaging return level set probabilities over time. Assuming unbiased estimation of return level sets, we have that

Equation 5 can be approximated empirically via the technique described in Murphy-Barltrop et al. (2024a). In particular, given observations \(\varvec{x}_t, \; t=1,\hdots ,T\) with T large, we set \(\hat{p} = (1/T)\sum _{t=1}^T \mathbbm {1}( \Vert \varvec{x}_t \Vert \le \hat{r}^p( {{\,\textrm{atan2}\,}}(x_{2,t},x_{1,t}),t))\), where \(\mathbbm {1}(\cdot )\) is the indicator function and \(\hat{r}^p(\cdot ,\cdot )\) denotes the estimated p-quantile function of \((R_t \mid \Phi _t = \phi )\) (i.e., the radial values of \(\hat{\mathcal {A}}^p_t\)) derived from the truncated gamma model fit. The corresponding diagnostic describes how well the fitted model captures the observed structure in the data, and indicates the suitability of the model for inferring risk measures.

In this section, we explore the performance of the framework introduced in Sections 2 and 3. Section 4.1 introduces a range of simulated examples exhibiting non-stationary dependence. In Section 4.2, we discuss metrics for assessing model performance over \(\mathcal {D}\). Section 4.3 discusses estimation of the penalty parameters from the frameworks introduced in Sections 3.4 and 3.5. Model formulation is discussed in Section 4.4 and in Section 4.5, we present results for each of the simulated examples. To estimate the quantile and gauge functions, we employ the qgam (Fasiolo et al. 2021b) and evgam (Youngman 2022) software packages, respectively, in the R programming language. Example code for fitting our framework is available at https://github.com/callumbarltrop/NSGE, and wrapper functions are included into the GitHub version of evgam; see https://github.com/byoungman/evgam.

Given an observation window T, we define a range of non-stationary structures in terms of the time covariate t. All samples are simulated on standard Laplace margins. The first two examples are obtained using the bivariate Gaussian copula, for which dependence is characterised by the correlation coefficient \(\rho \in [-1,1]\). For the first case, we set \(\rho (t):= 0.2 + 0.6(t-1)/(T-1)\), giving \(\rho (1) = 0.2\) and \(\rho (T) = 0.8\), i.e., moving from weak to strong positive dependence. For the second example, we define the functions \(a(t):= 0.45(t-1)/(T-1)\) and \(b(t):= 0.5\sin (2.5\pi (t-1)/(T-1))\), and set \(\rho _2(t) = a(t) + b(t)\); this corresponds to a harmonic, increasing correlation function. The dependence trend is complex in this case, moving from positive, to negative, then back to positive dependence. An illustrative figure of \(\rho _2(t)\) over t is provided in Appendix A.

For the third example, we use the inverted bivariate extreme value copula with the logistic family (Ledford and Tawn 1997). Dependence is classified through the parameter \(\alpha \in (0,1)\), with positive dependence increasing as \(\alpha \) approaches 0. We set \(\alpha (t):= 0.3 + 0.4(t-1)/(T-1)\), resulting in a trend that moves from strong to weak positive dependence.

For the fourth example, we consider the student-t copula family, for which dependence is quantified via a parameter \(\sigma \in [-1,1]\) representing the degree of collinearity in the data, alongside the degrees of freedom parameter \(\nu > 0\). Lower values of \(\nu \) correspond with stronger tail dependence. In this case, we fix \(\sigma = 0.5\) for all t and set \(\nu (t):= 0.5 + 1.5(t-1)/(T-1)\); this corresponds to a trend whereby the collinearity parameter remains constant over time, but the tail dependence gradually decreases.

Finally, for the fifth example, we consider the copula model introduced by Huser and Wadsworth (2019). This model was first proposed on one-tailed margins, and consequently we reformulate it slightly to define it with standard Laplace margins. A version of this model with asymmetric tails has also been considered by Gong and Huser (2022). Let S denote a standard Laplace variable and \(\varvec{V}\) denote a bivariate random vector with Gaussian copula on standard Laplace margins with correlation coefficient \(\rho \), with S and \(\varvec{V}\) independent from each other. Given any \(\delta > 0\), define the random vector

where addition is applied componentwise. Note that \(\varvec{Y}\) does not posses standard Laplace margins, and additional marginal transformation involving the parameter \(\delta \) is required to obtain standardised margins. The model in Eq. 6 can be interpreted as follows: for \(\delta < 1\), \(\varvec{V}\) is heavier tailed than \(\delta S\), inducing AI, while for \(\delta \ge 1\), S dominates \(\varvec{V}/\delta \), resulting in AD. One can observe that a Gaussian copula is recovered as we let \(\delta \) approach 0. The transition between the two dependence classes occurs at \(\delta \!=\! 1\), corresponding to an interior point of the parameter space; this is in contrast to most extreme value models, for which this transition occurs at the boundary of the parameter space(s) (Wadsworth et al. 2017). Consequently, we fix \(\rho \!\!=\!\! 0.5\) and set \(\delta (t) \!:=\! 0.7 \!+\! 1.8(t \!-\! 1)/(T \!-\! 1)\), corresponding to a smooth transition between AI and AD.

We simulate 200 samples from each of the copula structures outlined in this section. This allows us to assess both the bias and variability of our proposed modelling framework across the different dependence structures. For each copula example, illustrations of the resulting boundary sets over time are provided in Fig. 2. We note that the cases where ‘pointiness’ is observed in some, or all, of the quadrant corners (e.g., the student-t) correspond to copulas possessing strong tail dependence within such regions.

The overall goal of our modelling framework is to estimate non-stationary boundary sets, and consequently we assess performance in terms of set estimates. However, owing to the additional time dimension, assessing performance is not straightforward in the non-stationary setting. For example, one could plot all two-dimensional boundary set estimates on a time axes alongside the true sets, but the corresponding three-dimensional plot would be difficult to interpret, and it may be unclear in which regions the model performs well. Furthermore, many performance metrics proposed in the stationary setting, such as mean integrated squared error of the gauge function over \([0,2\pi )\) (Murphy-Barltrop et al. 2025) or bias in extremal coefficients (Simpson and Tawn 2024a; Majumder et al. 2025), cannot be easily generalised to the non-stationary setting.

Consequently, we follow Murphy-Barltrop and Wadsworth (2024), who propose fixing either the time index or the polar angle, then evaluating performance visually over the domain of the other component. Firstly, given a set of fixed time points \(\{1, \lfloor T/2 \rfloor , T \}\), corresponding to the start, middle and end of the observation window, respectively, we compute estimates of the non-stationary boundary set at each \(t \in \{1, \lfloor T/2 \rfloor , T \}\) and plot these alongside the true boundary sets. This allows us to visually assess what boundary set shapes can be represented by the framework proposed in Section 3. Evaluating the joint tail at the boundaries of the observation window also provides a good indicator of the overall model performance, since there is less data available for regression at these points. Secondly, we fix angles \(\{\pi /4, 3\pi /4, 5\pi /4, 7\pi /4 \}\), corresponding to the main diagonals in all four quadrants. For each \(\phi \in \{\pi /4, 3\pi /4, 5\pi /4, 7\pi /4 \}\), we estimate the corresponding Cartesian point on the boundary set, denoted \(\hat{\varvec{x}}^{\phi }_t:= \varvec{v}(\phi )/\hat{m}(\phi ,t) \in \widehat{\partial \mathcal {G}}_t\), for every t. We then compute the Euclidean distance of \(\hat{\varvec{x}}^{\phi }_t\), i.e., \(\hat{r}^{\phi }_t:= \sqrt{ (\hat{x}^{\phi }_{1,t})^2 + (\hat{x}^{\phi }_{2,t})^2 }\), and plot this estimate against t, alongside the corresponding true values. This allows us to visually assess how well our framework is capturing the simulated dependence trends in different regions. Note that we evaluate the Euclidean distance of the Cartesian point to allow for comparison between different norm definitions of the radial component; plotting the estimates of \(\hat{m}(\phi ,t)\) for different norm definitions, for instance, would be meaningless, since the coordinates at which the gauge function is evaluated are dependent on the choice of norm.

We note that both steps of our modelling procedure, as outlined in Sections 3.4 and 3.5, require smoothing parameters for penalising loss functions used for optimisation. Such parameters can be estimated automatically by making distributional assumptions directly on the smoothing parameters (Wood 2003). Moreover, for the Bayesian quantile regression framework, a learning rate parameter must be estimated prior to sampling from the posteriors of the GAM parameters. Estimation of, or obtaining samples from, the smoothing parameters is computationally expensive. Moreover, across different data sets from the same copula and sample size, it is unlikely the same smoothing parameters will be estimated twice, implying one cannot easily compare GAM coefficient estimates for different copula samples. Therefore, prior to running our large simulation study, we opt to fix model formulations to speed up computation and ensure comparability across sampling iterations.

For each of the possible combinations of model tuning parameters and features (see Section 3.5), we fit our modelling framework to 10 unique samples from each copula and use the procedures outlined in Fasiolo et al. (2021a) and Youngman (2019) to automatically select the relevant parameters. For the Bayesian quantile regression procedure, we compute the median learning rate across the 10 samples and then fix the learning rate to this estimate for the remainder of the simulation study. Note that being in Bayesian framework, this still requires one to sample from the posteriors of the GAM parameter distributions, and there is no guarantee such distributions will be the same across samples. For the frequentist REML framework, we compute the median smoothing parameters for the tensor product spline and then fix the smoothing parameter for any given combination to these medians. This ensures the same penalised likelihood function is used for each setup.

As noted in Section 3.5, our proposed framework requires us to specify a range of tuning parameters and features. Trying to select all of these terms simultaneously would be very difficult, and it is unlikely that the same model formulation would be optimal across all simulated examples. Consequently, we select each model component separately, loosely assuming the optimal choices of each model feature can be selected independently of others.

For the norm function \(\Vert \cdot \Vert \), we consider three popular choices: \(\Vert \varvec{x} \Vert _1 = |x_1| + |x_2|\), \(\Vert \varvec{x} \Vert _2 = \sqrt{x^2_1 + x^2_2}\) and \(\Vert \varvec{x} \Vert _{\infty } = \max \{|x_1|,|x_2|\}\), corresponding to the \(L^1\), \(L^2\) and \(L^{\infty }\) norms, respectively. As observed in Murphy-Barltrop et al. (2025), the unit ball associated with each norm directly influences the shape of the estimated boundary set. To see this, consider Fig. 3: in this plot, we have defined a cyclic cubic spline in polar coordinates, and for each norm, the corresponding Cartesian coordinates have been computed. While the spline is smooth on the polar scale, the resulting shapes in Cartesian coordinates possess noticeable kinks for the \(L^1\) and \(L^{\infty }\) norms. This occurs due to the shapes of corresponding unit balls, which are also illustrated in right panel of Fig. 3. Consequently, certain norm choices may be more appropriate for scaled sample clouds which appear particularly ‘pointy’ within certain regions. We stress here that in practice, it is unlikely the same spline would be estimated across all three norms, and Fig. 3 is meant to merely illustrate the difference in shapes when splines are transformed to the Cartesian scale.

For the quantile level, we consider \(\tau \in \{0.5,0.6,0.7,0.8,0.9\}\). This choice represents a typical bias-variance trade-off that is often observed in extreme value theory (see, e.g., Coles 2001, de Carvalho et al. 2026), and for the framework proposed in Section 3, this equates to assessing the rate of convergence of Eq. 3 to the true asymptotic tail distribution of the conditional radii variables. We wish to select \(\tau \) as high as possible without compromising the robustness or reliability of the modelling framework.

For the time basis dimension, we consider \(\kappa _t \in \{5,10,15 \}\) and define spline knots at equally spaced points along the observation window, i.e., for \(j = 1,\hdots ,\kappa _t\), set \(k^t_j:= 1 + (T-1)(j-1)/(\kappa _t - 1)\). As noted in Wood (2017), the exact location of spline knots is not so important, so long as the basis dimension is sufficiently high enough to provide adequate flexibility. For the polar angular basis dimension, we consider \(\kappa _{\phi } \in \{9,17\}\) and again define knots at equally spaced points on the angular window, i.e., for \(j = 1,\hdots ,\kappa _{\phi }\), set \(k^{\phi }_j:= 2\pi (j-1)/(\kappa _{\phi } - 1)\). These choices are practically motivated by the fact such basis dimensions result in knots that include the axes and principal diagonals in \(\mathbb {R}^2\); these regions are often where the boundary set exhibits key features related to extremal dependence properties. As such, it is important that our modelling framework can capture the behaviour around these regions. An illustrative figure of the polar angular knot locations is provided in Appendix A, and we refer to Simpson and Tawn (2024b) for further discussion.

To begin, we consider the effect of the quantile level \(\tau \). For this, we first fix the basis dimensions to their maximal values, i.e., \(\kappa _t = 15\), \(\kappa _{\phi } = 17\), so as to provide maximal flexibility, and just consider the \(L^2\) norm. For each example in Section 4.1, we also set \(T = 25,000\) to provide a large, but not unreasonable, sample size. Note that smaller T values are subsequently considered once all tuning parameters have been selected. For each copula example, we then apply the framework introduced in Section 3 across 200 examples and compute the performance metrics discussed in Section 4.2. These results are given in Appendix C.1. Encouragingly, the model appears remarkably robust to the choice of \(\tau \), accurately capturing the range of dependence structures over \(\mathcal {D}\). Naturally there is higher variability at higher \(\tau \) levels, as one would expect, though we also observe slightly better coverage for such values. Consequently, we fix \(\tau = 0.8\) for the remainder of this article; this value appears to offer a reasonable trade-off between bias and variance.

Next, we assess the effect of the basis dimensions \(\kappa _t\) and \(\kappa _{\phi }\). For the former, we fix \(\kappa _{\phi } = 17\) and again restrict attention to the \(L^2\) norm, then vary \(\kappa _t \in \{ 5,10,15\}\). Using the performance metrics described in Section 4.2, the results are given in Appendix C.2. For most of the copula examples, there appears to be little difference among the basis dimensions. However, for the complex dependence trend exhibited by copula 2, \(\kappa _t = 5\) does not appear sufficiently flexible to capture the underlying trends. Maintaining the same setup with \(\kappa _t = 15\) and allowing \(\kappa _{\phi }\) to vary in \(\{ 9,17\}\), we evaluate the effect of basis dimension for the angular component, with the results again given in Appendix C.2. It is clear that the additional flexibility arising from setting \(\kappa _{\phi } = 17\) allows the framework to better capture the complex variety of boundary set shapes. Taking these combined results into account, we fix \(\kappa _t = 10\) and \(\kappa _{\phi } = 17\) for the remainder of the article.

Finally, we consider the impact of the norm choice on the modelling framework. For this, we restrict attention to the first and fourth copula examples, where the underlying boundary sets are smooth and pointy respectively. With all other components of the model specified, we vary the norm across the three choices (i.e., \(L^1\), \(L^2\) and \(L^{\infty }\)). These results are given in Appendix C.3. It is clear that the choice of norm does not appear to significantly alter the rate of convergence for the modelling framework described in Eq. 3, but unsurprisingly affects the types of Cartesian shapes that can be represented by the framework. Visually, it appears that the \(L^2\) and \(L^{\infty }\) norms would be preferable for the first and fourth examples, respectively. However, without knowledge of theoretical boundary sets, it is not clear how one should select the norm choice in practice; we refer to Section 6 for further discussion. Consequently, for the remainder of this simulation study, we restrict attention to the \(L^2\) norm, acknowledging that this will impose smoothness on the resulting boundary set estimates.

Using the selected tuning parameters and model features, we present our simulation results. Fixing \(T=25,000\), Fig. 4 illustrates the boundary set estimates, along with the theoretical boundary sets, for \(t=1\) and \(t = T\) with the second, fourth and fifth copula examples; the results for the remaining copulas and time points are given in Appendix C.4. As can be observed, the estimated boundary sets capture the structures of the theoretical sets at the start and end of the time frame.

Setting \(\phi =\pi /4\) and \(\phi = 3\pi /4\) and again considering the second, fourth and fifth copula examples, the estimated radii of the corresponding boundary set points over time are illustrated against the truth in Fig. 5, with the figures for the remaining angles and copulas given in Appendix C.4. The framework appears to capture the dependence trends in most cases, especially for the Gaussian copulas. However, the model can fail to capture specific dependence features: for example, the abrupt transition between AI and AD for the copula described by Eq. 6, and the ‘pointiness’ of the student-t copula in the corners. Given the limited amount of knowledge available generally within the non-stationary setting, estimating dependence features exactly is seldom possible, and we already know that the \(L^2\) norm will not be able to represent the pointedness that arises in certain quadrants. On the whole, however, the slight bias in radii observed at certain angles is secondary to the fact that the estimated Cartesian sets provide reasonable approximations of the true limit sets for each of the copula examples.

Finally, we vary the observation period \(T \in \{5,000, \; 10,000, \; 25,000 \}\) and assess the corresponding boundary set estimates. These results are given in Appendix C.5. Unsurprisingly, increasing the sample size improves bias and reduces variability. However, even for \(T=5,000\), the model still appears able to approximate the shapes of the true boundary sets, and thus is somewhat robust to lower sample sizes.

Overall, our proposed framework is able to accurately capture a wide range of extremal dependence trends, demonstrating the robustness and versatility of the model. Furthermore, the choices of tuning parameters and features generally appear to have little overall influence on the model fits, provided an adequate amount of flexibility is permitted.

We now apply our modelling framework to the NASDAQ stock market. Specifically, we consider data from Apple, Microsoft, Amazon, Google, and Nvidia; these firms are all part of the ‘Magnificent Seven’ (M7), a term which refers to a group of seven highly influential big tech companies. As of September 9th, 2025, these firms constitute the five largest companies (in terms of market capitalisation) in the world. We do not consider the remaining M7 companies (Meta and Tesla) due to limited observational periods of the corresponding time series. In recent years, M7 companies have significantly impacted the U.S. stock market, driving a substantial portion of its growth. Their dominance stems from their leadership in key sectors like consumer electronics, cloud computing, social media, artificial intelligence, and electric vehicles. See Peprah et al. (2025) and de Carvalho and Palacios Ramirez (2025) for further discussion on M7 and on extreme losses of ‘big tech’ stocks, respectively. Owing to the fact these companies began publicly training on different days, the observational dates, and hence sample size, vary depending on which stocks are being considered. The listing days range from December 12th, 1980 (Apple) to August 19th, 2004 (Google). Daily prices data was accessed from NASDAQ on August 25th, 2025, resulting in time series ending on August 22nd, 2025.

For each company, we take the daily opening price \(P_t\) and compute log-returns \(Q_t:= \log ( P_{t}/P_{t-1} )\), for all \(t > 1\). This quantity represents a mathematically convenient way of representing price changes in financial markets (Tsay 2010). To avoid computational issues, we remove any time indices t for which \(Q_t = 0\) for any stock; this corresponds to two-day periods where there was no change in the opening price of at least one stock. Including such values would result in mixed discrete-continuous time series, which would violate the assumptions of our modelling framework.

In the context of financial risk management, many stakeholders are interested in understanding the joint tails of log-returns across different assets. For example, one may wish to analyse whether an extreme increase or decrease in one stock price will also result in similar changes across other stocks. This motivates the use of multivariate extreme value theory to evaluate joint tail behaviour in regions where variables are simultaneously extreme—either in the lower or upper tails. Considering the bivariate setting, this context represents an application for which tail behaviour across all quadrants is of interest; hence, statistical techniques restricted to a single quadrant (i.e., where both variables are simultaneously large) offer less utility in this setting.

In Section 5.2, we describe our procedure for modelling the marginal distributions of each time series in order to transform to standard Laplace margins. We apply our modelling framework in Section 5.3; the diagnostics subsequently presented in Section 5.4 indicate reasonable model fits across all pairs. Finally, in Section 5.5, we summarise dependence trends over the observation period and introduce further practical use cases for our proposed framework.

To apply the framework introduced in Section 2.2, we first require marginally IID observations on standard Laplace margins. This requires careful consideration for the log-returns data in question. In particular, log-returns series are unlikely to be independent; see McNeil et al. (2015) for detailed discussion. By computing the autocorrelation function (ACF) of \(Q_t\) and \(Q^2_t\) at different lags for each stock, it is clear that our data exhibit non-negligible temporal dependence; see Appendix D for the corresponding plots.

To account for the observed temporal dependence, we model the log-returns using a GARCH(1,1) process and compute the standardised residuals, which we then take to be independent. For this, we assume \(Q_t = \mu + \sigma _t \varepsilon _t\), where \(\mu \in \mathbb {R}\) denotes a mean term, \(\varepsilon _t\) denotes the residual process (with zero mean and unit variance), and the conditional volatility \(\sigma _t^2\) satisfies

with parameters \(c > 0\), \(\alpha , \beta \ge 0\) and \(\alpha + \beta < 1\). After fitting the model, the standardised residuals \(\hat{\varepsilon }_t = (Q_t - \hat{\mu }) / \hat{\sigma }_t\), which we term the filtered log-returns, should be approximately independent and identically distributed if the GARCH model is well specified. This works because the GARCH filter removes serial dependence in the conditional variance, capturing volatility clustering and leaving residuals that behave more like white noise (Engle 1982; Bollerslev 1986).

The plotted time series of the filtered log-returns are given in Appendix D, alongside ACF plots for \(\hat{\varepsilon }_t\) and \(\hat{\varepsilon }^2_t\). The latter plots indicate that the standardised residuals are independent in time. Furthermore, the plotted series do not seem to suggest any obvious periodicity or long term trends are present, suggesting the data may exhibit marginal stationarity. We additionally tested this property by applying Augmented Dickey–Fuller (Dickey and Fuller 1979) and Kwiatkowski–Phillips–Schmidt–Shin (KPSS, Kwiatkowski et al. 1992) tests to each filtered time series. In almost all cases, both tests at a \(5\%\) significance level indicated a stationarity assumption was reasonable for the data; the only exception was the KPSS test on the Apple time series for the Apple–Google pairing. Given the same test indicated stationarity for Apple series in all other pairings (for which the observation dates in Apple–Google pairing are included), and the computed p-value was just below the \(\alpha =0.05\) significance level, it is therefore reasonable to assume that this is a type 1 error that can be ignored. We henceforth assume that the filtered log-returns are marginally IID.

As an additional check, we compute empirical bivariate extremograms via the extremogram R package (Frolova and Cribben 2016) for all pairs of squared residuals. These plots are illustrated in Appendix D and indicate vanishing joint tail dependence, suggesting one can assume that extremes between filtered series are approximately independent across different time lags. We refer to Davis and Mikosch (2009) for further discussion.

We now employ standard techniques to model the individual marginal distributions. In particular, we employ the semi-parametric approach of Coles and Tawn (1991), where a generalised Pareto distribution (GPD) is fitted to the tails while the body is modelled empirically. In the case of returns data, we are interested in understanding both the extreme high and low values. Letting \(\mathcal {E}\) denote a filtered log-returns series with observations \(\{\varepsilon _1,\hdots ,\varepsilon _T \}\) and given some \(\alpha \) close to 0, define l and h to be the empirical \(\alpha \) and \((1-\alpha )\) quantiles, respectively, of \(\mathcal {E}\). We use the following marginal model

where \((\sigma _1,\xi _1)\) and \((\sigma _2,\xi _2)\) are the estimated GPD parameters from the conditional tail variables \((\mathcal {E} - h) \mid \mathcal {E} > h\) and \((l - \mathcal {E}) \mid l > \mathcal {E}\), respectively, \(\{x \}_+:= \max \{0,x\}\), and \(\tilde{F}_\mathcal {E}\) is an empirical rank transform given by \(\tilde{F}_\mathcal {E}(\varepsilon ) = \sum _{t=1}^T \mathbbm {1}(\varepsilon _{t} \le \varepsilon )/(T+1)\).

The selection of \(\alpha \) is equivalent to selecting a threshold for the GPD, which represents a bias-variance trade-off in practice. For this, we initially applied a leading automated threshold selection technique (Murphy et al. 2025) which outperforms many existing approaches. However, we found the model fits from the Murphy et al. (2025) approach to be inadequate for several time series. Therefore, we instead manually tested a range of \(\alpha \) values and eventually fixed \(\alpha = 0.03\); this appeared to give reasonable model fits for all of the tails. The resulting GPD QQ plots are given in Appendix D.

With marginal models specified, the probability integral transform is applied to obtain each series on standard Laplace margins. As an additional check, the series are split into non-overlapping rolling windows, and maximum likelihood techniques are used to estimate the location and scale parameter of the Laplace distribution for each window. These rolling window estimates remain close to the true Laplace parameters (i.e., zero location, unit scale) in all cases, indicating the margins have been successfully standardised. These plots are given in Appendix D.

Time-varying joint modelling in stock markets has been recognised as fundamental by many works (Poon et al. 2003; Patton 2006; Salvatierra and Patton 2015; Castro-Camilo et al. 2018), and we now use the framework proposed in Section 3 to this end. Since our model is bivariate, there are a total of 10 pairs of variables we could consider. For the sake of brevity, we restrict our analysis to consider the extremal dependence between Apple and all of the remaining stocks. Apple has been listed for longer than any of the remaining companies, thus representing the most ‘mature’ company, and the stock movements of Apple are often seen as an indicator of broader market trends, particularly within the technology sector. The size of observation windows varies from \(T = 6,598\) for the Apple–Nvidia pair to \(T = 9,432\) for Apple–Microsoft. As demonstrated in Section 4, such sample sizes are sufficient for applying the modelling framework proposed in Section 3.

Figure 6 illustrates scatterplots of each bivariate data set on Laplace margins, processed as described in Section 5.2. We just show observations in the joint tail (i.e., observations exceeding the estimated threshold function). A colour gradient is used to illustrate the change in joint behaviour over time, suggesting subtle variations in dependence. Ideally, one would wish to formally test for this feature; however, there are no best practices established for such testing within the literature, and we found a previous technique, which involves computing dependence coefficients across overlapping rolling windows (Castro-Camilo et al. 2018; Murphy-Barltrop and Wadsworth 2024), to be unreliable and misleading. In particular, for simulated IID data, the rolling windows approach could suggest spurious trends when using point estimates together with pointwise confidence intervals. We refer to Section 6 for further discussion.

To apply the framework, the same tuning parameters and formulations as suggested in Section 4.4, combined with the \(L^2\) norm, are used for model fitting. Figure 7 illustrates the time-varying limit sets for each pair of stocks; the corresponding threshold quantile sets, i.e., \(\widehat{\mathcal {R}}^{\tau }_t = \{ \varvec{v}(\phi )\hat{r^{\tau }}(\phi ,t): \phi \in [0,2\pi ) \}\) are given in Appendix D. These estimates suggest a range of potentially complex, non-linear extremal dependence trends are present for each pair, with the type of trend varying over quadrants. For example, within the positive quadrant, the extremal dependence varies moderately for the Apple–Google pair but only very slightly for the remaining pairs. Our analysis also reveals the dynamics governing joint extreme losses (i.e., third quadrant). For example, the pairs Apple–Amazon and Apple–Nvidia, shown in Fig. 7, exhibit strengthening extremal dependence in the lower-left (third) quadrant, suggesting that joint losses between these stocks have become more frequent in recent periods. These instances illustrate the benefits of applying the geometric approach on Laplace margins: one can evaluate the dependence structure across all quadrants simultaneously, providing a more complete picture of trends in extremal dependence.

To assess the models fits, we apply the diagnostics introduced in Section 3.6. Firstly, PP and QQ plots for the truncated gamma model fits were computed with 95% confidence intervals; Fig. 8 illustrates the QQ plots for all pairs, while the two variants of PP plots are given in Appendix D. One can observe generally good agreement between the observed and theoretical quantiles in all cases, although there is some deviation in the tail for the Apple–Nvidia pairing and the bodies of the Apple–Amazon and Apple–Google pairs. Overall, these plots indicate the truncated gamma model assumption of Eq. 3 is reasonable for capturing the radial tails across all pairs.

The return level set diagnostic is also computed for each model fit. Specifically, we define a set of 200 equally spaced probabilities from 0.8 to 0.99, and evaluate the average empirical probabilities over time. Letting p and \(\hat{p}\) be the true and empirical probabilities, respectively, we compare the pairs \((-\log (1-p),-\log (1-\hat{p}))\); this allows us to evaluate performance for reasonably extreme return level set probabilities. The corresponding probability plots are given in Appendix D. In all cases, we observe very good agreement between the pairs, indicating the estimated return level sets (and model fits) correctly represent the structure in the data.

We now demonstrate further potential use cases for our framework. To begin, we summarise the trends in dependence via estimates of coefficients. As detailed in Appendix B, coefficients of tail dependence, which summarise the form of joint tail behaviour, can be defined separately for each quadrant. As noted in Section 2, these coefficients can be approximated directly from the estimated limit sets using the theoretical results introduced in Nolde (2014) and Nolde and Wadsworth (2022). For each pair, the time-varying \(\eta _t\) estimates in quadrants 1 and 3 (i.e., the joint upper and lower tails, respectively) are shown in Fig. 9, with the remaining quadrants considered in Appendix D. Such coefficients provide interpretable summaries of the extremal dependence behaviour over time for different stock market dynamics. For example, these estimates suggest that the probability of joint high or low days for Apple and Nvidia returns are increasing throughout the time period, whereas the same probabilities for Apple and Google vary non-linearly over the observation period. It is also clear that the estimated trends can be very slight in some cases, such as the Apple and Microsoft pairing.

Setting \(p = 0.999\), we apply the framework to obtain estimates of return level sets, denoted \(\widehat{\mathcal {A}}^p_t\), which can be used as practical risk measures to inform decision making; see Haselsteiner et al. (2019); Papastathopoulos et al. (2023), and Simpson and Tawn (2024b) for further discussion. We consider three recent dates of significance in the financial context; namely, the onset of the COVID 2020 pandemic (March 16th, 2020, Mazur et al. 2021), the beginning of the Russian invasion of Ukraine (February 24th, 2022, Lo et al. 2022), and the 2025 ‘Liberation Day’ tariffs (April 2nd, 2025, Ignatenko et al. 2025). While the onset of the pandemic cannot be expressed by a single day, stock markets experienced significant volatility in mid-March, and thus this time period is of interest for studying joint tail behaviour. Return level sets for each date were computed across all pairs. The estimated sets were then transformed back to the original scale of the log-returns data, i.e., inverting the modelling procedure described in Section 5.2 and effectively ‘de-filtering’ the data. The resulting sets are illustrated in Fig. 10. These plots make it clear that the behaviour of this joint risk measure varies noticeably not only between different pairings, but also over different observation dates. We note here that differences across observations dates arise from both fluctuations in the dependence structure and marginal behaviour. For example, we observe significantly more market volatility around the onset of COVID, which helps to explain the larger return level set magnitudes on this date. Since all considered dates are relatively close in proximity, we observe only subtle differences in the extremal dependence structure; this is illustrated in Appendix D, where we show the same return level set estimates on the standard Laplace scale.

Finally, we demonstrate how our framework can be used to simulate from the conditional variable \((\varvec{X}_t \mid R_t>r^\tau (\phi ,t))\), i.e., the random vector in regions where the radial component is relatively large. Given \(t \in \{1,\hdots ,T \}\), this involves two steps:

1. 1.

   Simulation from \((\Phi _t \mid R_t>r^\tau (\phi ,t))\).

2. 2.

   Simulation from \((R_t \mid \Phi _t = \phi , R_t>r^\tau (\phi ,t))\).

Given a simulation \(\phi \) from \((\Phi _t \mid R_t>r^\tau (\phi ,t))\), step 2) simply involves applying inverse transform sampling using Eq. 3. Therefore, all that remains is to specify a model for step 1) from which angular observations can be simulated. For this, we apply a non-parametric kernel density approach known as the Nadaraya-Watson estimator (De Gooijer and Zerom 2003). Here, we set \(\hat{f}_t(\phi ) = [\sum _{t^* \in \mathcal {T}^\tau } K^\phi _{h_1}(\phi - \phi _t) K^t_{h_2}(\Vert t^* -t \Vert )]/[\sum _{t \in \mathcal {T}^\tau }K^t_{h_2}(\Vert t^* -t \Vert )]\), where \(\mathcal {T}^{\tau }:= \{t: \Vert \varvec{x}_t\Vert > \hat{r}^\tau ({{\,\textrm{atan2}\,}}(x_{2,t},x_{1,t}),t) \}\), and \(K^\phi _{h_1}\) and \(K^t_{h_2}\) denote kernel functions with bandwidth parameters \(h_1>0\) and \(h_2 > 0\) respectively. We use a circular von-Mises kernel for \(K^\phi _{h_1}\) and a Gaussian kernel for \(K^t_{h_2}\); the former accounts for the periodicity in the angular component. Bandwidth parameters are selected manually by comparing the estimated density functions to local histograms defined on time windows. Numerical integration is then applied to simulate directly from \(\hat{f}_t(\cdot )\) for any t. See Holmes et al. (2012) and Castro-Camilo et al. (2018) for similar applications of this modelling technique.

For each time point \(t \in \{1, \hdots , T \}\), we can simulate from the entire joint tail of \((\varvec{X}_t \mid R_t>r^\tau (\phi ,t))\). This presents a wide range of potential uses, and we present one specific use case here. When managing portfolios, practitioners often wish to assess the (tail) relationships between different stocks (McNeil et al. 2015). Joint risk measures are often used to aid with this assessment. One popular measure is known as the conditional value-at-risk (CoVaR) (Girardi and Ergün 2013; Adrian and Brunnermeier 2016). Given some probability p close to 0, the downside CoVaR at level p is defined in the non-stationary, bivariate setting via

where \({\text {VaR}}_{t,p}(X_{1,t}):= F^{-1}_{X_{1,t}}(p)\) is termed the downside value-at-risk (VaR) of \(X_{1,t}\). One can view \({\text {CoVaR}}_{t,p}\) as the resulting stress on \(X_{2,t}\) given an extreme loss for \(X_{1,t}\) (Nolde and Zhou 2021). However, as noted in Section 5.1, we are not just interested in extreme losses within our analysis. Thus, we adapt Eq. 7 to also consider the joint upper tail and define the upside CoVaR at level p via

where \({\text {VaR}}^*_{t,p}(X_{1,t}):= F^{-1}_{X_{1,t}}(1-p)\) is termed the upside VaR of \(X_{1,t}\). Conversely, \({\text {CoVaR}}^*_{t,p}\) describes the effect on \(X_{2,t}\) given an extreme increase for \(X_{1,t}\). We note the downside and upside naming conventions follow from similar nomenclature in the literature (e.g., Reboredo et al. 2016).

Observe that simulation can be used to approximate both \({\text {CoVaR}}_{t,p}\) and \({\text {CoVaR}}^*_{t,p}\) by empirically estimating the probabilities denoted in Eqs. 7 and 8 and using a rootfinder to solve for equality. For each date, we simulate 100, 000 observations in the joint tail and estimate the CoVaR measures at \(p = 0.01\) relative to the Apple returns (i.e., \(X_{1,t}\) denoting the Apple series). Figure 11 illustrates the VaR and CoVaR estimates for each pair of stocks on the observation scale; the corresponding estimates on the Laplace scale are given in Appendix D. In these plots, we also illustrate recent US business cycles relevant to the returns data in question. Each cycle corresponds to a period of recession, beginning at the peak (P) and ending at the trough (T) (see, e.g., de Carvalho et al. 2012). One can observe that the extremes in the VaR estimates for \(X_{1,t}\) do not always correspond to the extremes in CoVaR, illustrating the time-varying nature of the extremal dependence within the data. This provides evidence that classical multivariate extreme value modelling approaches would not provide adequate flexibility for capturing the structure in this data.

Additionally, in Appendix D, we illustrate VaR estimates of \(X_{2,t}\) against CoVaR estimates on the observation scale. We observe clear correlation between these estimates, yet CoVaR always dominates, owing to the impact of the extremal dependence between the stocks.

Interestingly, we also observe in Fig. 11 that both joint extreme gains and losses appear possible during recession periods. However, it is not just recession periods where such fluctuations are possible; considering the two years before and after 2020 (i.e., the COVID pandemic) we also observe non-recession periods where spikes are present in both the upper and lower tails, as illustrated in Fig. 12. In some cases, certain phenomena are recognisable (e.g., the Russian invasion of Ukraine in February 2022), while other spikes cannot easily be identified with an event or recession period. This demonstrates the complex and volatile nature of evolving extremal dependence in stock markets, alongside the importance of developing flexible modelling approaches that can capture fluctuations in this dependence.

As a final remark, we note that we observe significant trends in the CoVaR estimates on both the Laplace and observed scales, even though the estimated limit sets suggest weak trends in dependence. This occurs due to the fact our simulation procedure involves models for both the angular and radial components, and we observe significant variation in the estimated angular density \(\hat{f}_t(\phi ), \phi \in [0,2\pi )\), over time. This emphasises the importance of modelling both components accurately in the geometric setting.

In this article, we have introduced a novel geometric modelling framework for approximating non-stationary extremal dependence structures. This framework addresses certain limitations from existing works and, crucially, allows for more general and flexible modelling within the non-stationary setting. Our proposed model can capture a wide range of dependence structures, and is not limited to either AD or AI data sets. Furthermore, while we have proposed our framework in the non-stationary setting, the techniques given in Section 3 could be adapted to obtain stationary limit set estimates, adding to a growing literature in this area. As a by-product, our framework can be used to generate random vectors in regions where the radial component is large, thereby contributing to the fast-emerging literature on generative models for multivariate extremes, which aim to simulate data across all joint tail regions; see, e.g., De Monte et al. (2025) or Wessel et al. (2025). The ability to generate, or simulate, extremal data in different tail regions is desirable beyond the financial setting; for example, in the environmental setting to assist with flood risk assessment and catastrophe modelling (Keef et al. 2013; Quinn et al. 2019).

We close the paper with some comments on open problems and future work. Firstly, as noted in Section 3.5, the boundary set estimates obtained from our approach are not required to satisfy all of the theoretical properties of boundary sets outlined in Section 1; one could in future work explore whether this rescaling could be directly incorporated into the model fitting procedure, as was the case in Murphy-Barltrop et al. (2024a) and Campbell and Wadsworth (2024).

The proposed methodology introduced in Section 3 is specifically adapted for the case where extremal dependence varies smoothly in time. However, our framework is not restricted to this setting since GAMs are highly flexible and, with appropriate formulations, can capture a wide range of covariate interactions, including discrete covariates and abrupt regime changes (e.g., changepoints). Moreover, our code is written using the evgam R package, meaning it can easily be adapted to capture alternative covariate interactions. Adjusting our framework to capture more general trends in extremal dependence represents an important area for future work.

Currently, the proposed approach can only be applied the bivariate setting. Flexibly modelling non-stationary extremal dependence in the general multivariate setting remains an open problem. In the general d-dimensional setting, one can still transform to polar angles; however, fitting GAMs over such angles requires more than cyclic basis splines to achieve the correct periodicity properties. Alternatively, one could try to define models directly on the \((d-1)\)-sphere, as in Papastathopoulos et al. (2023) and Murphy-Barltrop et al. (2024a); however, this requires frameworks that can handle the collinearity that exists between the corresponding angular variables, which is not the case for many regression techniques. We also note flexible and accurate modelling of the angular distribution in higher dimensional setups remains an open problem in the stationary setting (Kakampakou and Wadsworth 2025), although some progress has been made in recent works (e.g., De Monte et al. 2025, Wessel et al. 2025).

There is no guarantee that our proposed estimators for the non-stationary gauge functions or boundary sets will converge to the corresponding true values. While it would be desirable to have theoretical proofs of consistency, say, such results generally necessitate strict and unreasonable modelling assumptions, which themselves can be difficult to verify. We have therefore opted for a more practical treatment of our proposed estimators, noting that for real data sets, one can only evaluate model performance using diagnostic procedures, such as those presented in Sections 3.6 and 5.

One notable observation from this analysis is that the choice of norm (i.e., the selected coordinate system) greatly affects the shapes of boundary set that can be represented by the modelling framework. In particular, when the boundary set exhibits ‘pointiness’ in certain regions, it may be appropriate to select a norm for which the corresponding unit ball is also pointy in these regions. However, it remains unclear how one can select this norm in practice, and future work could therefore explore whether model selection techniques could help one to perform this selection in a robust and systematic manner.

As stated in Section 5.3, best practices for establishing the existence of time-varying extremal dependence within data are yet to be developed. Constructing tests with a reasonable statistical power is implicitly more complicated in the extreme setting, since one achieves better convergence to asymptotic models as one moves further into the tail, but this vastly reduces the available sample size and hence power. Developing robust and reliable techniques for diagnosing trends in extremal dependence represents an important line for future work.

Finally, an obvious omission from the case study of Section 5 is the lack of uncertainty quantification for the fitted models and computed statistics. Standard bootstrapping techniques are not applicable here, since resampling such data samples would not preserve a time-varying dependence structure. Uncertainty quantification is further complicated by the fact our modelling framework contains many components, each of which is fitted independently; the quantile function, the gauge function, and the angular model. Jointly quantifying the uncertainty arising from each of these components is non-trivial and would require careful consideration. Being able to address the significance of trends represents an important line for future research. Noting that the trends we found in Section 5 were very slight, we also fit a stationary version of our model. The diagnostics for this also appear satisfactory, suggesting that the observed trends may indeed not be statistically significant in this case. A starting point for considering uncertainty quantification could be to adapt established approaches for non-stationary, univariate time series, such as the local block bootstrap (e.g., Politis et al. 1999, Chapter 4).

Authors

1. Callum J. R. Murphy-Barltrop
2. Jennifer L. Wadsworth
3. Miguel de Carvalho
4. Benjamin D. Youngman

$$\begin{aligned} \Pr ( \varvec{o}\varvec{X} > u ) \sim L(e^u)\exp (-u/\eta _{\varvec{o}}) \; \; \text {as} \; u \rightarrow \infty , \end{aligned}$$

(9)

Murphy-Barltrop, C.J.R., Wadsworth, J.L., Carvalho, M.d. et al. Modelling non-stationary extremal dependence through a geometric approach. Extremes (2026). https://doi.org/10.1007/s10687-026-00543-w

Download citation

• Received: 26 September 2025

• Revised: 07 April 2026

• Accepted: 10 June 2026

• Published: 24 July 2026

• Version of record: 24 July 2026

• DOI: https://doi.org/10.1007/s10687-026-00543-w