Statistics of extreme ocean environments: Non-stationary inference for directionality and other covariate effects
Matthew Jones David Randell Kevin Ewans Philip Jonathan






	




dur]Matthew Jones
man]David Randell
sar]Kevin Ewans
man]Philip Jonathancor1
[cor1]Corresponding author. Email: philip.jonathan@shell.com
[dur]Department of Mathematical Sciences, Durham University, Durham DH1 3LE, United Kingdom.
[man]Shell Projects & Technology, Manchester M22 0RR, United Kingdom.
[sar]Sarawak Shell Bhd., 50450 Kuala Lumpur, Malaysia.








Numerous approaches are proposed in the literature for non-stationarity marginal extreme value inference, including different model parameterisations with respect to covariate, and different inference schemes. The objective of this article is to compare some of these procedures critically. We generate sample realisations from generalised Pareto distributions, the parameters of which are smooth functions of a single smooth periodic covariate, specified to reflect the characteristics of actual samples from the tail of the distribution of significant wave height with direction, considered in the literature in the recent past. We estimate extreme values models (a) using Constant, Fourier, B-spline and Gaussian Process parameterisations for the functional forms of generalised Pareto shape and (adjusted) scale with respect to covariate and (b) maximum likelihood and Bayesian inference procedures. We evaluate the relative quality of inferences by estimating return value distributions for the response corresponding to a time period of 10 × the (assumed) period of the original sample, and compare estimated return values distributions with the truth using Kullback-Leibler, Cramer-von Mises and Kolmogorov-Smirnov statistics. We find that Spline and Gaussian Process parameterisations estimated by Markov chain Monte Carlo inference using the mMALA algorithm, perform equally well in terms of quality of inference and computational efficiency, and generally perform better than alternatives in those respects.



extreme covariate non-stationary smoothing non-parametric spline Gaussian process mMALA Kullback-Leibler












§ INTRODUCTION



Accurate estimates of the likely extreme environmental loading on an offshore facility are vital to enable a design that ensures the facility is both structurally reliable and economic. This involves estimating the extreme value behaviour of meteorological and oceanographic (metocean) variables that quantify the various environmental loading quantities, primarily winds, wave, and currents. Examples of such parameters are significant wave height, mean wind speed and mean current speed. These characterise the environment for a given short period of time within which the environment is assumed to be stationary.

The long-term variability of these parameters is however non-stationary, in particular with respect to time, space and direction. From a temporal point of view metocean parameters generally have a strong seasonal variation, with an annual periodicity, and longer term variations due to decadal or semi-decadal climate variations. At any given location, the variability of a particular parameter is also dependent on the direction; for example, wind forcing is typically stronger from some directions than others, and fetch and water depth effects can strongly influence the resulting magnitude of the waves. Clearly these effects will vary with location: a more exposed location will be associated with longer fetches, resulting in a more extreme wave climate. 

When estimating the long-term variability of parameters, such as significant wave height, the non-stationary effects associated with e.g. direction and season can be incorporated by treating direction and season as covariates. The common practice is to perform extreme value analysis of hindcast data sets, which include many years of metocean parameters, along with their associated covariates. Such data sets have all the information needed for input to covariate analysis.

From a design perspective, the metocean engineer is often required to specify return values for directional sectors such as octants centred on the cardinal and semi-cardinal directions. These directional return value estimates must be consistent with the estimated omnidirectional return value. In a similar manner, return values may be required corresponding to particular seasons or months of the year, consistent with an all-year return value. Clearly, therefore, efficient and reliable inference for non-stationary extremes is of considerable practical interest, requiring estimation of (a) the rate and (b) the size of rare events. This work addresses the latter of these objectives.

A non-stationary extreme value model is generally superior to the alternative “partitioning” method sometimes used within the ocean engineering community. In the partitioning method, the sample is partitioned into subsets corresponding to approximately constant values of covariate(s); independent extreme value analysis is then performed on each subset. For example, in the current work we might choose to partition the sample into directional octants, and then estimate (8 independent stationary) extreme value models for each of the octants. There are two main reasons for favouring a non-stationarity model over the partitioning method. Firstly, the partitioning approach incurs a loss in statistical efficiency of estimation, since parameter estimates for subsets with similar covariate values are estimated independently of one another, even though physical insight would require parameter estimates to be similar.  This problem worsens as the number of covariates and covariate subsets increases, and the sample size per subset decreases as a result. In the non-stationary model, we require that parameter estimates corresponding to similar values of covariates be similar, and optimise the degree of similarity during inference. For this reason, parameter uncertainty from the non-stationary model is generally smaller than from the partitioning approach. Secondly, the partitioning approach assumes that, within each subset, the sub-sample for extreme value modelling is homogeneous with respect to covariates. In general it is difficult to estimate what effect this assumption might have on parameter and return value estimates (especially when large intervals of values of covariates are combined into a subset). In the non-stationary model, we avoid the need to make this assumption.

Numerous articles have reported the essential features of extreme value analysis (e.g. ) and the importance of considering different aspects of covariate effects (e.g. ). <cit.> considers estimation of annual maxima from monthly data, when the distribution functions of monthly extremes are known.  <cit.> describes directional modelling of extreme wind speeds using a Fourier parameterisation. <cit.> models the long-term time series of significant wave height with non-linear threshold models. <cit.> reports that estimates for 100-year significant wave height from an extreme value model ignoring seasonality are considerably smaller than those obtained using a number of different seasonal extreme value models. <cit.> describes smooth extreme value models in finance and insurance. <cit.> provides a straight-forward description of a nonhomogeneous Poisson model in which occurrence rates and extreme value properties are modelled as functions of covariates. <cit.> uses a Bayesian hierarchical model to characterise extremes of lichen growth. <cit.> considers identification of changes in peaks over threshold using Bayesian inference. <cit.> uses a hierarchical model to identify location and seasonal effects in marginal densities of hourly maxima for wind speed. <cit.> considers seasonal non-stationarity in extremes of NOAA buoy records. <cit.> discusses estimation for return values for significant wave height in the South China Sea using a directional-seasonal extreme value model.

