Temporal evolution of the extreme excursions of multivariate k th order Markov processes with application to oceanographic data

We develop two models for the temporal evolution of extreme events of multivariate k th order Markov processes. The foundation of our methodology lies in the conditional extremes model of Heﬀernan and Tawn (2004), and it naturally extends the work of Winter and Tawn (2016, 2017) and Tendijck et al. (2019) to include multivariate random variables. We use cross-validation-type techniques to develop a model order selection procedure, and we test our models on two-dimensional meteorological-oceanographic data with directional covariates for a location in the northern North Sea. We conclude that the newly-developed models perform better than the widely used historical matching methodology for these data.


Introduction
Farmers, stock brokers and sailors have one thing in common: they or their businesses are most heavily affected by extreme events like droughts and rainfall, stock market crashes, or extreme winds and waves, respectively.Understanding the statistical behaviour of such events as a whole is crucial for risk analyses.To make this more precise, if we let (X t ) t∈Z be a stationary d-dimensional random process of interest, then we seek to model excursions of the process in and out of a set E ⊂ R d in time, i.e., the behaviour of where E is associated with extreme events of the random variable X which is identically distributed to any X j , j ∈ Z.Moreover, we assume that the random process consists of multiple components that can be extreme.To solve this task, we assume that the multivariate random process is a realisation of a kth order Markov chain.We use extreme value theory, a subfield of statistics, to characterise excursions.There is considerable attention to this area in the literature, but most of extreme value theory for stationary Markov chains dates back over 20 years.Rootzén (1988) and Perfekt (1997) develop limiting results for univariate Markov chains and multivariate Markov chains, respectively.Smith (1992) calculates the extremal index (Leadbetter et al., 1983) for a univariate Markov chain and Smith et al. (1997) use parametric bivariate transition distributions to model the extremes of a univariate first order Markov process.Finally, Yun (2000) develops asymptotic theory for functionals of univariate kth order Markov extreme events.All of these authors derive results under the assumption of asymptotic dependence (Joe, 1997), i.e., for a stationary process (X t ) t∈Z satisfying suitable long-range mixing conditions, under the assumption that for any lag l = 1, 2, . . .lim u→x * P(X t+l > u|X t > u) > 0, where x * is the right upper end point of the distribution of X t .This early work doesn't consider what happens when asymptotic independence is present, i.e., when this limiting probability converges to 0 for some l.The first paper which considers such processes is Bortot and Tawn (1998) who assume a first order Markov model, with Ledford and Tawn (2003) considering a general framework for the modelling of asymptotic independent processes, and key recent probabilistic developments given by Papastathopoulos et al. (2017) and Papastathopoulos et al. (2023).Randell et al. (2015) speculate that a statistical model for the evolution of (multivariate) trajectories would be a valuable enhancement of description of ocean storm events.The first statistical work the current authors are aware of, that defines a model for the distribution of all observations during an excursion is Winter and Tawn (2016), who assume a flexible univariate first order Markov process exhibiting either asymptotic independence or asymptotic dependence across lags.Winter and Tawn (2017) incorporate higher order dependence model to give kth order Markov processes with k > 1.Finally, Tendijck et al. (2019) extend that model to a kth order univariate Markov process with a directional covariate.We remark that their work cannot be considered to model the extremes of bivariate Markov processes since the associated directional covariate does not take on extreme values.Feld et al. (2015) use a sophisticated covariate model for the most extreme observation (the most extreme value of the dominant variable) in an excursion, combined with a historical matching approach for the intra-excursion trajectory; in Section 3.4 we adopt a version of this methodology as a benchmark for our case study.Finally, we mention well-established literature on multivariate time series, e.g., Tiao and Tsay (1989), which is not directly applicable to modelling environmental extremes because such models are only designed to model typical behaviours.Financial timeseries models, e.g., Bauwens et al. (2006), are also not applicable because these are specifically tailored to model data exhibiting volatility, with tail switching during extreme events (Bortot and Coles, 2003).
In this work, we present a natural extension to Tendijck et al. (2019) by defining two multivariate kth order Markov models that exhibit both asymptotic (in)dependence across variables and/or at some lags.The work is motivated by our case study in which we model excursions of meteorological-oceanographic (met-ocean) data: significant wave height, wind speed, and their associated directions, for a location in the northern North Sea.
We use the following set up.Assume that at each time t ∈ Z, the distribution of the d-dimensional random variable X t is stationary through time; that is, X t has the same distribution as some X = (X 1 , . . ., X d ) with distribution function F X .For 1 ≤ j ≤ d, write F Xj as the jth marginal distribution of F X .The distribution functions F Xj are unknown and must be estimated.For extreme arguments of F Xj , we use univariate extreme value theory to motivate a class of parametric tail forms.More precisely, we assume that for each 1 ≤ j ≤ d, the excesses tail above some high level u j ∈ R of the marginal distribution F Xj are approximated with a generalised Pareto distribution (Davison and Smith, 1990).For non-extreme arguments x < u j of the function F Xj , an empirical model usually suffices.
In multivariate extreme value theory, it is common to consider the marginals and the dependence of random variables separately, such that the usually-dominant marginal effect does not influence the modelling of a possibly complex dependence structure.So given the marginal models as discussed above, we transform the random process (X t ) t∈Z onto standard Laplace margins (Y t ) t∈Z using the transformation: , where F −1 L is the inverse of the standard Laplace distribution function.Here the choice of Laplace margins is made to allow for the modelling of potential negative dependence at certain lags or across components (Keef et al., 2013).
For multivariate random processes, there are many ways of defining an extreme event.In our case study, we take the met-ocean variable significant wave height H S as the excursion-defining component.We follow Winter and Tawn (2017) and Tendijck et al. (2019) in adopting the conditional extremes model of Heffernan and Tawn (2004), see also Section 2.2, as the foundation of our approach.Without loss of generality, we first define the component X 1 of X as the defining variable for the extreme events.So, we set our excursion set for some high threshold u ∈ R + and rewrite our definition of an excursion as for a, b ∈ Z, indices for the start and the end time points of the excursion, respectively.In shorthand, the excursion is then Y a:b .We remark that in this definition, we accept that multiple excursions can occur close together in time, and thus these cannot be considered independent.The reason for this choice is that imposing a minimal separation of excursions would complicate the modelling significantly.We recognize that this is a feature of the current approach which can be improved upon in future work.
The remaining part of this paper is organised as follows.In Section 2, we present our strategy for modelling excursions by defining time intervals corresponding to so-called "pre-peak", "peak" and "postpeak" periods, and we present our kth order Markov models for each of these time periods.In Section 3, we apply the two Markov model forms we propose to met-ocean data for a location in the northern North Sea.We compare the model performance with a baseline historical matching approach by assessing their respective performance in estimating the tails of the distributions of complex structure variables (Coles and Tawn, 1994), corresponding to approximations of the response of hypothetical offshore or coastal facilities to extreme met-ocean environments.We find that in general the new models are preferred.2 The models

