covXtreme : MATLAB software for non-stationary penalised piecewise constant marginal and conditional extreme value models

The covXtreme software provides functionality for estimation of marginal and conditional extreme value models, non-stationary with respect to covariates, and environmental design contours. Generalised Pareto (GP) marginal models of peaks over threshold are estimated, using a piecewise-constant representation for the variation of GP threshold and scale parameters on the (potentially multidimensional) covariate domain of interest. The conditional variation of one or more associated variates, given a large value of a single conditioning variate, is described using the conditional extremes model of Heffernan and Tawn (2004), the slope term of which is also assumed to vary in a piecewise constant manner with covariates. Optimal smoothness of marginal and conditional extreme value model parameters with respect to covariates is estimated using cross-validated roughness-penalised maximum likelihood estimation. Uncertainties in model parameter estimates due to marginal and conditional extreme value threshold choice, and sample size, are quantified using a bootstrap resampling scheme. Estimates of environmental contours using various schemes, including the direct sampling approach of Huseby et al. 2013, are calculated by simulation or numerical integration under fitted models. The software was developed in MATLAB for metocean applications, but is applicable generally to multivariate samples of peaks over threshold. The software and case study data can be downloaded from GitHub, with an accompanying user guide.


Introduction
Reasonable estimation of the characteristics of extreme environments is essential to quantify natural hazards, and assure the reliability and safe operability of man-made infrastructure.For example, the extreme ocean environment can be thought of as a multivariate random process in space and time, describing interacting wind, wave, current, tide and surge phenomena.Accurate theoretical modelling of extremes from the ocean environment is problematic, because of the complexity of the numerical calculations involved.Practical risk assessment for marine and coastal structures and operations therefore typically requires the development of a statistical model from observational data for the joint behaviour of multiple meteorological-oceanographic ("metocean") variables.
Physical intuition and previous work has demonstrated that the marginal distribution of a variable such as ocean storm severity (quantified in terms of significant wave height, H S ) is dependent on covariates such as the direction in which the storm evolves, and the time of year in which the storm occurs (e.g.Randell et al. 2016).Over the long term, there is evidence that storm severity also evolves slowly in time due to anthropogenic climate change as well as natural climate cycles (e.g.Ewans and Jonathan 2023).It is important therefore that any statistical model for the extreme environment captures covariate non-stationarity.The nature and extent of covariate dependence varies across metocean variables.Experience also suggests that it is reasonable to assume that typical metocean variables are drawn from so-called "max-stable" distributions (e.g.Coles 2001).Therefore, typically conditional on covariates, the conditional distribution of threshold exceedances will be well approximated by the generalised Pareto (GP) distribution, provided that the threshold level is sufficiently high.This suggests that the marginal characteristics of metocean variables can be described reasonably using a non-stationary GP model (e.g.Davison and Smith 1990, Chavez-Demoulin and Davison 2005, Randell et al. 2016).
In addition we require that the statistical model also describes the joint tail of all metocean variables in general.However the nature of extremal dependence between different metocean variables is generally unknown.The specification of the statistical model therefore needs to be sufficiently general to admit different extents of extremal dependence, the specifics of which are then estimated by fitting the model to data.The conditional extremes model of Heffernan and Tawn (2004) is an attractive candidate, because it admits different classes of extremal dependence, it has a simple form and is relatively easily interpretable.There is also evidence that the nature and extent of extremal dependence also varies systematically with covariates (e.g.Jonathan et al. 2013, Ross et al. 2018, Shooter et al. 2021).
The sizes of data samples available (e.g. from direct observation of the environment) for modelling are typically small relative to the time-scales of the inference task.For example, in engineering design we typically need to characterise events which occur on average once in 1000 or 10000 years based on a sample of at most the order of 100 years of observation.In these circumstances, our estimate of the joint tail of the distribution of metocean variables is typically uncertain.It is critical that this uncertainty is captured and incorporated appropriately in our risk assessment.
For individual metocean variables, the size of a rare occurrence is often quantified in terms of a return value associated with some return period T (=1000 years, say).The return value is defined straightforwardly as the quantile of the distribution of the annual maximum of the metocean variable with exceedance probability 1/T , or (for large T , see Jonathan et al. 2021) the quantile of the distribution of the T -year maximum with non-exceedance probability exp(−1).In engineering design, the joint tail of the distribution of two (or sometimes three) metocean variables is often quantified in terms of an environmental design contour ; points on the contour are "equally rare" with respect to some criterion, related to the joint cumulative distribution function or density of the variables (e.g.Ross et al. 2020, Haselsteiner et al. 2021).Therefore estimates for the distribution of the annual maximum or that of the T -year maximum for individual metocean variables, and environmental contours for pairs or triplets of variables, are key outputs from the statistical analysis.
There are many software tools available for applied extreme value analysis; Stephenson and Gilleland (2006), Gilleland et al. (2013), Gilleland (2016) and Belzile et al. (2023) provide reviews.In the statistical software literature, the R packages texmex (Southworth et al. 2024), evgam (Youngman 2022 andmgcv (Wood 2023) provide useful functionality for extreme value analysis.The discussion of Belzile et al. (2023) notes the lack of good "off the shelf" software for practitioners, particularly incorporating more recent advances in methodology.Some software has been published for environmental contour estimation (e.g.Haselsteiner et al. 2019, Mackay andde Hauteclocque 2023).The purpose of the covXtreme software introduced here is to provide the environmental practitioner with a straightforward means to estimate marginal and joint tails of distributions of random variables, and quantify extremes in terms of return values and environmental design contours.covXtreme accommodates covariate non-stationarity in both marginal and dependence behaviour, provides flexibility in estimating extremal dependence, and careful uncertainty quantification in inferences.
The sophistication of the methodology underpinning covXtreme has been set deliberately, based on the authors' experience of metocean applications, to accommodate necessary effects in a pragmatic manner, whilst avoiding unnecessary complexity.Extreme value analysis is in some senses a less mature research area (e.g.Davison and Huser 2015) with numerous pressing challenges.For example, there are many approaches to modelling non-stationary marginal extremes; the work of Biller (2000), Wood (2004), Brezger and Lang (2006), Huser and Genton (2016), Wood et al. (2016), Wood et al. (2017), Youngman (2019), Youngman (2020), Shao et al. (2022) provide examples.With sufficiently rich data, complex extreme value models can be well estimated.However, in practice, typical samples tend only to support the adoption of relatively simple covariate models; for example, Zanini et al. (2020) show that a simple piecewise constant covariate representation (used by covXtreme ) provides marginal inferences competitive with those from more sophisticated P-spline (e.g.Eilers and Marx 2010) and Bayesian adaptive regression spline models (e.g.Biller 2000) for a metocean application.
Elements of the covXtreme methodology, such as marginal and dependence modelling for extreme values, uncertainty quantification and environmental risk decision support, are reflected in recent developments in the environmental software and modelling literature.These typically address modelling of extremes of rainfall, river levels, temperatures and heatwaves, but sometimes less common environmental variables such as urban black carbon (Maciejewska et al. 2015).For example, regarding marginal extreme value analysis, De Luca and Napolitano (2023) introduce the EXTRASTAR software for modelling of time-series of annual maxima, applied to rainfall.Zhao et al. (2022) use the generalised extreme value distribution to estimate return values of rainfall.Diez-Sierra and del Jesus (2017) use extreme value analysis of block maxima within their MENSEI-L software.Souaissi et al. (2023) exploit extreme value analysis of block maxima in their model for river temperatures.Cobos et al. (2022) offer extreme value functionality within their MarineTools.temporalpackage.Minguez et al. (2010) discuss non-stationary extreme value analysis of block maxima, using the Akaike Information Criterion (AIC) for optimal covariate parameterisation.Zhang et al. (2022) consider non-stationary extreme value analysis of block maxima for heat wave tracking.Regarding joint modelling, Cheng et al. (2023) consider both marginal extreme value analysis and dependence modelling using copulas in their development of a drought index, and Liu et al. (2021) consider vine copulas in their model for water level prediction.Hao et al. (2017) develop the R package drought incorporating marginal and joint extreme value modelling using copulas, for applications to modelling and assessment.Regarding uncertainty quantification, Villa et al. (2023) present a stochastic model for future extreme temperatures, for decision support in infrastructure analysis, and the ValueDecisions software of Haag et al. (2022) considers uncertainty quantification for decision making, including some elements of extreme value analysis.More generally, Čampulová et al. (2022) discuss extreme value analysis in the context of identifying outliers in environmental time-series, and Vrugt (2016) offers extreme value functionality in their DREAM software for Bayesian inference using Markov chain Monte Carlo.
An early version of the covXtreme software was published alongside Ross et al. (2018) for modelling of storm surge, and Ross et al. (2020) for estimation of environmental contours.Bore et al. (2019) suggest that covXtreme is useful for analysis of extreme ocean current profiles with depth.Vanem et al. (2020) use covXtreme to compare different contour methods with response-based methods for extreme ship response analysis.Guerrero et al. (2023) have used their enhancement of an earlier version of covXtreme for analysis of electrical signals in the human brain.The software has motivated the development of analogous prototype software PPL for marginal extreme value modelling with penalised piecewise-linear covariate representations (Barlow et al. 2023).The software has also been used as a pre-processor for transformation of data to standard marginal scale, allowing joint and conditional extreme value analysis (e.g.Shooter et al. 2021;Shooter et al. 2022), and in metocean consultancy work.
Target audience for covXtreme covXtreme was developed for oceanographic applications, but it can and has been used more widely.We believe that the reader will find covXtreme potentially useful, if (a) there is interest in quantifying extremes using statistical analysis of a data set, and possibly estimation of extreme quantiles or "return levels"; (b) the data are not stationary, and vary with known covariates; (c) there is interest in quantifying extremes of multiple variables at the same time; and (d) there is a need to quantify uncertainty carefully.Under these conditions, covXtreme offers a pragmatic but statistically sound analysis.

Objectives and layout of article
The objective of the current article is to provide motivation and description of the covXtreme software, and illustrations of its use in the development of design conditions for ocean engineering.The layout of the article is as follows.Section 2 provides an overview of the software and the statistical methodology on which it is based.Sections 3 and 4 present case studies, involving a bivariate response and single covariate (Section 3), and trivariate response with 2-D covariate (Section 4).An accompanying user guide (available at Towe et al. 2023b) provides a detailed step-by-step description of developing a covXtreme model for ocean engineering data sets provided with the software.
2. Overview of software and statistical methodology covXtreme provides simple software for multivariate non-stationary extreme value analysis of peaks over threshold.The modelling framework of covXtreme (a) is statistically straightforward to understand and use without sacrificing rigour, (b) is computationally efficient, (c) provides good estimation of key quantities of interest to the extreme value practitioner, and (d) provides realistic quantification of modelling uncertainties.
Extreme value analysis of peaks over threshold of metocean variables requires the satisfactory completion of a number of analysis steps.Inference typically involves using a sample of time-series data for a number of variables (corresponding to some period P years of observation) as the basis for estimating the characteristics of extreme environments.These might include return values, associated values (e.g.Towe et al. (2023a)) and design contours (e.g.Haselsteiner et al. 2021) corresponding to a return period of T years, T ≫ P ).For reasonable extreme value inference of a single variable (such as significant wave height), key features of the data must be accommodated carefully in the analysis, including the serial correlation of time-series, the dependence of values of metocean variables on covariates such as direction and season.Appropriate model forms for tails of distributions should also be adopted.For peaks over threshold analysis of a single variable, choice of threshold level can be source of considerable uncertainty.Moreover, different metocean variables (such as wind speed and significant wave height) are likely to be intercorrelated, and the nature of this dependence must be characterised carefully especially for extremes of one or more variables, using appropriate model forms; in these models, choice of threshold is again not always straightforward.Further, since samples of peaks over threshold are generally small (e.g. for high thresholds), it is critical to quantify uncertainty in return values and associated extremal quantities thoroughly using well-understood estimators.
The structure of the covXtreme software is written to reflect the sequential nature of practical extreme value analysis, as a set of MATLAB functions (see Section 2.6), illustrated in the workflow of Figure 1.For a given application, the five stages are executed sequentially, leading the user through the inference, addressing important analysis considerations in order.A typical analysis stage involves the specification of tuning parameters, the appropriate choice of which is informed by diagnostic information generated at that stage; generally it will therefore be necessary to repeat the analysis stage until satisfactory diagnostics are achieved, before moving on to the next analysis stage.
Given a sample of multivariate time-series of metocean variables, a typical analysis would proceed as described below.Note, for ease of reference, that subsection numbering 2.1-2.5 here corresponds to the numbering of stages 1-5 of the analysis procedure, and also to the numbering of sections in the covXtreme user guide.