<cit.> explores the directional characteristics of hindcast storm peak significant wave height with direction for locations in the Gulf of Mexico, North-West Shelf of Australia, Northern North Sea, Southern North Sea, South Atlantic Ocean, Alaska, South China Sea and West Africa. Figure <ref> illustrates the essential features of samples such as these. The rate and magnitude of occurrences of storm events vary considerably between locations, and with direction at each location. There are directional sectors with effectively no occurrences, there is evidence of rapid changes in characteristics with direction and of local stationarity with direction. Any realistic model for such samples needs to be non-stationary with respect to direction.





The objective of this article is to evaluate critically different procedures for estimating non-stationary extreme value models. We quantify the extent to which extreme value analysis of samples of peaks over threshold exhibiting clear non-stationarity with respect to covariates, such as those in Figure <ref> or simulation case studies in Section <ref> below, is influenced by a particular choice of model parameterisation or inference method. The 6 simulation case studies introduced in Section <ref> are constructed to reflect the general features of the samples in Figure <ref>, with the advantage that the statistical characteristics of the case studies are known exactly, allowing objective evaluation and comparison of competing methods of model parameterisation and inference. Our aim is that the results of this study are generally informative about any application of non-stationary extreme value analysis. We generate sample realisations from generalised Pareto distributions, the parameters of which are smooth functions of a single smooth periodic covariate. Then we estimate extreme value models (a) using Constant, Fourier, B-spline and Gaussian Process parameterisations for the functional forms of generalised Pareto parameters with respect to covariate and (b) maximum likelihood and Bayesian inference procedures. We evaluate the relative quality of inferences by estimating return value distributions for the response corresponding to a time period of 10 × the (assumed) period of the original sample, and compare estimated return values distributions with the truth using Kullback-Leibler (e.g. ), Cramer-von Mises (e.g. ) and Kolmogorov-Smirnov statistics. We cannot hope to compare all possible parameterisations, but choose four parameterisations useful in our experience. Similarly, there are many competing approaches for maximum likelihood and Bayesian inference, and general interest in understanding their relative characteristics. For example, <cit.> compares maximum likelihood and Bayesian inference for the three-parameter Weibull distribution. In this work, we choose to compare frequentist penalised likelihood maximisation (see Section 2.3) with two Markov chain Monte Carlo (MCMC) methods of different complexities. Non-stationary model estimation is a growing field. There is a huge literature on still further possibilities for parametric (e.g. Chebyshev, Legendre and other polynomial forms) and non-parametric (e.g. Gauss-Markov random fields and radial basis functions) model parameterisations with respect to covariates. Moreover, in extreme value analysis, pre-processing of a response to near stationarity (e.g. using a Box-Cox transformation) is preferred (e.g. 

The outline of the paper is as follows. Section <ref> outlines the different model parameterisations and inference schemes under consideration. Section <ref> describes underlying model forms used to generate samples for inference, outlines the procedure for estimation of return value distributions and their comparison, and presents results of those comparisons. Section <ref> provides discussion and conclusions.




§ ESTIMATING NON-STATIONARY EXTREMES
 


Consider a random variable Y representing an environmental variable of interest such as significant wave height. The characteristics of Y are dependent on covariates such as (wave) direction, season, location and fetch. In this work we assume that a single periodic covariate θ (typically direction, or season) is sufficient to characterise the non-stationarity of Y. That is, we assume that Y|θ has a stationary distribution. For exceedances Y-μ(θ) of some high threshold μ(θ), extreme value theory suggest that the conditional distribution of Y-μ(θ) given that Y>μ(θ) can be approximated by the generalised Pareto distribution


    Pr( Y>y | Y>μ(θ) ) = 1/σ(θ)( 1 + ξ(θ)/σ(θ)( y - μ(θ) ) )^-1/ξ(θ)   .


For design purposes, return value distributions are typically estimated using peaks-over-threshold of significant wave height. The characteristics of these peaks, each corresponding to a different storm event, vary typically with respect to wave direction. Conditional on wave direction, the peaks are reasonably assumed to be independent of one another. To evaluate the relative performance of different approaches to non-stationary extreme value analysis, it is therefore natural to use random simulation from generalised Pareto models with known directional characteristics.




 §.§ Generalised Pareto model



We assume we observe a sample = {1,…,} of peaks over threshold drawn independently from a generalised Pareto distribution, the parameters of which are functions of corresponding observed covariate values = {1,…,}. The sample likelihood is a product of generalised Pareto (GP) likelihoods for each of the observations


    |,,,    = ∏_i=1^i|i,i,i
        = ∏_i=1^1/i(1+i(i-i)/i)^-1/i-1


where  and  are the shape and scale parameters as functions of covariate. We do not attempt to estimate the threshold function , assuming it is 0 for all covariate values. We also assume that the rate of occurrence  of exceedances of  varies with covariate, but that  is known (see Section <ref>). It is computationally advantageous (e.g. , ) to transform variables from (,) to the asymptotically independent pair (,), where = (1+). Inference therefore amounts to estimating the smooth functions  and , although we usually choose to illustrate the analysis in terms of  and . In practical application, estimation of μ(θ) is itself generally also problematic (e.g. ), particularly in the presence of non-stationarity (), but as necessary for inference as reliable estimation of ξ(θ) and σ(θ). We choose to focus on the latter in this work.




 §.§ Covariate parameterisations



To accommodate non-stationarity, we parameterise  and  as linear combinations of unknown parameters  and  respectively, where

    =  ,  and   =


and  and  are row vectors of basis functions evaluated at . We consider four different forms of basis function, corresponding to Constant (stationary), Fourier, Spline and Gaussian Process parameterisations for  and , as described below.

Physical considerations suggest that we should expect GP model parameters to vary smoothly with covariate. In general, the basis parameterisation introduced here permits estimation of functional forms for GP model parameters which are too variable (or too rough) with respect to covariate. We therefore need a mechanism to restrict the roughness of functional forms, such that their roughness is optimal given the evidence in the data. For each parameterisation, we therefore specify roughness matrices  (for = ,) to regulate the roughness of  with respect to  during inference. This ensures that the elements of  weight the individual basis functions in such a way that the resulting estimate is optimally smooth in some sense. The form of the roughness penalty term  is 1/2', for some roughness coefficient . The penalty is incorporated directly within a penalised likelihood for maximum likelihood inference, and within a prior distribution for  in Bayesian inference, as described in Section <ref>.




  §.§.§ Constant (stationary) parameterisation



In the Constant parameterisation, the values of  and  do not vary with respect to . We therefore adopt a scalar basis function which is constant across all values of covariate, so that =  = 1, and corresponding roughness matrices =  = 1. We do not expect the return value distributions estimated under this parameterisation to fare well in general in our comparison, since samples are generated from non-stationary distributions. Quality of fit is expected to be poor, at least in some intervals of covariate. However, many practitioners continue to use stationary extreme value models, perhaps with high thresholds to mitigate non-stationarity, in applications; inclusion of a stationary parameterisation provides a useful point of reference for comparison, therefore.




  §.§.§ Spline parameterisation



Under a Spline parameterisation, the vector of basis functions for each of  and  is made up of  local polynomial B-spline functions with compact support, joined at a series of knots evenly spaced in the covariate domain (e.g. )


    =  = [ 1 ⋯   ]   .


We specify roughness matrices =  = DD which penalise squared differences between adjacent elements of the coefficient vectors, where


    D = [ -1  1  0  ⋯  0;  0 -1  1     0;  ⋮        ⋱  ⋮;  0  0  0  ⋯  1 ]


is a (-1)× difference matrix. In this work we set p to 50.




  §.§.§ Fourier parameterisation



We use basis vectors composed of sine and cosine functions of  different periods


    =  = [      1  sin() sin(2)      ⋯  sin()  cos()      ⋯  cos() ]   .


The roughness matrix is computed by imposing a condition on the squared second derivative of the resulting parameter function. If we write


    () = ∑_k=1^( kcos(k) + ksin(k) )


where = or , and k and k are the parameters from  corresponding respectively to the sine and cosine functions of period k. The roughness criterion (from ) becomes


    = ∫_0^2π (”())^2 d = ∑_k=1^ k^4(k^2+k^2)


such that the penalty matrix can be written in matrix form as


    = diag(0,1,2^4,…,k^4,…,^4,1,2^4,…,k^4,…,^4)


for the =2+1 Fourier parameters (0, 1, …,, 1, …,). In this work, we set the value of n_p to 25 so that p=51.




  §.§.§ Gaussian Process parameterisation



For a set of  nodes {1, 2, …, p} on the covariate domain, we use a Gaussian Process parameterisation <cit.> , and relate each covariate input to a knot using the following basis vectors


    =  = [ 1 ⋯   ]


where the indicator functions j. are defined as


    j = 
    1    if  |-j|<|-k| ∀k≠j
    
    0    otherwise   .


Roughness matrices  (where = ,) are defined by the coefficient correlation matrix  via = ^-1, and the elements of  generated by a periodic squared exponential covariance function <cit.>


    jk = exp(-2/^2sin(j-k/2)^2)


where  are correlation lengths for each of the parameters, fixed to likely values by comparison with the covariate functions used to generate the data; a value of r_η for 0.6 was used throughout. ^-1 penalises on the angular difference between the j^th and k^th nodes, reducing in value from exp(0) at angular difference zero to exp(-2/^2) at angular difference 180^∘. Estimating the Gaussian Process parameterisation on the partitioned covariate domain, as opposed to fitting it to each of the data inputs, greatly reduces the number of parameters to estimate, and is physically reasonable provided that p is sufficiently large. In this work, we use p=50 equally-spaced nodes. Estimating a parameter for each data point would have made the computational burden for the Gaussian Process parameterisation significantly greater than that for any of the other parameterisations. For example, computations with roughness matrix  (of dimension p × p on the partitioned covariate domain) are more efficient that those using the N × N version of  (N ≈ 1000) defined per data input.




  §.§.§ Model complexity



We assume it is known from physical considerations that the extremal characteristics of the environmental variable of interest (e.g. significant wave height) vary smoothly with directional covariate. Specifically, we expect the form for the tail of the distribution to be homogeneous within a narrow directional sector of width ≈ 10^∘. This in turn suggests a suitable minimum complexity for the non-stationary model parameterisations considered in this work. For the Spline parameterisation, we set p=50 corresponding to 50 spline basis functions equally spaced on [0,360), with a distance between peaks of adjacent spline basis functions of 7.2^∘. We estimate the Gaussian Process on a regular partition of the covariate domain into p=50 bins; the distance between the centres of adjacent bins is again 7.2^∘. For the Fourier parameterisation, we use a Fourier order of n_p=25 (and p=51), such that half a wavelength of the highest frequency Fourier component again corresponds to 7.2^∘. In this way, we expect the directional resolution of these three model parameterisations to be comparable. Similarly, the number p of basis coefficients to be estimated in each of the three parameterisations is comparable. In this way, we hope to focus fairly in the analysis below on the different inferential challenges presented by Spline, Fourier and Gaussian Process parameterisations and different estimation schemes for problems of comparable complexity. The stationary Constant parameterisation is clearly less complex (with p=1), but a useful baseline for comparison: we expect inferences from the Constant parameterisation to be relatively poor, and therefore make obvious the need to consider non-stationarity.




 §.§ Inference procedures



We consider two methods for estimating parameters and return value distributions for the models and parameterisations described above, namely (a) maximum penalised likelihood estimation with bootstrapping to quantify uncertainties, and (b) (two forms of) Bayesian inference using Markov Chain Monte Carlo (MCMC). These are discussed below.




  §.§.§ Maximum likelihood estimation



We use an iterative back-fitting optimisation (see Appendix) to minimise the penalised negative log likelihood -|, ;λ_,λ_ with respect to  and  for given roughness coefficients λ_ and λ_, where


    -|, ;λ_,λ_   =    -|,  +  + 
       =    -|,  + 1/2'  + 1/2'    .


Here, -|, is the negative log sample GP likelihood from Section <ref> expressed as a function of  and , and  and  are additive roughness penalties. The values of λ_ and λ_ are selected using cross-validation to maximise the predictive performance of the estimated model, and bootstrap resampling is used to quantify the uncertainty of parameter estimates. The original sample is resampled with replacement a large number of times, and inference repeated for each resample. We use the empirical distributions of parameter estimates and return values over resamples as approximate uncertainty distributions.

We refer readers interested in further information on penalised likelihood methods to the work of <cit.>, <cit.>, <cit.>, <cit.>, outlined in applications to metocean by <cit.>.



  §.§.§ Bayesian inference



From a Bayesian perspective, all of , ,  and  are treated as parameters to be estimated. Their joint posterior distribution given sample responses  and covariates  can be written


    ,,,|,   ∝|,,,|||,|,


where |,,, is the sample GP likelihood from Section <ref> and prior distributions |, |, |, and |, are specified as follows. Parameter smoothness of GP shape and (modified) scale functions is encoded by adopting Gaussian priors for their vectors   of basis coefficients (for = ,), expressed in terms of parameter roughness 


    |∝^1/2exp(-/2)   .


The roughness coefficient  can be seen, from a Bayesian perspective, as a parameter precision for . It is assigned a Gamma prior
distribution, which is conjugate with the prior Gaussian distribution for . The values of hyper-parameters are set such that Gamma priors are relatively uninformative; a_η and b_η takes values of 10^-3 throughout this study for all parameterisations. The Bayesian inference can be illustrated by the directed acyclic graph show in Figure <ref>.





Estimates for , ,  and  are obtained by sampling the posterior distribution above using MCMC. We choose to adopt a Metropolis-within-Gibbs framework (e.g. ), where each of the four parameters is sampled in turn conditionally on the values of others. The full conditional distributions | and | of precision parameters are Gamma by conjugacy, and are sampled exactly in a Gibbs step. Full conditional distributions for coefficients  and  are not available in closed form; a more general Metropolis-Hastings (MH) scheme must therefore be used. 

There are a number of potential alternative strategies regarding the MH step for  (= ,). We choose to examine two possibilities: (a) a straightforward MH sampling of correlated Gaussian proposals for , and (b) the mMALA algorithm of <cit.>, exploiting first- and second-derivative information from the log posterior to propose candidate values for the full vector of coefficients in high-probability regions. Implementations are described in the Appendix. Henceforth we refer to these two schemes as MH and mMALA respectively for brevity. The MH approach is simple to implement, but is likely to generate MCMC chains which mix relatively poorly. The mMALA scheme is expected to explore the posterior with considerably higher efficiency; however, its implementation requires knowledge of likelihood derivatives.




  §.§.§ Comparing uncertainties



Parameter uncertainty is estimated by bootstrap resampling for ML inference and by sampling from the posterior distribution of parameters for Bayesian inference. These two approaches seek to estimate parameter uncertainty, but in different ways. Since Bayesian priors are chosen to be relatively uninformative, we expect - at least naively - that inferences concerning parameter estimates and return value distributions from the two approaches will be similar, but not the same. In general, the relationship between bootstrap uncertainty estimates and those from Bayesian inference is complex, and an open topic in the literature (e.g. ). A thorough theoretical analysis of this relationship in the current application is beyond the scope of the current work.




§ EVALUATION OF METHODS
 


This section describes evaluation of relative performance of different model parameterisations and inference schemes introduced in Sections <ref> and <ref>. We assess performance in terms of quality of estimation of distributions of return values corresponding to long return periods, estimated under models for large numbers of replicate samples of data from pre-specified underlying models.

We simulate 100 sample realisations, each of size 1000 from three different underlying models, described below and referred to henceforth as Cases 1, 2 and 3, and further simulate 100 sample realisations of size 5000 from the same triplet of underlying models, referring to these as Cases 4, 5 and 6 respectively. We next estimate extreme value models for all sample realisations, model parameterisations and inference schemes. We assume that any sample realisation (for any Case) corresponds to a period 𝒯 years of observation.  We then simulate 1000 replicates of return period realisations, each replicate consisting of observations of directional extreme values corresponding to a return period of 10×𝒯, and estimate the distribution of the maximum observed (the 10𝒯-year maximum return value) for all model parameterisations and inference methods, by accumulation from the 1000 replicates. Return value simulations under models estimated using Bayesian inference proceed by sampling a different vector of model parameter estimates at random from the estimated joint posterior distribution of parameters and simulating the appropriate return period of events from the corresponding distribution, for each of the 1000 replicates. For ML inference, for each replicate, we sample the vector of model parameter estimates at random from a set of 100 parameter estimate vectors of generated by the bootstrap analysis. Return value distributions are estimated omnidirectionally (that is, including all directions) and for 8 directional octants centred on the cardinal and semi-cardinal directions (by considering only those observations from the return period realisation with the appropriate directional characteristics). For each sample realisation from Cases 4, 5 and 6, we estimate return value distributions for all parameterisations but for only mMALA inference, since as will be discussed in Section <ref> below, the computational effort associated with any of mMALA, MH and MLE (with bootstrap resampling) for these Cases is large.

We quantify the quality of return value inference by comparing the empirical cumulative distribution function generated under the fitted model for each sample realisation with that from simulation under the known underlying Case. We quantify the discrepancy between empirical distribution functions by estimating Kullback-Leibler, Cramer-von Mises and Kolmogorov-Smirnov statistics. We visualise relative performance by plotting the empirical cumulative distribution function of the test statistic over the 100 sample realisations, for each combination of Case, model parameterisation and inference method. We also compare performance in terms of prediction of the 37.5^th percentile of the 10 𝒯-year return value distribution, since this is often used in metocean and coastal design applications; it corresponds approximately to the location of the mode of a Weibull distribution with shape parameter ≈ 2. However, we are not only interested in quality of inference, but also in computational efficiency. This is evaluated and illustrated in Section <ref>.

To complete the full analysis described here, comprised of 100 random samples of each of 6 Cases, required running 3 dedicated workstations (exploiting each of 48 cores and 196GB RAM per workstation) for approximately 10 weeks. All assessments of return value distributions, in terms of e.g. Kullback-Leibler divergence or a central percentile, are therefore based on 100 independent estimates per Case. We note that the precision with which the median and quartiles of the distribution of Kullback-Leibler divergence, or the distribution of a central percentile (or is bias), is therefore considerably higher than that with which extreme quantiles are estimated.

In practical application to metocean design using a sample of measured or hindcast significant wave height data, for example, it is critical to demonstrate that an iterative simulation algorithm such as MH or mMALA has produced a chain which has itself converged to the stationary distribution, and whether that distribution has been adequately explored. This can be achieved, for example, by comparing inferences from multiple independent chains (e.g. ). This diagnosis is critical to such applications, and is the main evidence that MCMC inference is valid. In the current work, visual inspection of trace plots was used to confirm that an adequate period of burn-in had been specified for all combinations of parameterisation and inference scheme. However since the underlying true models are known, we choose to use comparison with the truth as the basis to assess the relative performance of different model parameterisations and inference schemes. The finding of this work is that mMALA as implemented here behaves more reasonably as a default approach to inference (requiring less user intervention and fine tuning) than MH. There is no doubt that chain convergence diagnostics would have indicated this also, and might also have prompted refinement of the MH scheme in particular to improve its performance.




 §.§ Case studies considered
 


First, we describe model Cases 1-6 used to generate sample realisations for extreme value modelling. For each Case, potentially all of Poisson rate  of threshold exceedance, GP shape  and scale  of exceedance size vary as a function of covariate . The extreme value threshold  is fixed at zero throughout.



  
Case 1: For extreme value threshold = 0, we simulate 1000 observations with a uniform Poisson rate = 1000/360 per degree covariate, and a low order Fourier parameterisation of GP shape = sin() + cos(2) + 2 and scale = - 0.2 + (sin(-30))/10.



  
Case 2: For extreme value threshold = 0 and the same Fourier parameterisation of GP shape and scale as in Case 1, a non-uniform Poisson rate = max(sin()+1.1,0) × 1000/c_ρ, where c_ρ=∫_0^360max(sin()+1.1,0) d is used to simulate 1000 observations.



  
Case 3: For extreme value threshold = 0, the forms of each of ,  and  are defined by mixtures of between one and five Gaussian densities, as illustrated in Figure <ref>. Sample size is 1000.



  
Cases 4, 5 and 6: These cases are identical to Cases 1, 2 and 3 respectively, except that Poisson rate ρ is increased by a factor of five. Sample size is therefore 5000.

Figure <ref> illustrates typical sample realisations of Cases 1, 2 and 3. Parameter variation of GP shape ξ and scale σ with direction θ are identical in Cases 1 and 2. Poisson rate ρ is constant in Case 1 only. In Cases 2 and 3, ρ is very small at θ≈ 270^∘ leading to a sparsity of corresponding observations. ξ is largest (but negative) at θ≈ 120^∘ for Cases 1 and 2, leading to larger observations here. For Case 3, ξ is largest (and positive) at θ≈ 30^∘ leading to the heaviest tail in any of the Cases considered. Figure <ref> shows parameter estimates for ξ and σ, corresponding to the sample realisation of Case 2 shown in Figure <ref>, for different model parameterisations using mMALA inference. Visual inspection suggests that estimates of similar quality are obtained using all of Spline, Fourier and Gaussian Process parameterisations, but that the Constant parameterisation is poor. It is also apparent that identification of ξ is more difficult than σ. Corresponding plots (not shown) for  maximum likelihood and Metropolis-Hastings inference show broadly similar characteristics, as do plots for other realisations of the same Case, and realisations of different Cases. (Posterior) cumulative distribution functions of return values based on models for the sample realisations of Case 2 illustrated in Figures <ref> and <ref>, corresponding to a return period of ten times the period of the original sample, are shown in Figure <ref> for different model parameterisations and mMALA inference.

It can be seen omnidirectionally that the Constant parameterisation provides best agreement with the known return value distribution, despite the fact that parameter estimates in Figure <ref> do not reflect the directional non-stationarity present. In <cit.>, it is demonstrated that a stationary extreme value model may produce good estimates for omni-directional return value distributions in at least two situations. Firstly, when the extreme value threshold is set sufficiently high that only observations from the most extreme interval of the covariate domain are modelled, therefore threshold exceedances and estimated models will be effectively stationary. Secondly, it may be that different bias effects introduced by the stationarity assumption compensate for each other in estimation of return values at a certain return levels. For example, with reference to the sample illustrated in Figure <ref>, the constant model underestimates the value of GP shape in general, but does a relatively good job of identifying the maximum values of shape and scale in directions around 120^∘. The constant model does a very poor job for directions around 270^∘. The uncertainty in parameter estimates from the constant model is also lower since the number of parameters to estimate is lower for that model. Reference to Figure <ref> for Case 2 shows also that the rate of occurrence of threshold exceedances is larger for directions where the constant model performs well. These characteristics are reflected in the corresponding return value estimates in Figure <ref>. For directional octants, the constant model does relatively well for directions around 120^∘, but very poorly around 270^∘. However, the omnidirectional estimate is dominated by directions around 120^∘. The relatively good performance of the directional model around 120^∘, together with the large rate of occurrence of events there and the relatively small parameter uncertainties for the constant model, result in its providing good estimation (in this case) of the omnidirectional return value. This occurs despite the fact that the directional bias of the constant model is greater than that of any of the non-stationary models. In general however, it is not possible to know a priori how a constant model will perform in estimating the omnidirectional value. It is clear however that its directional bias will be larger than that of an appropriate directional model. Reference to Figure <ref> below, assuming that mMALA provides more satisfactory inference as implemented here, suggests that values of KL divergence are somewhat larger and more variable for the constant model compared to the Spline or Gaussian process parameterisations. Reference to Figure <ref> shows that the Constant model does very poorly for certain directional sectors. Omnidirectionally, and for 8 directional octant sectors, non-stationary model parameterisations perform similarly. However, it is clear that the Constant parameterisation does particularly poorly for the western and north-western sectors, for which the rate of occurrence of events is relatively low, and both ξ and σ are near their minimum values. Figure <ref> illustrates uncertainty (over all 100 sample realisations) in the cumulative distribution function of return value for Case 2, corresponding to a return period of ten times the period of the original sample, using the Spline parameterisation and mMALA inference. The median estimate for the return value distribution (over all 100 sample realisations of Case 2) is shown in solid grey, with corresponding point-wise 95% uncertainty band in dashed grey. The true return value distribution is given in solid black. There is good agreement in all sectors. We explore differences in inferences for return value distributions more fully in Section <ref>.




















 §.§ Assessing quality of inference
 


The criteria used to compare distributions of return values are now described. Since, for comparison only, we only have access to samples from distributions, where necessary we project empirical distributions onto a linear grid using linear interpolation, and evaluate grid-based approximations to facilitate comparison. Then we compare empirical return value distributions using each of the following three statistics. The Kolmogorov-Smirnov criterion compares two distributions in terms of the maximum vertical distance between cumulative distribution functions, as 0,1 = sup_x|1x-0x| . The Cramer-von Mises criterion evaluates the average squared difference of one distribution from a second, reference distribution, using 0,1 = ∫_-∞^∞(1x-0x)^20x dx. The Kullback-Leibler divergence compares distributions using the average ratio of logarithms of density functions 0,1 = ∫_-∞^∞log(0x/1x)0x dx; in this work, we use the approximation of <cit.>. The general characteristics of differences in return value inference due to model parameterisation and inference method were found to be similar for each of the three statistics. Only comparisons using Kullback-Leibler (KL) divergence are therefore reported here. We note that perfect agreement between 1x and 0x yields a minimum KL divergence of zero.

For illustration, Figure <ref> shows empirical cumulative distribution functions for the KL divergence between return value distributions (corresponding to a return period of ten times that the original sample) estimated under the true model and those estimated under models  of sample realisations with different parameterisations and mMALA inference for Case 2. The distributions of KL divergence from all non-stationary model parameterisations appear to be very similar, as might be expected from consideration of figures similar to Figure <ref>. However, the Constant model yields the best performance omnidirectionally in this case (since the corresponding distribution of KL divergence is shifted towards zero). In stark contrast, the Constant model does particularly badly in the eastern, south-western, western and north-western sectors. Figure <ref> gives empirical cumulative distribution functions of the KL divergence between return value distributions (corresponding to a return period of ten times that of the original sample) estimated under the true return value distribution and those estimated under models of sample realisations with Spline parameterisations and different inference procedures for the same Case. There appears to be little to choose between mMALA and MLE inference methods for this Case, with MH somewhat poorer. 

Figure <ref> summarises the characteristics of distributions for KL divergence corresponding to the omnidirectional return value distribution for all Cases, model parameterisations and inference methods considered in this work. In general, we note that all non-stationary parameterisations perform well with mMALA inference. With MH inference, performance is generally poorer, especially for Fourier parameterisation. MLE does better than MH. We note that the Constant parameterisation generally performs well for the omnidirectional return value, but there is some erratic behaviour, notably for Case 4. Figure <ref> is the corresponding plot for the (generally sparsely populated) western directional sector. The Spline and Fourier parameterisations with mMALA inference perform best. We note that the Fourier parameterisation does less well using MH and MLE, and that the Constant parameterisation behaves very erratically. We also note that, somewhat surprisingly, the Gaussian Process model performs considerably less well than the Spline and Fourier parameterisations. We surmise that this is due to mean-reversion in the absence of observations, compared with the other parameterisations which prefer interpolation to reduce parameter roughness. We note that the Constant parameterisation also performs badly for Case 4.

















From a practitioner's perspective, it is also interesting to quantify the performance of different model parameterisations and inference methods in estimating some central value (e.g. the mean, median, of 37.5^th percentile) of the distribution of the return value. We choose the 37.5^th percentile, since this percentile is near the mode of the distribution, and commonly used in the met-ocean community.


Figure <ref> illustrates this comparison in terms of a box-whisker plot. In each panel, the true value, estimated by simulation under the true model, is shown as a black disc. The distribution of estimates from 100 different sample realisations of each Case is summarised by the median (white disc with black central dot), the interquartile range (grey rectangular box) and the (2.5%, 97.5%) interval (grey line).

The performance of non-stationary parameterisations with mMALA is good, with the possible exception of Case 3 (and Case 6); yet MH inference tends to produce greater bias and variability in estimates of the 37.5^th percentile. We may surmise that this difference may be due to the fact that mMALA exploits knowledge of likelihood gradient and curvature. The Constant model again performs more erratically. Corresponding box-whisker plots for the eastern directional octant (not shown) show similar characteristics: for a given inference scheme, all non-stationary models yield similar performance, but the Constant model overestimates throughout. It is interesting that MLE shows some bias in return value estimation for the omnidirectional 37.5^th percentile, but this is not the case in general. True values of the 37.5^th percentile for the western sector (see Figure <ref>) are considerably lower than for the eastern sector, and lower again than the omnidirectional values. In this sector, the rate of occurrences of events is generally lower in all cases. Nevertheless, Figure <ref> has many similar features to Figure <ref>. However, we note that MH struggles in combination with the Fourier parameterisation, probably since the latter has the whole of the covariate domain as its support; intelligent proposals (like those used here in MLE and mMALA) are necessary. Estimates using the Constant model are erratic. Overall, we note that MLE and mMALA inference for all of Spline, Fourier and Gaussian Process parameterisations perform relatively well, and equally well.












 §.§ Assessing efficiency of inference
 


The effective sample size (m^*, e.g. ) gives an estimate of the equivalent number of independent iterations that a Markov chain Monte Carlo represents, and is defined by m^* = m/(1+2∑_k=1^∞ c_k), where c_k is the autocorrelation of the MCMC chain at lag k, and m is the actual chain length. The effective sample size per hour is defined by m^*/T, where T is the elapsed computational time (in hours) for m steps of the chain. For maximum likelihood inference with bootstrap uncertainty estimation, since bootstrap resamples are independent of one another, we estimate effective sample size per hour as m_BS/T where m_BS is the number of bootstrap resamples used and T is now the total elapsed computational time (in hours) to execute analysis of the m_BS bootstrap resamples. Comparison of effective sample sizes per hour for different cases, parameterisations and inference methods gives some indication of relative computational efficiency, although objective comparison is difficult. In particular we note that software implementations in MATLAB exploiting common computational structures between different approaches have been used; these are almost certainly to the detriment of computational efficiency for some of the approaches, particulary the Constant parameterisation. For this reason, we do not report effective sample size per hour for the Constant parameterisation. Computational run-times are also of course critically dependent on software and hardware resources used. We note that the focus of this work is primarily quality of inference, rather than its computational efficiency. Specific modelling choices, such as the set-up of the cross-validation strategy adopted for MLE and choice of burn-in length and proposal step-size for MH and mMALA within reasonable bounds may not influence inferences greatly, but will obviously however affect run times. Similarly the Spline, Fourier and Gaussian Process parameterisations used were chosen to be of similar complexity, but small differences may again influence relative computational efficiency of inference. With these caveats in mind, Figure <ref> illustrates the distribution of estimated effective sample size (ESS), and effective sample size per hour (ESS/hr) for different cases, model parameterisations and inference procedures.

The left hand side of Figure 12 illustrates log_10(ESS) for all combinations of parameterisations and inference schemes. For MLE inference, we choose to report the number of bootstrap resamples used, as described in Section 3. The effective sample sizes for mMALA inference are considerably larger than for MH for B-spline and Fourier parameterisations. For the Gaussian process parameterisation, mMALA still provides a larger ESS, but the difference between mMALA and MH is smaller. For mMALA inference, ESS is largest for the B-spline parameterisation; the Gaussian process parameterisation provides the smallest ESS on average. For MH inference, ESS for B-splines and Fourier parameterisations is near 10, suggesting that the posterior density has not been sufficiently explored due to poor MCMC mixing using the MH algorithm as implemented. The value of ESS is approximately constant across the different Cases examined for all combinations of model parameterisation and inference method. The right hand side of Figure 12 indicates that for B-spline parameterisation, ESS/hr is also higher for mMALA than for MH. For Fourier and Gaussian process parameterisations, ESS/hr is comparable for mMALA and MH. The right hand side of Figure <ref> shows that, for mMALA inference, the ESS/hr is considerably lower for Cases 4, 5 and 6, indicating that inference using large sample sizes is slower. For this reason, in this work, we do not provide results for Cases 4, 5 and 6 using MLE and MH. Overall, comparing non-stationary parameterisations, ESS/hr is larger for Splines and Gaussian Processes than for Fourier. There is little difference in ESS/hr for different model parameterisations using MLE.








§ DISCUSSION



Adequate allowance for non-stationarity is essential for realistic environmental extreme value inference. Ignoring the effects of covariates leads to unrealistic inference in general. The applied statistics and environmental literature provides various competing approaches to modelling non-stationarity. Adoption of non- or semi-parametric functional forms for generalised Pareto shape and scale in peaks over threshold modelling is far preferable in general in real world applications, than the assumption of a less flexible parameterisation. We find that B-spline and Gaussian Process parameterisations estimated by Bayesian inference (using mMALA) perform well in terms of quality and computational efficiency, and generally outperform alternatives in this work.

The Gaussian Process parameterisation is computationally unwieldy for larger problems, unless covariate gridding on the covariate domain is performed. The Fourier parameterisation, utilising bases with global support (compared to Spline and Gaussian Process basis functions whose support is local on the covariate domain), is generally somewhat more difficult to estimate well in practice, showing greater instability to choices such as starting solution for maximum likelihood estimation. The Constant parameterisation performs surprisingly well in estimating the omnidirectional return values distribution in some cases, but is generally very poor in estimating directional variation. 

Various choices of methods of inference are also available. Competing approaches include maximum (penalised) likelihood optimisation and Bayesian inference using Markov chain Monte Carlo sampling. It appears however that the major difference, in terms of practical value of inference, is not between frequentist and Bayesian paradigms but rather the advantage gained by exploiting knowledge of likelihood gradient and curvature. In addition, it appears that inference schemes which sample from a negative log likelihood surface randomly, rather that seeking its minimum deterministically, are more stable, and therefore more routinely implementable and useable. Moreover, Bayesian inference gives a more intuitive framework for statistical learning and communication of uncertainty, particularly to a non-specialist audience.

Here, we have focussed on the estimation of non-stationary shape and scale parameters for the conditional distribution of independent peaks over threshold. In practical application, careful estimation of non-stationary extreme value threshold is at least as important for reliable inference. We emphasise that for practical application, any non-stationarity of extreme value threshold should be examined and identified either before or alongside non-stationarity of GP parameters. Anderson et al. (2001) notes that the combination of non-stationary threshold and stationary GP shape and scale is sufficient for modelling a sample of significant wave heights in the North Sea. Physical and statistical intuition suggest, when considering a an extreme value model for a quantity such H_S, that non-stationary estimates should be sought in order for each of (a) extreme value threshold, (b) then GP scale, and (c) finally GP shape, and adopted only if justified statistically. Note however that physically plausible oceanographic examples corresponding, for instance, to stationary extreme value threshold but non-stationary GP parameters are also conceivable. In the current work, we assume effectively that any non-stationarity in extreme value threshold as already been identified perfectly and its effect removed from the sample cases considered.

It appears that specification of a piecewise model for the whole sample (e.g. , , ) incorporating the appropriate tail form, rather than a model for threshold exceedances in isolation is a promising route, since then threshold and tail parameters can be estimated together. However, even the simplest form of a whole sample model requires estimation of additional parameters for the model below the threshold; these parameters are also typically non-stationary with respect to covariates. Efficient and reliable non-stationary estimation using parameterisations and inference schemes similar to those presented here is key.




§ ACKNOWLEDGEMENT


We are grateful to Kathryn Turnbull of Lancaster University, UK for useful comments.






§ APPENDIX






§ MAXIMUM LIKELIHOOD ESTIMATION



For maximum likelihood estimation (MLE), we use a back-fitting (or iteratively re-weighted least-squares, IRLS) algorithm to estimate vectors of basis coefficients  (= ,) derived in <cit.>. For fixed value of smoothness parameter , we initialise coefficients to starting value 0, and then iterate the following step until convergence


    i+1 = ( i  + )^-1( i +  i i)


where


    = ∇_|   and    = -∇_∇_|


are derived at the end of this Appendix. This algorithm is similar to the mMALA algorithm (see below) used to generate proposals for the corresponding Metropolis-Hastings step in MCMC, in that both exploit first- and second-derivative information to move towards regions of high probability. The back-fitting iteration is of course deterministic, whereas the mMALA step is stochastic. In practice we use the expected values of  for ease of computation.




§ MCMC SAMPLING ALGORITHMS



Denoting the set of parameters to be estimated by = {,,,}, inference proceeds by sampling from the full conditional distributions k|,,k, for each parameter in  in turn, where = {,,,} is the set of fixed hyper-parameters for prior distributions. The form of the full conditional distribution varies depending on the type of parameter being estimated, as explained below.




 §.§ Full conditional distributions for basis coefficients



Vector  (= ,) has the following conditional distribution


    |,,,∝|,|


which is not available in closed form, and therefore cannot be sampled directly in a Gibbs step. Instead, we generate samples using the Metropolis-Hastings algorithm: given current state i, we propose new parameter  from proposal distribution
|i and evaluate the acceptance ratio


    ,i = |,,i|/i|,,|i


accepting the proposal with probability q = min(1,A), setting i+1 =. Otherwise we reject the proposal and set i+1 = i. As outlined in Section <ref> and detailed below, we consider two different methods for generating multivariate proposals for . In the first approach (referred to as MH), we make an entirely stochastic Gaussian random walk proposal using a fixed covariance matrix; in the second (referred to as mMALA), we make a proposal which is partly deterministic and partly stochastic, accounting for local curvature of the likelihood surface.

For Metropolis-Hastings (MH) inference, we generate Gaussian random walk proposals of the form


    = i + ( + κ_)^-1ϵ


where ϵ is a vector of independent standard Normal random variables, and the values of step size  and scale factor κ_ are adjusted to achieve reasonable acceptance rates of approximately 0.25. For inference using the Riemann manifold Metropolis-adjusted Langevin algorithm (mMALA, as implemented by ) we propose using derivatives of the target distribution at the current sample. This promotes proposals in regions of higher probability, at the additional computational cost of computing necessary derivatives and matrix inverses. At iteration i of the sampling algorithm, where the current sample of the coefficients is i, proposals are made as


    = i + ^2/2 -1i i + √(-1i) ϵ


where ϵ is a vector of independent standard Normal random variables,  is (adjustable) step size, and


    i = ∇_|_i   and   i = -∇_∇_|_i


are the negative gradient and negative Hessian of the log density, with


    = log|,,,   and   ∇_ = (∂/∂η1,…,∂/∂η)   .


Computation of likelihood derivatives is described at the end of this Appendix. In practice we use the expected values of i for ease of computation.




 §.§ Full conditional distributions for prior precisions



Prior precision parameter  (= ,) has the following conditional distribution


    |,,,∝||,   .


By construction, since the Gamma distribution is a conjugate prior for the precision of a Gaussian distribution, we know that the full conditional distribution is also Gamma, with updated parameters


    =  + /2   and    =  + 1/2   .






§ DERIVATIVES OF THE POSTERIOR DISTRIBUTION



Here we find the derivatives of the log posterior distribution, required for maximum likelihood and mMALA inference. The log likelihood of the observed data under the generalised Pareto distribution is


    | = ∑_i=1^[-log(i/1+i) - (1/i+1)log(1+i/i(1+i)i)]    for i≠ 0 
    ∑_i=1^[-log(i/1+i) - (1+i)i/i]    for i = 0    .


The log conditional distribution for the vector of basis coefficients  (= ,) is then the sum of this likelihood plus a contribution from the prior distribution


    = log(|,,,)
        = | - /2   .


We note the equivalence between this expression and the penalised (negative log) likelihood used for maximum likelihood inference. The gradient of the log conditional distribution for vector of coefficients  (= ,) is


    ∇_ = ∇_| -    .


Using the chain rule, the likelihood gradient can be computed as


    ∇_|    = ( ∇_ ())( ∇_|)
        = (∇_|)   .


The components of ∇_| are computed as


    ∂/∂i = 
    -1/i^2G_i(1-2i)(G_i-1) + 1/1+i + 1/ilog(G_i)     for i≠ 0 
    
    -i/i + 1/1+i    for i = 0


where G_i = 1 + i/i(1+i)i, and the components of ∇_| are


    ∂/∂i = 1/i(1 - (1/i+1)G_i-1/G_i)    for i≠ 0 
    1/i(1 - G_i-1/i)    for i = 0    .


Differentiating ∇_ (= ,) again gives the Hessian matrix


    ∇_∇_ = ∇_∇_| -    .


Applying the chain rule


    ∇_∇_| = (∇_∇_|)   .


Note that the components of ∇_| and (∇_∇_|) are computed separately for = and =. Further, the expected values of likelihood second derivatives with respect to  and  are


    -𝔼_Y [∂^2/∂i∂j|] = 1/(1+i)^2    for  i = j 
    
    0     for  i ≠ j


and


    -𝔼_Y [∂^2/∂i∂j|] = 1/^2(1+2i)    for  i = j 
    
    0     for  i ≠ j


such that Hessian matrices are diagonal. Moreover, the expectations of all of the cross derivatives ∂^2/∂i∂j| are zero, since estimates of  and  are asymptotically independent by construction (e.g. ).




elsarticle-harv


§ REFERENCES