Modelling strategy
To model excursions as in definition (2), two types of approaches have been proposed in the literature of univariate extremes: a forward model (Rootzén, 1988) and a peak model (Smith et al., 1997).Both of these are two-step approaches by nature.The forward model first describes the distribution of a random exceedance Y t > u with a univariate extremes model and a conditional model for the distribution for any . ., j) where y t > u.Even though this approach does not directly model the univariate equivalent of excursions in formulation (2), estimates of some extremal properties of the process (Y t ) t≥1 , such as the extremal index (Leadbetter et al., 1983), can still be obtained by allowing the excursion threshold to be significantly lower than the cluster threshold used in extremal index estimators.Notably, Winter andTawn (2016, 2017) use the forward approach in their work.The peak model, on the other hand, does model excursions as defined here.This method relies on a univariate extremes model for the largest observation of an excursion, e.g., Eastoe and Tawn (2012), and a conditional model for observations before and after the excursion maximum.Winter and Tawn (2016) use this approach for their first order model but not for their kth order model (Winter and Tawn, 2017).They avoid this method explicitly because of difficulties that arise in preserving model characteristics in forward and backward simulations near the excursion maximum (i.e., the time point at which the defining variate X 1 achieves its maximum value during the excursion).Tendijck et al. (2019) use the peak method, but they do not address the issues associated with forward and backward simulation under the method.Because the excursion maximum is usually the most important observation of an excursion for risk assessments, we also use the peak method in the current work, but with consideration of backward and forward models.We separate the modelling of excursions into three stages: the modelling of the period of the peak, and the modelling of the pre-peak and post-peak periods; see Figure 1 in which the three time periods are illustrated for k = 3.Without loss of generality, let t = 0 be the time point at which the first component Y t,1 takes its maximum value within an excursion such that Y 0,1 > u for the threshold u.The period of the peak P k 0 of an excursion of a kth order model is then defined as the set of 2k − 1 observations: The pre-peak P pre and post-peak P post periods are defined as the sets of observations that include the excursion maximum and the observations before and after, respectively: P pre := {Y t : t ≤ t ≤ 0, with t = min{s < 0 : min i=s,...,0 {Y i,1 } > u}} and P post := {Y t : 0 ≤ t ≤ t , with t = max{s > 0 : min i=0,...,s {Y i,1 } > u}}, so each of them intersects with P k 0 .The length of P k 0 can be longer or shorter than the length of an excursion if the excursion ends within the period of the peak.We choose to define the period P k 0 in this manner so that the pre-peak and post-peak parts of the excursion are both initialized with k observations.
We then model an excursion as follows: (i) we model the excursion maximum Y 0,1 using a generalised Pareto distribution; (ii) we model the period of the peak P k 0 conditional on the storm maximum Y 0,1 using the model described in Section 2.2; (iii-a) if min j=1,...,k−1 Y j,1 < u (min j=1,...,k−1 Y −j,1 < u), then the period P post (P pre ) of the excursion has ended; (iii-b) if min j=1,...,k−1 Y j,1 ≥ u (min j=1,...,k−1 Y −j,1 ≥ u), then the remaining part of the excursion is modelled with our time-series models from Sections 2.3-2.4 until there exist a j 1 , j 2 > 0 such that Y j1,1 < u and Y −j2,1 < u; (iv) if max −j2≤i≤j1 Y i,1 > Y 0,1 , then the model for the excursion contradicts the definition of the period of the peak of an excursion, and so we reject such occurrences.
In the next sections, we discuss forward models that are applicable to model the post-peak period P post .We model the pre-peak period P pre using the forward models applied to (Y −t ) t∈Z (with potentially different parameters, although these would be the same if the process was time reversible).Importantly, we do not impose consistency in the forward and backward models to yield a kth order Markov chain, e.g., in the case of asymptotic dependent Markov chains the precise dependence conditions between the forward and backward hidden tail chains are given by Janßen and Segers (2014).We make this choice for two reasons: (i) for environmental applications, such as in this work, the pre-peak and post-peak period have different distributions, see for example the asymmetry in Figure 5, which is due to different physics in the growth and decay of a storm; (ii) the assumption of a kth order Markov process is an approximation for the process that generates our data.Thus, imposing forward and backward consistency for a kth order Markov chain is likely to yield worse results for our application.So, we consider the violating of this assumption as a benefit more than a limitation as it can yield more flexible descriptions of excursions.