Data preparation
The main aspect of data preparation is the isolation of storm peak events.Storm peak events correspond to peaks over threshold for a dominant variable (e.g.significant wave height for a wave-dominated metocean application), associated values of all other variables of interest (e.g.wind speed) per event, and storm peak covariate values for all relevant covariates (e.g.direction and season).Storm peak events are isolated using the procedure described by Ewans and Jonathan (2008), and discussed in Section 1 of the covXtreme user guide.Storm peak events are then assumed to be conditionally independent given (storm peak) covariates.
Mathematically, consider a sample Ÿ1 (t i ), Ÿ2 (t i ), ..., ŸD (t i ) of multivariate time-series for D metocean variables observed at regularly-spaced time points t i , i = 1, 2, 3, ... over some period.Data preparation requires the isolation of a sample of values { ẏi1 , ẏi2 , ..., ẏiD } N i=1 of N storm peak events (by convention for the first variable) and associated events (for the remaining variables) summarising the peak characteristics of each of N storms observed, and corresponding storm peak covariate values {x i1 , x i2 , ..., x iC } N i=1 for C storm peak covariates X 1 , X 2 , ..., X C defined on some domain X .(Note that the "double dot" notation (e.g.Ÿ ) indicates time-series variables from which storm peak events, indicated by "single dot" notation (e.g.Ẏ ) must be isolated.)Storm peak events are identified as local maxima (between successive up-and down-crossings of a given level) of Ÿ1 .Associated values for a storm are the values of the other random variables at the time of occurrence of the storm peak event.Storm peak and associated values are taken to be conditionally-independent given covariates, in the sense that ẏid can be viewed as an independent draw from Ẏd |(X Note that covXtreme also provides functionality to simulate data with known characteristics for checking of the performance of the statistical methodology.

Specification of covariate bins
covXtreme adopts a piecewise-constant parameterisation for the distribution of peaks over threshold on a partition of the covariate domain comprised of elements referred to as "covariate bins".Covariate bins are not estimated as part of the analysis, but must be specified carefully by the user prior to the analysis, for reasonable inference.Briefly, the user is required to partition the domain of each covariate (independently) into a set of intervals.The resulting covariate bins over all covariates are then simply Cartesian products of these intervals, and hence a rectangular grid of covariate bins.covXtreme assumes, within a given covariate bin, that the statistical properties of a storm peak or associated variable are no longer dependent on covariates.Diagnostic plots illustrating features such as the variation of storm peak and associated variables with individual covariates, are provided to inform reasonable choice of marginal partitions.Illustrations of the specification of covariate bins are given in Section 2 of the covXtreme user guide.
Mathematically, each observation ẏi1 , ẏi2 , ..., ẏiD is allocated to one of B covariate bins by means of an allocation vector A, with All observations within a specified covariate bin are assumed to have common extreme value characteristics, specified in terms of the parameters of the marginal model for peaks over threshold from that bin for each of the D components of the observation.Hence in particular, all threshold exceedances of storm peak variable Ẏd , d = 1, 2, ..., D from covariate bin b can be viewed as independent draws from a generalised Pareto distribution with common shape, scale and threshold parameter values.

Marginal modelling
A two-stage modelling procedure is used to describe the marginal distribution of storm peak and associated variables.First a three-parameter gamma distribution is fitted independently to all the data for each covariate bin in turn, providing a good description of the bulk of the distribution of storm peak and associated variables on the covariate domain.The extreme value threshold for each covariate bin is then set to the quantile of the corresponding gamma distribution with pre-specified non-exceedance probability.Finally the distribution of threshold exceedances per covariate bin is estimated using a non-stationary generalised Pareto distribution.Since generalised Pareto shape parameter is typically more difficult to estimate than scale, its value is assumed fixed but unknown across all covariate bins.Moreover, variation of the generalised Pareto scale parameter on the covariate domain is regulated using roughness penalisation, set to maximise the predictive performance of the generalised Pareto model on a hold-out sample within a cross-validation scheme.We judge from prior experience that these assumptions regarding the variation of generalised Pareto parameters are reasonable, and confirm using diagnostics during analysis that the assumptions give reasonably well-fitting models.By reducing the number of degrees of freedom for model fitting, the assumptions are a possible source of estimation bias, but also of a reduction in variance in estimated model parameters and return values.Further details and illustrations of marginal modelling are given in Section 3 of the covXtreme user guide.
Mathematically, for each variable Ẏd , d = 1, 2, ..., D, and covariate bin B = b, b = 1, 2, ..., B independently, we estimate a three-parameter gamma distribution using maximum likelihood estimation with sample likelihood , where f Gmm is the density of the gamma distribution with shape ω bd ∈ R, scale κ bd > 0 and location l bd ∈ R given by where Γ(•) is the gamma function.In practice for fitting the gamma model, the location parameters {l bd } are first estimated using a low empirical quantile of the sample per variate per covariate bin, and the remaining shape and scale parameters estimated using maximum likelihood.Next we calculate the extreme value threshold as ψ bd = F −1 Gmm (τ d ; ωbd , κbd , lbd ) using the estimated gamma parameters, for the pre-specified non-exceedance probability τ d , where F Gmm is the cumulative distribution function of the gamma distribution.The marginal sample likelihood for threshold exceedances of ψ bd for variable Ẏd over all covariate bins can therefore be written where f GP is the density of the generalised Pareto distribution with shape ξ d ∈ R and scale ν bd > 0 given by where [A] + = A when A > 0, and = 0 otherwise.The corresponding cumulative distribution functions of the gamma and generalised Pareto distributions are where γ(•, •) is the lower incomplete gamma function, and Optimal predictive performance for the marginal generalised Pareto model is achieved using roughnesspenalisation for the scale parameters {ν bd } B b=1 across covariate bins.The corresponding penalised negative log likelihood takes the form The smoothness penalty λ d controls the extent to which the generalised Pareto scale varies across covariate bins, and the penalty is proportional to the variance of ν over covariate bins.Parameters can be estimated to minimise the penalised negative log likelihood for each variable Ẏd , d = 1, 2, ..., D in turn for each of a set of pre-specified values for λ d .The optimal value λ • d of λ d is chosen to maximise predictive likelihood for a hold-out sample within a k-fold cross-validation procedure.Gamma parameter estimates per covariate bin, and generalised Pareto parameter estimates evaluated using the full sample for λ d = λ • d are carried forward to subsequent inference.
Since the numbers of parameters in the various marginal models above are relatively small, a simplex search procedure provides a straightforward approach to parameter estimation by minimisation of (penalised) negative log likelihoods.
Threshold selection for extreme value analysis of peaks over threshold is an important consideration (e.g.Northrop and Jonathan 2011, Scarrott and MacDonald 2012, Wadsworth 2016).Within a covXtreme analysis, multiple marginal models based on different random choices of threshold non-exceedance probabilities τ d ∈ I τ d ⊆ [0, 1) are constructed, d = 1, 2, ..., D. The user's task is to specify the interval I τ d for each variable Ẏd .As explained in Section 3 of the user guide, this choice is aided by numerous diagnostic plots, including examination of the stability of the estimated value of ξ d as a function of τ d .
Uncertainty quantification for marginal inference is performed using a non-parametric bootstrap procedure.The original sample of storm peak and associated values is resampled with replacement.For each variable Ẏd , d = 1, 2, ..., D, the full marginal extreme value analysis is then repeated using the bootstrap resample together with a new random selection of τ d .The outcome of the complete marginal inference is then quantified in terms of sets of parameter estimates ( {ξ r d , {ν r bd , ψ r bd , ω r bd , κ r bd , l r bd ) for each of R (typically = 100 or 250) bootstrap resamples (where superscript r indicates an estimate from a resample).Typically, the value of λ • d estimated using the original sample is adopted for all bootstrap resamples also, although the option to estimate a new optimal roughness coefficient per bootstrap resample is provided.
The distribution of the T -year maximum event per covariate bin, and over all covariate bins can then be estimated (or sampled) under the fitted marginal model.From above, the full marginal model for storm peak or associated variate Ẏd in covariate bin B = b is Then under the model, the distribution of a random occurrence of Ẏd from any covariate bin is where p b is an empirical estimate of the probability of observing a storm event in covariate bin b.If we further assume that the number N of storms in T -years is Poisson-distributed with mean T ρ, where ρ is an empirical estimate for the number of storms per annum, and suppressing parameter dependence for brevity, the distribution of the T -year maximum is simply We can use a similar approach to estimate and sample from the distribution of the T -year maximum of Ẏd for any combination of covariate bins, by restricting the set of covariate bins over which the summation is performed in Equation 9 (and linearly scaling the values of {p b } so that they sum to unity).

Extremal dependence modelling
Given non-stationary marginal models for storm peak and associated variables, we next seek to describe the nature of the dependence between them for large values of the storm peak variable.This is achieved using the conditional extremes model of Heffernan and Tawn (2004).Under the fitted conditional extremes model, we can then estimate the characteristics of the joint distribution of all associated variables given a large storm peak, and thereby estimate joint environmental design conditions and design contours.The conditional extremes model is specified for sets of variables on a common standard Laplace marginal scale, rather than their original physical scales.For this reason, a necessary first step is to transform the sample of storm peak and associated values to this scale, using the fitted marginal models.Incorporation of covariate effects is generally important for extremal dependence modelling, and these can be accommodated using conditional extremes models of different complexities with respect to covariates.The covXtreme software allows any number of the parameters of the conditional extremes model to vary with covariates.Nevertheless, experience suggests that the data indicate the need for covariate non-stationarity of just one model parameter (the "slope" parameter, see below) for typical metocean applications.We therefore describe this specific model form as a recommended "default" approach here.Further details and illustrations of extremal dependence modelling are given in Section 4 of the covXtreme user guide.

Marginal transformation to standard Laplace scale
The marginal transformation to standard Laplace scale is achieved using the probability integral transform such that for i = 1, 2, ..., N , d = 1, 2, ..., D, b = 1, 2, ..., B F GmmGP ( ẏid ; ω bd , κ bd , l bd , ξ d , ν bd , ψ bd ) = F Lpl (y id ) for where F GmmGP is the marginal cumulative distribution function of storm peak variable Ẏd , for sets of parameters of marginal gamma and generalised Pareto models from Equation 8, and F Lpl is the cumulative distribution function of the standard Laplace distribution given by F Lpl (y) = 0.5 exp(−|y|) for y ≤ 0 and = 1 − 0.5 exp(−|y|) otherwise.

Conditional extremes modelling
The Laplace-scale sample {y i1 , y i2 , ..., y iD } N i=1 from random variables {Y 1 , Y 2 , ..., Y D } is next characterised using the conditional extremes model, for values y of the conditioning variate Y 1 above a dependence threshold ϕ(τ ) = F −1 Lpl (τ ), for which the conditional extremes model is assumed to hold, for carefully specified non-exceedance probability τ .For y > ϕ(τ ), in covariate bin indexed B = b, the recommended non-stationary model takes the form with linear slope parameters α bd ′ ∈ [−1, 1], d ′ = 2, 3, ..., D varying across covariate bins, scalar exponent parameters ] common to all covariate bins, and residual random variable For estimation of slope and exponent parameters only, we assume that each component Z d ′ of Z is independently distributed according to for mean µ d ′ ∈ R, scale σ d ′ > 0 and random variable W d ′ ∈ R following a generalised Gaussian (or delta-Laplace) distribution with zero mean, unit variance and shape δ d ′ .µ d ′ and scale σ d ′ are treated as nuisance parameters in the estimation, and are not needed for subsequent inferences under the fitted model.For general mean m and variance s 2 , the corresponding density f GG of the generalised Gaussian distribution is where The conditional dependence likelihood Ld ′ for the model Y d ′ |(Y 1 = y) and y > ϕ(τ ) can then be written As for marginal models, we regulate the smoothness of {α bd ′ } B b=1 optimally on the covariate domain using a crossvalidation procedure, selecting a value λ• d ′ for roughness coefficient λd ′ in roughness-penalised negative log likelihood to maximise predictive performance on cross-validation hold-out samples.Further, it is important to confirm reasonable choice of the dependence threshold τ by inspection of diagnostic plots, for example of stability of estimated parameters as a function of threshold.Two algorithms are provided to minimise the penalised negative log likelihood in Equation 16.The default is a Newton-Raphson approach exploiting gradients; the alternative is a simplex search procedure.
Once parameter estimates are available, we represent the distribution of the (D − 1)-dimensional residual random variable Z using the set E b = {e i ′ bd ′ }, for all d ′ = 2, 3, ..., D, and all i ′ = 1, 2, ..., N such that y i ′ 1 > ϕ(τ ), b = A(i ′ ), of residuals from the fit per covariate bin, with elements During subsequent simulation under the fitted conditional extremes model, these residuals are resampled jointly as {e i ′ b2 , e i ′ b3 , ..., e i ′ bD }, thereby preserving the dependence between them, and hence the dependence between variables Y 2 , Y 3 , ..., Y D |(Y 1 = y, y > ϕ(τ ), B = b) per covariate bin.Similarly, simulation below threshold ϕ(τ ) is achieved simply by resampling from the original Laplace-scale storm peak sample.We note that covXtreme also facilitates estimation of variants of this dependence model, for which any number of conditional extremes model parameters α, β, µ and σ are allowed to vary between covariate bins, and their overall roughness penalised for good performance by extending the penalty term in Equation 16to include all "non-stationary" parameters; see Section 4.2 of the covXtreme user guide.It is also possible to pool estimates of residuals across covariate bins, useful when covariate bins with low occupancy are present.
Uncertainty quantification for extremal dependence inference is performed by extending the bootstrap procedure described in Section 2.3.Conditional extremes models are estimated for each of the bootstrap resamples based on which marginal models were estimated.The combined marginal and conditional extremes analysis therefore produces the set {{ξ r d , {ν r bd , ψ r bd , ω r bd , κ r bd , l r bd } of parameter estimates and residuals for each of R bootstrap resamples indexed by superscript r.As for marginal modelling, the value of roughness coefficient λ• d ′ is typically estimated using the original sample only, and adopted for all bootstrap resamples also, although the option to estimate a new optimal roughness coefficient per bootstrap resample is provided.
A sampling procedure is used to estimate the conditional return value distribution of the associated value for Ẏd ′ given a T -year maximum event for Ẏ1 per covariate bin, and over unions of covariate bins.The approach combines importance sampling from the marginal T -year maximum distribution of Ẏ1 (exploiting Equation 10) with Monte Carlo sampling under the fitted conditional extremes model, managing transformations between marginal and standard Laplace scales.Further discussion of conditional return values can be found in Towe et al. (2023a).

Estimation of design contours
Using the estimated marginal and conditional extremes models, covXtreme allows estimation of quantities such as return values for the dominant variable, conditional return values for associated variables, and environmental design contours.Haselsteiner et al. (2017), Ross et al. (2020) and Haselsteiner et al. (2021) provide recent reviews of contour estimation, encompassing a range of approaches.Algorithms to estimate three types of design contour are implemented in covXtreme .These are (a) constant exceedance contours; (b) direct sampling contours; and (c) conditional extremes (or "HT", in acknowledgement of the authors of the conditional extremes model, Heffernan and Tawn) constant density contours.Each of the three contour methods estimates a line (in 2-D) or surface (in 3-D) on which certain characteristics of the distributions of Ẏ2 |( Ẏ1 = y) (in 2-D, or ( Ẏ2 , Ẏ3 )|( Ẏ1 = y) in 3-D) for large y are preserved.For example, the constant exceedance contour in 2-D is a line on which the probability of exceedance in an "outward" sense (e.g. on the contour; see Jonathan et al. (2014) for details.The direct sampling contour of Huseby et al. (2015) is similar, except that now the probability in the half plane defined by the tangent to the contour at any point of interest is constant, and the contour itself must be convex.In 2-D, the conditional extremes constant density contour defines a line on which the joint density of the pair ( Ẏ1 , Ẏ2 ) is constant.By construction, all contours pass through a so-called "lock point", defined as an extreme quantile (e.g. the T -year return value) of dominant variable Y 1 and the corresponding conditional median of associated variate Y 2 |(Y 1 = T -year value).The lock point defines the value of the distributional characteristic preserved on the contour.Further discussion and illustrations are provided in Section 5 of the covXtreme user guide.The contour methods discussed in this section are intended to illustrate the kinds of inferences which can be made routinely using simulation under a fitted covXtreme model.All the approaches discussed provide means to quantify the extent of some or all of a cloud of realisations.No claims are made regarding their usefulness in any particular field of application (e.g.Mackay andHaselsteiner 2021, Hafver et al. 2022).The user interested in making specific inferences for particular structure variables, or choices of design contour, is free to implement these as alternatives or extensions to the existing Stage 5 of the analysis procedure Mathematically, in 2-D, the constant exceedance contour ) can be defined by the equation where 2 ) is a reference location (see e.g.Jonathan et al. 2014), and , 2, for some small probability p > 0. We set the value of r o 2 equal to the conditional mean Ẏ2 |( Ẏ1 =T-year value), and the value of r o 1 to the minimum value of Ẏ1 in the sample, to ensure continuity of the T -year contour when Ẏ1 exceeds the sample minimum.
To estimate the direct sampling contour in two-dimensions for probability level α, following Huseby et al. (2015), we first find the function C(θ), the (1 − α)-quantile of the distribution of the projection Ẏ1 cos(θ) + Ẏ2 sin(θ) for each value of θ ∈ [0, 2π) Then we estimate the contour and potentially further smooth C as a function of θ.Following Winterstein et al. (1993), for at T -year return period, it is recommended that the value of α be set to 1/T .The conditional extremes constant density contour consists of the set of points ζ ∈ R 2 for which the joint density f Ẏ (ζ) = c for some small value c > 0. Algorithms for estimation of the three contour types exploit importance sampling (e.g.Section 3.4 of Towe et al. 2021) for computational efficiency where possible.Contours can be estimated per covariate bin (e.g. as "directional" or "seasonal" contours), or integrated over covariate bins (e.g. to provide "omni" estimates over all bins).
The direct sampling contour is convex by definition, in the sense that it encloses a convex region (of R 2 for 2-D contours), possibly with the help of the coordinate axes.Depending on the nature of the joint density, the conditional extremes constant density "contour" could exist in the form of a set of disjoint curves, some of which may be closed.The constant exceedance and direct sampling contours are invariant to transformations of variables, whereas the conditional extremes constant density contour is not (e.g.Ross et al. 2020).

Accessing the software, and performing a typical analysis
The covXtreme is available for download from GitHub at Towe et al. (2023b).Software is written using MAT-LAB object-oriented programming.Full specifications of classes, properties and methods are therefore provided, and available for user inspection.
A typical analysis involves executing each of the MATLAB scripts Stage1 PeakPicking, Stage2 SetBinEdges, Stage3 FitMargin (for each variable marginally), Stage4 FitHeffernanTawn and Stage5 Contour in sequence, as illustrated by the workflow in Figure 1.As explained in detail in the covXtreme user guide, each stage of the analysis involves the specification of control parameters for that stage.Default values for control parameters are provided, but it is necessary for the user at each stage to assess whether values are appropriate by confirming that diagnostic plots generated have reasonable characteristics.It may be necessary to repeat a given stage multiple times until acceptable diagnostic characteristics are found, before proceeding to the next stage.
Sections 3 and 4 below provide brief descriptions of the application of covXtreme to a pair of random variables (specifically significant wave height and spectral peak period) varying with directional covariate (Section 3) and to four random variables (significant wave height, spectral peak period, wind speed and overturning momemt) varying with 2-D directional-seasonal covariate.

Case study : Single covariate, bivariate response
This section outlines to application of covXtreme to estimation of design contours for extreme storm peak significant wave height and associated spectral peak period, both varying with storm direction.The analysis follows the five stages discussed in Section 2. The data correspond to approximately 35 years of time-series output from a hindcast simulator for a location in the northern North Sea.
Following isolation of storm peak significant wave height H S , associated spectral peak period T P and storm peak direction using Stage 1, Figure 2 illustrates the output of Stage 2 of the analysis, showing storm peak H S and associated T P as a function of the direction from which storms emanate, in degrees clockwise from north.The figure also shows the user-input bin edges at 0 • , 20 • , 60 • , 225 • , 270 • and 315 • , specified so that the variation of H S and T P is approximately independent of direction within each bin, but different between bins.Thus, for example, the bin including 45 • corresponds to the land shadow of Norway, with resulting low values of storm peak H S , whereas the bins including 250 • and 340 • contain large storms from the Atlantic Ocean and Norwegian Sea. Figure 3  plots of associated T P on storm peak H S for each of the six directional bins.The rates of occurrence of storms, and the marginal characteristics of H S and T P are clearly different between bins.In directional bin [225,270), the dependence between T P and H S appears to be particularly strong.

provides scatter
Figure 4 illustrates marginal extreme value models for storm peak H S , the panels describing the variation of the estimates of generalised Pareto scale ν, gamma shape ω and scale κ with direction, in terms of means (solid) and 95% uncertainty bands (dashed) over bootstrap resamples and random choices of marginal threshold non-exceedance probabilities drawn from the interval [0.7, 0.9].The empirical density of corresponding (stationary) estimates of generalised Pareto shape ξ is shown in Figure 9 of the user guide to be Gaussian-like with mean at approximately -0.2.The variation of ν and κ with direction appears consistent with expectations given Figures 2 and 3.In particular, given estimated ν, the largest return values for storm peak H S would be expected to emanate from the Atlantic and Norwegian Sea.Indeed, inspection of directional maxima for storm peak H S for return periods of 10 and 100 years in Figure 5 confirms this: the largest contributors to the omni-directional maximum (in black) are the covariate bins corresponding to the Norwegian Sea ([315, 0), magenta) and Atlantic ([225, 270), cyan).
It is critical to assess the diagnostic plots generated by covXtreme to confirm that model fit is adequate.A number of illustrative diagnostic plots corresponding to this case study are shown in the user guide.For the current application, a plot of the predictive negative log likelihood from the cross-validation procedure (see Equation 7) suggests that the optimal choice λ • of roughness penalty lies at around three.The goodness of fit of the marginal model is assessed by examining the stability of the estimate for generalised Pareto shape parameter ξ as a function of extreme value threshold τ .Comparison of empirical tails directly from the sample, with corresponding tails (and their uncertainty) estimates under the extreme value model, per covariate bin and omni-directionally, also suggests the marginal model   is reasonable.The corresponding full marginal analysis must also be performed for associated T P with direction.
Using the marginal models, Figure 6 illustrates parameter estimates for the conditional extremes model for associated T P given storm peak H S on standard Laplace margins.The non-stationary estimate for slope parameter α suggests that dependence between the variables is high (near the maximum possible of unity), especially in the covariate bin corresponding to Atlantic storms.The uncertainty band for the slope parameter included unity, admissible at sub-asymptotic levels (Tendijck et al. 2023).The value of exponent parameter β is slightly larger than zero, indicating that the sizes of residuals "y β Z" from the conditional extremes model fit (Equation 12) grow very slowly with the conditioning value y on Laplace scale.The empirical densities for residual parameters µ and σ are typical in our experience; unimodal densities, with approximately Gaussian shape.
Diagnostic plots are again assessed to confirm reasonable goodness of fit for the conditional extremes model.For the current application, parameter estimates from α per covariate bin are reasonably stable as a function of dependence threshold τ on the interval [0.7, 0.85] (see Figure 20 of the user guide).Further, the distribution of residuals E from the model fit do not appear to be obviously dependent on the directional covariate (see Figure 18 of the user guide).The overall distribution of residuals from the model fit (see Figure 17 of the user guide) is typical; a Gaussian-like density with positive skew.
We simulate under the fitted models to generate the environmental design contours (see Section 2.5) shown in Figure 7 omni-directionally and Figure 8 per directional bin.The constant exceedance, direct sampling and HT density contours are labelled as "Exc", "Hus" and "HTDns" in the figures.The three different contour methods produce estimates which have similar characteristics in terms of describing the "extent" of the data cloud, reflecting the positive dependence between H S and T P , and passing through the appropriate lock points (shown in green in the figures).However the methods also use different definitions of the environmental contour; it is not surprising therefore that the contours estimates do not agree fully.The HTDns contour in particular produces a more variable estimate, especially when sample size is small (e.g. in the [20,60) directional bin in Figure 8).For engineering design, points on the contour with large values of H S would typically be adopted to test the integrity of models for the offshore or coastal structure of interest.

Case study : 2-D covariate, multivariate response
The second case study is an extension of that reported in Section 3, which incorporates the over-turning moment (OTM) experienced by an offshore structure subjected to wind and wave loading.Specifically, we seek to characterise the joint distribution of (storm peak) H S and wind speed (WS) conditional on extreme values of OTM, subject to directional and seasonal covariate effects.The analysis again follows the 5 stages described in Section 2.
Stage 1 of the analysis is isolation of storm peak events from the underlying storm H S time-series data, following the same procedure as in Section 3.For Stage 2, Figure 9 shows the directional and seasonal variation of the conditioning variate OTM, and associated variates H S and WS together with the directional and seasonal bin edges specified.Note that a smaller number of directional bins is used for this case study, to limit the total number of directional-seasonal bins and hence the number of parameters to be inferred in the analysis.Nevertheless, we are careful to allow for possible directional effects due to storms from the Atlantic and Norwegian Sea.There are clear directional and seasonal effects present.
Figure 9: Directional (left) and seasonal (right) variation of overturning moment (OTM, in mega Newton metres, MNm, top), storm peak significant wave height (H S in m, second row) and associated wind speed (WS in metres per second, ms −1 , bottom).Also shown are bin edges for three directional bins coupled with two seasonal bins (and hence a total of 6 = 4 × 2 directional-seasonal bins).The variation of OTM, H S and WS is approximately independent of covariates within bins, but different between bins.For the directional convention, see caption of Figure 2. Season (or seasonal degree) is defined by mapping the calendar year linearly into (0,360], with 0 indicating midnight on the 1st January. The resulting scatter plots of WS on OTM per covariate bin are shown in Figure 10 (and the corresponding plot for H S on OTM in Figure 28 of the user guide).Illustrative diagnostic plots for the estimation of marginal models for each of H S , OTM, and WS with direction and seasonal covariates are given in Figures 30-32 of the user guide.Marginal model parameters show clear directional and seasonal variation, and comparison of empirical and model-based tails suggest reasonable model fit.The fitted marginal models are then used to estimate the cumulative distribution functions for T -year maxima of interest; estimates for the 10-and 100-year maximum of WS are given in Figure 11.Unsurprisingly, the "omni" distribution (estimated over all directional and seasonal bins) is dominated by winter storms from the Atlantic sector, and the most probable 100-year maximum WS is approximately 29 ms −1 .
Following transformation to standard Laplace scale, a non-stationary conditional extremes model is estimated for H S and WS jointly, conditional on large OTM, given directional and seasonal covariates.Figure 12 gives parameter estimates for WS, showing directional and seasonal dependence of the conditional extremes slope parameter α; as would be expected perhaps, the dependence between WS and OTM is strongest for winter storms emanating from the North-West.The estimates for index β, common over covariate bins, is somewhat larger than found in Case Study 1.The empirical densities of nuisance parameters µ and σ are approximately Gaussian-shaped.Supporting diagnostic  plots for the conditional extremes model fit are given in Figures 33-38 of the user guide.In particular, Figure 37 indicates that the value of the α parameter is relatively stable for conditional extremes threshold levels τ ∈ [0.8, 0.9], and Figure 38 indicates that a minimum roughness penalty for predictive lack of fit is achieved at around 10. Figure 13 gives estimated conditional return value distributions for H S and WS given 10-year and 100-year maximum OTM for individual directional-seasonal covariate bins, and "omni" over all covariate bins.For H S , conditional return values are again largest for Atlantic winter storms.For WS however, we observe an interesting transition involving winter storms from directional sectors [275,315) to [230,275).The most probable conditional value of WS given a 100-year maximum OTM is approximately 26.5 ms −1 , lower than the marginal most probable maximum 100-year wind speed.
The resulting omni-covariate environmental design contours for H S and WS, given occurrences of the 10-year and 100-year maximum OTM, are shown in Figure 14.The general features of the contours are similar to those of Figure 7. Corresponding illustrative plots of contours per covariate bin for WS are given in Figure 40 of the user guide.

Discussion
This article introduces the covXtreme software for pragmatic multivariate extreme value analysis with covariate non-stationarity.The MATLAB software provides functionality to isolate temporal peaks from time-series, and to define an appropriate partition of the covariate domain.Using this partition, marginal generalised Pareto models are estimated for each variate independently assuming piecewise constant threshold and scale parameterisations within a penalised likelihood framework; optimal scale parameter roughness is estimated using cross-validation.Marginal models are then used to transform the sample to standard Laplace margins.A non-stationary conditional extremes model is then estimated with piecewise constant parameterisation for the slope ("α") parameter, again using penalised likelihood estimation with cross-validation to estimate optimal roughness.Simulations and importance sampling are then used to estimate the distribution of T -year maxima for each variate, and environmental design contours conditional on a single conditioning variate.Uncertainty due to marginal and dependence threshold selection is quantified by fitting multiple models with randomly chosen thresholds within user-specified plausible intervals of threshold non-exceedance probability.Uncertainties in parameter estimates and subsequent inferences from model fitting are estimated using bootstrap resampling.As noted in Section 1, the software has already been used in a number of studies, mostly but not exclusively metocean-related.
The covXtreme methodology makes a number of simplifying assumptions, motivated by the authors' experience of extreme value analysis applied to the ocean environment using a range of methodologies of difference complexities.For example, covXtreme relies on sensible user-specified partitioning of the covariate domain into bins within which it is reasonable to assume common marginal tails and a common dependence structure; this simplifies inference considerably compared with competitor approaches.Moreover, we believe that inferences using covXtreme with good partitioning are competitive with alternatives using more sophisticated tools.Specifically, covXtreme marginally is equivalent to a Voronoi set representation with pre-specified covariate partition, which was demonstrated by Zanini et al. (2020) to be competitive with P-spline and Bayesian adaptive regression spline covariate representations.covXtreme further assumes that the generalised Pareto shape parameter ξ in each marginal model is constant with covariate; because of the relative difficulty of estimating the shape parameter compared with the scale, this would appear reasonable in the absence of strong evidence to the contrary, especially for small samples of data.Likewise, the β, µ and σ parameters of the conditional extremes model are assumed stationary.This appears reasonable since β is an exponent, again difficult to estimate.Moreover, µ and σ are essentially nuisance parameters; any model misspecification caused by the assumption of stationarity will be accommodated to some extent by the adoption of residuals from model fitting for inferences under the model.It might be appropriate to relax some of the assumptions for specific applications, for example (a) when there is strong evidence that generalised Pareto ξ is unlikely to be constant (e.g.due to land shadow and fetch limitation effects on H S ), or (b) since parameter estimates for conditional extremes α and µ are highly correlated when β is close to unity.We also note that it is not clear whether non-stationary marginal extreme analysis is necessarily the best approach to estimate marginal return values from non-stationary data.Provided that sufficient sample is available, so that a sufficiently high threshold for peaks over threshold analysis can be set, often a stationary marginal analysis is at least competitive if not preferable; see Mackay and Jonathan (2020) for further discussion.
We anticipate that covXtreme might provide a pragmatic starting point to studies of spatial and temporal dependence of extremes, since the underlying methodology of both conditional spatial extremes (e.g.Shooter et al. 2019, Shooter et al. 2022, Wadsworth and Tawn 2022) and Markov extremal and similar time-series models (e.g.Winter and Tawn 2017, Tendijck et al. 2019, Tendijck et al. 2024) involves each of Stages 1-4 of the covXtreme approach.We also hope that covXtreme might be useful generally in estimating joint characteristics of extremes from multivariate time-series.

Figure 1 :
Figure 1: Schematic of covXtreme workflow.The five analysis stages are executed in sequence.Generally, at each stage, model tuning based on refinement of modelling hyper-parameters is necessary, guided by diagnostic information, before the next stage is attempted.Section numbers in the main article and user guide, providing relevant details for each stage of analysis, are given.

Figure 2 :
Figure 2: Directional variation of storm peak significant wave height (H S in metres, top) and associated spectral peak period (T P in seconds, bottom).Directions (in degrees) indicate the angle from which a storm emanates, defined clockwise from North.Also shown are directional bin edges (red) for 6 bins.The variation of H S and T P is approximately independent of direction within each bin, but different between bins.

Figure 3 :
Figure 3: Associated spectral peak period (s) on storm peak significant wave height (m) per directional bin.Panel titles indicate that the covariate is direction "D", and give the angular interval corresponding to the bin.It is apparent that the dependence between H S and T P varies between bins.

Figure 4 :
Figure 4: Marginal directional extreme value model for storm peak significant wave height (m).Variation of parameter estimates for GP scale ν, gamma shape ω and scale κ with direction.Mean estimates in solid black, and bootstrap 95% uncertainty bands in dashed black.Also shown are directional bin edges (red).Directional dependence in particularly clear for ν and κ.See Equations 2 and 4 for model forms and parameter definition.

Figure 5 :
Figure5: Cumulative distribution functions for the 10-year (left) and 100-year (right) maximum of storm peak significant wave height (m) per directional bin and over all bins ("omni", black).Horizontal dashed lines drawn at the exp(−1) quantile and median.

Figure 6 :
Figure 6: Parameter estimates from the conditional extreme value mode for associated peak period given storm peak significant wave height (m).Top left: directional variation of estimated linear slope parameter α summarised as mean and 95% bootstrap uncertainty band; top right: histogram of estimated exponent parameter β; bottom left: histogram of estimated residual mean µ; bottom right: histogram of estimated residual standard deviation σ.See Equations 12 and 13 for model form and parameter definition.

Figure 7 :
Figure 7: Omni-directional environmental contours of associated peak period (s) and storm peak significant wave height (m), for 10and 100-year maximum values of storm peak significant wave height shown as solid and dashed lines respectively, corresponding to the Exceedance (Exc, blue), Heffernan-Tawn density (HTDns, orange) and Huseby (Hus, yellow) contours.Lock points for the respective return periods are shown in green.

Figure 8 :
Figure 8: Environmental contours of associated peak period (s) and storm peak significant wave height (m), for individual directional bins.Contours shows for 10-and 100-year maximum values of storm peak significant wave height, shown as solid and dashed lines respectively, corresponding to the Exceedance (Exc, blue), Heffernan-Tawn density (HTDns, orange) and Huseby (Hus, yellow) contours.Lock points for the respective return periods are shown in green.

Figure 10 :
Figure10: Associated wind speed (ms −1 ) on storm peak overturning moment (MNm) per directional-seasonal bin.Panel titles indicate that the angular intervals of direction "D" and season "S" for the bin "D".It is apparent that the dependence between WS and OTM varies between bins, with obvious non-linearity in some bins.A similar plot is generated for associated significant wave height (m) on storm peak overturning moment.

Figure 11 :
Figure11: Cumulative distribution functions for the 10-year (left) and 100-year (right) maximum of associated wind speed (ms −1 ) per covariate bin and over all bins ("omni", black).Horizontal dashed lines drawn at the exp(−1) quantile and median.

Figure 12 :
Figure 12: Parameter estimates from the conditional extreme value mode for associated wind speed given storm peak overturning moment.Top left: directional and seasonal variation of α summarised as mean and 95% bootstrap uncertainty band; top right: histogram of β; bottom left: histogram of µ; bottom right: histogram of σ.Note the very high directional dependence for α at around 280 • .See Equations 12 and 13 for model form and parameter definition, and Figure 6 for comparison.

Figure 13 :
Figure 13: Conditional cumulative distribution functions of associated significant wave height (m, top) and associated wind speed (ms −1 bottom) per directional bin and over all bins ("omni", black), conditional on the 10-year (left) and 100-year (right) maximum of storm peak overturning moment (MNm).Horizontal dashed lines drawn at the exp(−1) quantile and median.

Figure 14 :
Figure 14: Omni-directional-seasonal environmental contours of associated significant wave height (m, left) and associated wind speed (ms −1 , right) and storm peak overturning moment (MNm).Contours shows for 10-and 100-year maximum values of storm peak significant wave height, shown as solid and dashed lines respectively, corresponding to the Exceedance (blue), Heffernan-Tawn density (HTDns, orange) and Huseby (yellow) contours.Corresponding contours per covariate bin are also generated.
Heffernan and Tawn (2004)ibution (with zero mean and unit variance) on W d ′ , appropriate when the dependence between Y d ′ and Y 1 is low.δd′ = 2 imposes a standard Gaussian distribution on W d ′ , the original assumption ofHeffernan and Tawn (2004), appropriate otherwise.For estimation purposes therefore, from the properties of the generalised Gaussian distribution, Y d ′ |(Y 1 = y) is assumed to follow a generalised Gaussian distribution with mean m bd ′