The conditional extremes model
We introduce the conditional extreme value model of Heffernan and Tawn (2004), henceforth denoted the HT model, with notation specific to modelling the period of the peak P k 0 .The HT model is widely studied and applied to extrapolate tails of multivariate distributions, e.g., in oceanography (Ross et al., 2020), finance (Hilal et al., 2011), spatio-temporal extremes (Simpson and Wadsworth, 2021), and multivariate spatial extremes (Shooter et al., 2022).The HT model is a limit model and its form was originally motivated by deriving possible limiting forms for numerous theoretical examples. Let be a random matrix on R (2k−1)×d with standard Laplace margins (Keef et al., 2013), and define the irregular random matrix Y to be Y −(k−1):(k−1) without the (k, 1)th element Y 0,1 .That is, we define the irregular matrix x ∈ R (2k−1)d−1 as follows: , such that x does not contain the (k, 1)th element.Equivalently, we can write x = x −(k,1) for x ∈ R (2k−1)×d .Additionally, we assume that the joint density of Y −(k−1):(k−1) exists.
In practice, we exploit these results by assuming they hold exactly above some high finite threshold u > 0. So, we approximate the conditional distribution of Y|Y 0,1 = y for y > u, y ∈ R (2k−1)d−1 as and we assume independence of (Y − αY 0,1 )Y −β 0,1 and Y 0,1 .There is no finite-dimensional parametric form for H, so non-parametric methods are typically applied.However, we remark that there are applications of the conditional extreme value model where the copula H is assumed to be Gaussian (Towe et al., 2019) or a Bayesian semi-parametric model is used (Lugrin et al., 2016).For inference, see Section 2.5.

Multivariate Markov extremal model
For ease of presentation, we present the multivariate Markov extremal model (MMEM) of order k only for a two-dimensional time-series (Y t ) t∈Z such that Y t = (Y t,1 , Y t,2 ) in the notation of Section 1, i.e., Y t has standard Laplace margins.We only describe a forward model that is applicable to the post-peak period P post , since the backward model has a similar construction.As mentioned in Section 2.1, we apply a different forward MMEM model to (Y −t ) t∈Z to yield the backward model for the pre-peak period P pre .Concisely put, the MMEM exploits the HT model to estimate the distribution for Y t+k conditional on (Y t , . . ., Y t+k−1 ) when Y t,1 > u for a large threshold u > 0. As in Section 2.2, for each t ∈ Z, we define xt ∈ R k × R k+1 to be an irregular matrix with k + 1 rows and 2 columns without the element that is on the first row and first column: .
Then, we assume that for a large threshold u > 0, there exist parameters α0 Winter and Tawn (2017), for t ∈ Z, j ≥ 1 when Y t+j,1 > u, we then get For inference, we refer to Section 2.5.

Extremal vector autoregression
Here, we introduce extremal vector autoregression (EVAR) for extremes of the process (Y t ) t≥1 .This model combines the HT model with a vector autoregressive model for the joint evolution of the time-series at high levels.Here we focus on the post-peak period, but note that the pre-peak period is modelled analogously.We define an EVAR model of order k with parameters with Y t,1 = y for y > u, where u > 0 is a large threshold and ε t is a d-dimensional multivariate random variable that has non-degenerate margins and is independent of (Y t , . . ., Y t+k−1 ).Usually for a vector autoregressive model, parameter constraints would be imposed so that the resulting process is stationary.
In the current extreme value context, stationarity is not of concern to us, since we reject trajectories that exceed the excursion maximum, and stop the process once the first component dips below threshold u.We define EVAR 0 as a special case of EVAR corresponding to B = 0. EVAR 0 therefore has clear similarities with a regular vector autoregressive model (Tiao and Box, 1981), yet we emphasise that there is considerable difference between the two, since the parameters of EVAR 0 do not need to yield a stationary process, and the parameters of EVAR 0 are estimated using only extreme observations.To estimate the EVAR model, we adopt the same approach as that used to estimate the HT model, see Section 2.5.As explained in Appendix A, the resulting parameter estimators Φ(i) are highly correlated.Hence a reparameterisation is introduced to reduce this correlation, and improve inference efficiency and computation.
For practical applications, an advantage of EVAR over MMEM is that it provides a lower-dimensional residual distribution when k > 1 (with dimensions d and kd, respectively).As a consequence, the EVAR residual distribution is less affected by the curse of dimensionality.A drawback of EVAR is that it might be insufficiently flexible to describe complex dependence well.

Inference for conditional models
We discuss inference for each of the conditional extremes, MMEM and EVAR models with parameter vector θ.We discuss these together because they can be summarized in the same form.Specifically, let W = (W 1 , . . ., W d ) be a d-dimensional random variable and assume that for some high threshold u > 0, for some parametric functions multivariate random variable that is non-degenerate in each margin and independent of W 1 .As an example, for MMEM, g 1,j (x) = α j x for some α j and g 2,j (x) = x βj for some β j .
3 Case Study -Northern North Sea

Overview
We apply MMEM, EVAR and a historical matching procedure (introduced in Section 3.4, henceforth referred to as HM) to characterise excursions of significant wave height H S and wind speed W s with directional covariates for a location in the northern North Sea.Our goal is to estimate parsimonious predictive models for the joint evolution of H S and W s time-series conditional on H S being large.In Section 3.2, we describe the available met-ocean data.In Section 3.3, we outline a model for the evolution of storm direction that is needed for our time-series models.Section 3.4 then summarises the HM procedure, and in Section 3.5, we introduce structure variable responses that approximate fluid drag loading on a marine structure such as a wind turbine or coastal defence.Finally, in Section 3.6, we compare the predictive performance of MMEM and EVAR (over a set of model orders) with the HM method in estimating structure variables for withheld intervals of time-series.

Data
We have 53 years of hindcast data of time-series for four three-hourly met-ocean summary statistics at a location in the northern North Sea (Reistad et al., 2009): significant wave height (H S,i in metres), wind speed (W s,i in metres per second), wave direction (θ H i in degrees) and wind direction (θ W i in degrees) for each i ∈ T .To use MMEM and EVAR, we transform significant wave height and wind speed onto Laplace marginals: , e.g., using directional marginal extreme value models for the tails (Chavez-Demoulin and Davison, 2005), but ignoring seasonality.This part of the analysis has been reported on numerous occasions, see for example Randell et al. (2015).Because the marginal transformation includes direction as a covariate and because direction is not constant during an excursion, we also establish a model for the directional evolution of excursions in order to transform them between standard and original margins, see Section 3.3.
Let D L be the collection of the transformed data To define excursions in D L , we set the excursion threshold u equal to the 95% percentile of a standard Laplace distribution, i.e., u ≈ 2.3, yielding 1, 467 observations of extreme excursions E u .This choice of threshold is not unusual as similar conclusions are drawn for excursion thresholds that are slightly different from our original choice.
Figure 2 shows four intervals of the time-series chosen to contain the observations corresponding to the 100%, 95%, 90% and 85% sample percentiles of the set of excursion maximum significant wave heights, on original and standard Laplace margins, with directional covariates.Excursions are centred around extreme events.There is a large dependence of H S and W s on both original and standard margins.Moreover, variables associated to significant wave height, i.e., H S , H L S and θ H , are much smoother than their wind speed counterparts.Additionally, the directional covariates θ H and θ W centre around each other with no large deviations during extreme events.
In Figure 3, we visualize the (across variable joint) dependence of key variables H L S and W L s on Laplace scale at time lags up to lag 4 using a series of scatterplots where a unit of lag corresponds to three hours of observation time.The figure illustrates the complex dependence of the bivariate time-series of significant wave height and wind speed on Laplace margins.As expected, we observe (slow) convergence to an independent variable model as lag increases.Most notably, we observe a similar level of dependence of (H L S,t , W L s,t+4 ) and (W L s,t , W L s,t+4 ) which suggests counter-intuitively that H L S,t would be a better predictor for W L s,t+4 than W L s,t .In Figure 4, we plot (cross) correlation functions for these variables, and also for the change in directional covariates at various lags.These show that the dependence of (H L S,t , H L S,t+τ ) decays relatively slowly as τ grows to 90 hours, and that indeed the cross dependence between (H L S,t , W L s,t+τ ) is larger than the dependence of (W L s,t , W L s,t+τ ) for large τ .Finally, the correlation plot of the change in directional covariates ∆θ H S,i := (θ H S,i+1 − θ H S,i , mod 360) and ∆θ W s,i := (θ W s,i+1 − θ W s,i , mod 360) on the right shows that a first order model for these covariates is appropriate since the correlations nearly vanish at lag 2 (for wind and wave) or 6 hours (for all other combinations).

Directional model
We model wave direction θ H i in a similar fashion as Tendijck et al. (2019), summarised as follows.Let I ⊂ T be the set of indices of the original data that correspond to all observations of any excursion.Next, let {d(θ H i+1 , θ H i ) : i ∈ I} be the set of changes in wave directions, where d(θ, θ ) = (θ − θ + 180; mod 360) − 180 ∈ [−180, 180) denotes the circular difference of θ and θ in degrees.In our application, the set of changes in wave directions during excursions do not contain values close to −180 or 180.In particular, all of the observed changes centre around 0.
For i ∈ I, we transform observations ) on Gaussian margins, where F denotes the empirical distribution function of the set of changes in wave directions.Assume that {δ H i : i ∈ I} are realisations of the random variables {∆ H i : i ∈ I}.We estimate the following autoregressive model for ∆ H t of order p 1 = 1, 2, 3, . . .with parameters ϕ H j ∈ R for j = 1, . . ., p 1 as where ε t is a standard Gaussian random variable, and standard error ζ(h) is given by with λ j > 0 for j = 1, 2, 3, see Tendijck et al. (2019).In particular, the standard error ζ(h) decays as h grows due to the significantly larger amounts of energy needed to change the direction of more severe sea states.The parameters of this model are inferred with maximum likelihood, and in contrast to the inference discussed in Section 2.5, we do not reject the assumption that ε t is a standard Gaussian.In practice, we use p 1 = 1 in line with Tendijck et al. (2019).
Given model ( 8), we propose the following model for wind direction θ W t conditional on wave direction θ H t , where γ t is a zero-mean stationary AR(p 2 ) process.That is, there exist parameters ϕ W j ∈ R, 1 ≤ j ≤ p 2 , and a non-degenerate residual distribution r t independent of γ t−j for j ≥ 1, such that and such that the polynomial 1 − p2 j=1 ϕ W j z j has roots outside the unit circle.The model parameters and the distribution of r t are inferred as described in Section 2.5 conditional on the model order p 2 , which is selected by investigating the correlation function in Figure 4 and the partial autocorrelation function of γ t (not reported).In our application, we conclude that p 2 = 1 is sufficient.

Historical matching
An empirical method for simulating excursions is described in Feld et al. (2015) and termed historical matching (HM) in this work.They model trajectories of significant wave height, wave direction, season and wave period during extreme events.The key assumption they make is that storm trajectory (or excursion) profiles are not independent of storm maximum conditions.Specifically, the HM approach is a composition of four models: (i) a model for storm maximum wave direction; (ii) a model for storm maximum significant wave height conditional on storm maximum wave direction; (iii) a model that selects at random a historical storm trajectory with similar storm maximum characteristics to that simulated; (iv) a model that adjusts the historical storm trajectory by matching storm maximum characteristics of simulated and historical storms.
Specific details of the individual models are as follows, but this level of detail is not required for understanding the impact of the core methodology developments in Section 3.For model (i), we simply sample at random from the observed wave directions associated with storm maximum significant wave height (excursion maximum).In model (ii), storm maximum significant wave height are modelled using a generalised Pareto distribution conditional on the sampled storm maximum wave direction using a generalised additive model with the parameters as B-splines conditional on directional covariates (Chavez-Demoulin and Davison, 2005).In model (iii), we use a distance measure to calculate the dissimilarity between pairs of storm maximum significant wave heights and storm maximum wave directions for simulated and historical trajectories.Here, we use the heuristic recommended by Feld et al. (2015) ensuring that a difference of 5 degrees in storm maximum wave direction corresponds to the same dissimilarity as 0.5m of difference in storm maximum significant wave height; one of the closest 20 matching storms is then selected at random for associated with the simulated storm maximum.In model (iv), we match the variables of the chosen historical trajectory as follows: (a) the historical significant wave height series are multiplied by the ratio of the simulated maximum significant wave height to the maximum of the historical significant wave height; (b) the historical wave directions are shifted such that the storm maximum wave directions of simulated and historical trajectories coincide; (c) the associated historical wind directions are rotated in the exact same way as wave direction; (d) for the full set of historical storm maxima, storm maximum associated wind speed W M s (namely the value of wind speed at the time point corresponding to the storm maximum event) conditional on storm maximum significant wave height H M S is described using linear regression with parameters β 0 , β 1 ∈ R, σ > 0: with ε a standard normal random variable; (e) wind speed for the selected historical trajectory is scaled linearly such that it agrees with the storm maximum associated wind speed from (d).
Perhaps the main deficiencies of the HM approach are (i) it does not provide a means for modelling the extremal temporal dependence characteristics of excursions, and the extremal dependence between different components of the time-series for excursions to levels beyond those observed in the historical sample, and (ii) it does not provide a model framework for the assessment of fit or uncertainty propagation.

Response variable
To measure the practical impact of extreme met-ocean excursions, we define structure response variables for a simple hypothetical marine offshore facility.A structure response variable is a function of the met-ocean variables, key to assessing the integrity of the design of a physical structure of interest.Specifically, we consider a structure in the form of a unit cube standing above the water, supported by thin rigid legs, with vertical cube faces aligned with cardinal directions.Only wave and wind impact on the cube itself is of interest to us, and we neglect the effects of other oceanic phenomena such as swell, surge, tide, and potential climate non-stationarity.For simplicity, we also assume that when H S < h, for some known value h > 0, the wave impact on the structure is negligible, and structural response is dominated by wind.When H S ≥ h, we assume that wave impact increases cubically with H S and quadratically with W s (see Morison et al. 1950 and Ma and Swan 2020 for supporting literature).Hence, the impact of an extreme excursion on the structure is defined by the instantaneous response variable R where ] is the exposed cross-sectional area of the cube, see below, and the parameter c > 0 is specified such that both significant wave height and wind speed have an approximately equal contribution to the largest values of R.
Here both c and h are values that can be changed by altering structural features.The exposed cross-sectional area A(θ) ∈ [1, √ 2] of the cube is given by for a given wave direction θ H .The inline wind-speed I W is the component of the wind speed in the direction of the wave given by To simplify notation, we write R i (c, h) := R(H S,i , W s,i , θ H i , θ W i ; c, h) for i ∈ T .To define a structure response for a complete excursion E u , we write for some a < b such that for a threshold u > 0 (on Laplace margins) H L S,i > u for a ≤ i ≤ b and H L S,a−1 , H L S,b+1 ≤ u.Next, let i * := i * (E u ) be the time of the excursion maximum, i.e., H S,i * is the maximum of H S,i over E u .We define two natural structure response variables representing the maximum impact of an excursion max {a≤i≤b} R i (c, h), and the cumulative impact of an excursion {a≤i≤b} R i (c, h), respectively.For our application, we consider slight alterations That is, we consider responses that do not depend directly on the characteristics of the excursion near to the excursion maximum, to exaggerate the dependence of the structure variables on pre-peak and post-peak periods compared to the period of the peak, and hence the importance of estimating good models for the pre-peak and post-peak periods using MMEM or EVAR.Moreover, we define R max (c, h) and R sum (c, h) as the random structure responses related to a random excursion.

Model comparisons
Here, we use our time-series models to characterise extreme excursions for the met-ocean data D of Section 3.2 with structure responses R max and R sum .First, we investigate the model fits, then we describe our model comparison procedure, and finally we assess model performance using a visual diagnostic.We fit EVAR, EVAR 0 and MMEM with model orders k = 1, 2, . . ., 6 to data D L .The fitting of these 18 models is a two-stage procedure.In the first part, we fit (six) conditional extremes models for the period of the peak P k 0 for each k.In the second part, we fit 2 • 18 = 36 models to the pre-peak P pre and post-peak P post periods.In Table 1, we report parameter estimates of the period of the peak model, and in Tables 2-3, we report parameter estimates of MMEM on P post and P pre , respectively.Finally, we report parameter estimates of EVAR on P post and P pre in Tables 4-5, respectively.These indicate that all models agree on some level of asymptotic independence at each lag (coefficients of α0 are less than 1) with decreasing levels of dependence as lag increases, which can be seen by decreasing coefficients of ã0 for entries further down the table.We remark that for EVAR(2) on P pre , the coefficient of H S at time t + 1 (0.96) is larger than the coefficient of W s at time t + 1 (0.50) for estimating W s at time t + 2. This has the interpretation that significant wave height might be a better predictor for wind speed than wind speed itself, also suggested by Figure 4.
In the supplementary material, we produce analogous plots for each of the 18 models considered and HM.We observe that EVAR(4) characterizes the period of the peak, and also the pre-peak and post-peak periods of the excursion well.Moreover, EVAR(4) also reproduces the observed excursion survival probability.
Next, in Figure 6, we plot estimates of conditional probabilities χ H (u, l) , and χ W (u, l) := P(W L S,t+l > u | W L S,t > u) using EVAR, MMEM and HM with model orders 1 and 4, and we compare these with empirical estimates.1We make the following observations: HM is significantly worse at characterizing each of χ H , χ W and χ HW compared to EVAR and MMEM.Moreover, estimates obtained from EVAR of large enough order, e.g., k ≥ 4, agree well with empirical estimates.MMEM, on the other hand, yields estimators that are slightly positively biased.In particular, larger model orders yield considerable improvements.
In Figure 6, we discuss goodness-of-fit of each of the models.To compare MMEM and EVAR with each other and with HM, we take a similar approach to Gandy et al. (2022), who adjust standard cross-validation techniques to extreme value applications by taking a small training set and a larger test set.We select at random 25% of the observed excursions for our training sample; the remaining 75% forms our test sample.Below, we calculate performance statistics for the response variables by averaging over 50 such random partitions of the sample.
For training, we fit EVAR, EVAR 0 and MMEM with model orders k = 1, 2, . . ., 6 as explained in the second paragraph of this section.For each of the 18 models and HM, we simulate 20, 000 excursions, calculate structure response variables R max and R sum , and compare distributions of simulated structure response variables with those corresponding to the withheld test data.This is achieved by defining a dissimilarity distance function D that measures the level of difference in tails of distribution functions.We select 20 equidistant percentiles p 1 , . . ., p 20 ranging from 97% to 99.9% corresponding to moderately extreme to very extreme levels with respect to the (smaller) training sample but not too extreme for the (larger) withheld data.We define the distance D of distribution functions F M (of model M ) and F E (an empirical distribution function) as the mean absolute relative error over these percentiles, i.e.,
We remark that in the above definition, we never divide by zero because we only use D to measure the dissimilarity of distributions of positive random variables.
In Figure 7, we show the results for the 50 random partitions of the original sample by plotting the average distance D (with 80% confidence intervals) for each model together with HM for four different structure response variables corresponding to two choices of c and h for each of R max and R sum .Note that similar studies for other values of c and h for R max and R sum were examined, and general findings are consistent with those illustrated in Figure 7.For legibility, we omit confidence bands for EVAR 0 since the and EVAR(4) (red) excursions: median (solid), and the 10% and 90% quantiles (dashed).In the bottom panel, we plot survival probabilities for observed (black) and EVAR(4) (red) excursions relative to the time of the excursion maximum, see equation (10).Figure 6: Estimates of measures of extremal dependence across time lags 1 and 4, and variables given by χ H , χ HW and χ W (left, middle, and right respectively) for each of the models: EVAR (red), MMEM (blue), HM (green), data (grey).For EVAR and MMEM, we plot these estimates for different model orders k = 1 and k = 4 with line types: one (solid), four (dotted).Moreover, the grey region depicts the confidence bounds for empirical estimates of these extremal dependence measures from the data.
difference with EVAR is minimal.Model selection now involves choosing the model that yields the smallest average dissimilarity D whilst keeping the model order as low as possible.
We make a number of observations.For the R max response, EVAR and MMEM clearly outperform HM regardless of model order.However, for the R sum response, high order (e.g., k = 4, 5, 6) EVAR and MMEM are necessary to be competitive with HM.We observe also that performance of EVAR and MMEM does not significantly improve or worsen for k > 4.This finding is further supported with an unpublished study with Markov model orders of k ≤ 10.We note that llustrations of excursions in the supplementary material demonstrate that MMEM(1) does not explain the variability of the pre-peak and post-peak periods well.
By looking at the average relative errors in R max and R sum of our proposed selection of methods, we conclude that a third or fourth order MMEM and a fourth order EVAR are competitive models within their class.Since these models have similar performance, we prefer EVAR(4) because of its simpler twodimensional residual distribution.

Conclusions and discussion
In this paper, we provide models for extreme excursions of multivariate time-series.Excursions are characterized by a three-stage modelling procedure for the period of the peak, the pre-peak and the post-peak periods.We model the period of the peak using the conditional extremes framework (Heffernan and Tawn, 2004), and for the pre-peak and post-peak periods, we define two classes of time-series models: MMEM, motivated by the Markov extremal model of Winter and Tawn (2017); and EVAR, an extreme-value extension of a vector autoregressive model.We compare these excursion models with a baseline historical matching method, motivated by Feld et al. (2015).We find that the excursion models -for a reasonably informed choice of k, the order of the Markov process -are at least competitive with historical matching and often outperform it in the estimation of the tail of a range of notional structure response variables for a met-ocean application in the northern North Sea.
Statistical modelling of extreme excursions of multivariate time-series is difficult as it requires the estimation of complex model forms.MMEM requires the estimation of the conditional distribution of highdimensional residual random variables and EVAR is highly parameterized.Nevertheless, for realistically sized directional samples of significant wave height and wind speed time-series, we found that MMEM(3), MMEM(4) and EVAR(4) perform well.Even when the empirical historical matching procedure is competitive, adoption of an excursion model is advantageous because it allows for rigorous uncertainty quantification.We expect that our excursion models are applicable more generally, e.g., for the modelling of higher-dimensional met-ocean time-series and spatial fields.
We model wind speed and significant wave height marginally conditional on directional covariates.However, we did not investigate the explicit effect of the directional components on the dependence models.Since, we remove the marginal effect of direction before modelling the dependence, we do not expect this covariate to have a significant impact on the dependence.However, it would be very interesting to adapt our models to be able to investigate this further in future research.

A Reparameterization of EVAR
As opposed to inference for vector autoregressive models, we cannot estimate the EVAR parameters by least squares due to the presence of the Y B t,1 term.Instead, we apply the inference methodology discussed in Section 2.5.Not surprisingly, the parameter estimates Φ(i) for i = 1, . . ., k are highly intercorrelated because of the linear dependence between the components of Y t−1 , . . ., Y t−k .Reparameterization to reduce the correlation between parameter estimators is therefore attractive.
where αi,j is the maximum likelihood estimate for α i,j .Under this reparametrization, estimators of Φ(i) j,k are less correlated, which we demonstrated in unreported experiments comparing the dependence of the original and the reparameterized parameters using adaptive MCMC methodology (Roberts and Rosenthal, 2009).

Figure 1 :
Figure 1: Illustration of the periods of the pre-peak, peak and post-peak periods for two excursions from a Markov model with order k = 3.

Figure 2 :
Figure2: Intervals of oceanographic time-series: (top) key variables: significant wave height H S,i and wind speed W s,i on original margins; (middle) on Laplace margins; (bottom) covariates: wave direction θ H i and wind direction θ W i .The four columns correspond to time periods that contain the 100%, 95%, 90% and 85% empirical percentiles of H S,i , respectively.

Figure 3 :
Figure 3: Matrix plot of observed H L S,i and W L s,i at various time lags up to lag 4 (corresponding to 12 hours in real time) including cross dependece.

Figure 4 :
Figure 4: Estimated correlation and cross-correlation at various time lags of: (left) the key variables on Laplace margins: H L S,i and W L s,i ; (right) the covariates: change in wave direction ∆θ H i := (θ H i+1 − θ H i , mod 360), change in wind direction ∆θ W i := (θ W i+1 − θ W i , mod 360) and γ i , see definition (9).

Figure 5 :
Figure5: Excursions of H S and W s from EVAR(4) model (left; black), and data (middle; right) on original margins such that storm peak significant wave height is in [11.5, 12.5]; (right) summaries of the data (black) and EVAR(4) (red) excursions: median (solid), and the 10% and 90% quantiles (dashed).In the bottom panel, we plot survival probabilities for observed (black) and EVAR(4) (red) excursions relative to the time of the excursion maximum, see equation (10).

Figure 7 :
Figure7: Average mean relative errors of HM, EVAR, EVAR 0 and MMEM (dashed/dotted) and 80% confidence regions (shaded) for estimating the distribution of structure responses using 25% of data for training and 75% of data for testing.For details, see the text.

Table 1 :
Estimates of model parameters α and β for the period of the peak P k 0 with model order k = 4. Also shown in parentheses are 90% bootstrap confidence intervals.The structure of the irregular matrix estimates of α and β is explained in Section 2.2.

Table 2 :
Estimates of MMEM model parameters α0 and β0 with model order k = 4 for P post .Also shown in parentheses are 90% bootstrap confidence intervals.The structure of the irregular matrix estimates of α and β is explained in Section 2.3.

Table 3 :
Estimates of MMEM model parameters α0 and β0 with model order k = 4 for P pre .Also shown in parentheses are 90% bootstrap confidence intervals.The structure of the irregular matrix estimates of α and β is explained in Section 2.3.

Table 4 :
Estimates of EVAR model parameters (Section 2.4) with model order k = 1 (left), 2 (right) for P post .Also shown in parentheses are 90% bootstrap confidence intervals.

Table 5 :
Estimates of EVAR model parameters (Section 2.4) with model order k = 1 (left), 2 (right) for P pre .Also shown in parentheses are 90% bootstrap confidence intervals.