Extreme value statistics theory

$%%\DeclareMathOperator{\MPM}{MPM}$ $\newcommand{\MPM}{M\!P\!M}$

The theory for the calculation of the extreme value statistics results provided by OrcaFlex depends on which extreme value statistics distribution is chosen:

Rayleigh distribution

The derivation and use of spectral moments to fit the Rayleigh distribution is described in detail by Ochi (he calls it the "exact method"). For the Rayleigh distribution to be appropriate for the peak responses, the response variable must be a stationary Gaussian process.

Ochi shows that, under the Gaussian assumption, the most probable maximum value, $\MPM$, occurring in storm duration $T$ is given by \begin{equation} \MPM = \mu + \sigma (2 \ln{n})^{1/2} \end{equation} where $n = T/T_z$ is the number of peaks and $\mu$, $\sigma$ and $T_z$ are the mean, standard deviation and mean up-crossing period, respectively, of the time history being analysed. OrcaFlex reports this most probable extreme value for the specified storm duration.

Ochi goes on to show that, for low values of risk parameter $\alpha$, the extreme value which will be exceeded with probability $\alpha$ is \begin{equation} \label{Exceedence} \mu + \sigma \left\{ 2 \ln{\left(\frac{n}{\alpha}\right)} \right\}^{1/2} \end{equation} Formula (\ref{Exceedence}) is however only approximate, and is valid only for low values of $\alpha$. So OrcaFlex instead calculates the extreme value which will be exceeded with probability $\alpha$ using the alternative formula \begin{equation} \label{BetterExceedence} \mu + \sigma \left[ 2 \ln{\left\{ \frac{-n}{\ln{(1-\alpha)}} \right\}} \right]^{1/2} \end{equation} Ochi derives this more accurate formula (\ref{BetterExceedence}) using Cramér's approximation method in his paper on prediction of extreme values, and it gives good results over the whole range of values of risk factor $\alpha$.

Note that all the above formulae are for maxima, i.e. analysing the upper tail extremes. Corresponding formulae for minima (lower tail) are trivially obtained by replacing $\mu+\sigma(\ldots)$ with $\mu-\sigma(\ldots)$ in these formulae.

Weibull and generalised Pareto distributions

The Weibull and generalised Pareto (GPD) distributions are both fitted using the maximum likelihood method. Both distributions are fitted to extremes of the data that are selected using the peaks-over-threshold method with (optional) declustering. They differ only in the actual statistical distribution that is fitted and then used to predict the extreme values.

Maximum likelihood

Maximum likelihood estimation is a standard general-purpose statistical technique for fitting a given parametric distribution to an arbitrary set of data. It fits the parameters of the distribution to the data by calculating the parameters that make the observed data most likely under the chosen distribution. The technique is widely described in statistical texts; in particular the book by Coles is particularly relevant, since he describes (section 4.3.2, parameter estimation) its application to the GPD. Coles also documents the motivation for, and implementation of, the peaks-over-threshold method (chapter 4, threshold models) and declustering (5.3, modelling stationary series).

Threshold models

When the object of primary interest is extrapolation into the tails of the data, it is safest to use a model that is fitted only to the tails of the data. This is because models that are fitted to the entire dataset tend to be driven by features in the body of the data, which may not be relevant to the tail behaviour. Mathematical theory (e.g. Coles, Theorem 4.1) tells us that as the size of the dataset and the threshold for fitting increase, the peaks over the threshold converge in distribution to a generalised Pareto distribution.

Declustering

An important assumption underpinning the maximum likelihood method is that the data are independent. However, this is often not the case. Consider daily temperature for example, where one cold day is likely to be followed by another: successive daily values are not independent, but values a week or a month apart might well be considered independent. For OrcaFlex time history results this is even more important, since unless the log interval is extremely long, then successive time history values will certainly not be independent.

We can get much closer to independence by sub-sampling the data, ensuring that the points we choose are sufficiently far apart as to be approximately independent. The usual way to do this is declustering: OrcaFlex uses either runs declustering as illustrated by Coles (5.3.3, Wooster temperature series), or defines a cluster to be the set of values between successive up-crossings of the mean value.

The complete model-fitting process is then summarised as:

  1. Use declustering to identify clusters of exceedences
  2. Determine the maximum excess within each cluster
  3. Assuming independence between these maxima, use maximum likelihood to fit the chosen distribution to them.

Weibull distribution

OrcaFlex fits the two-parameter Weibull distribution with distribution function \begin{equation} F_{\mathrm{W}_2}(y) = 1 - \exp\left\{-\left(\frac{y}{\sigma}\right)^\xi\right\} \end{equation} to the declustered excesses over threshold, $y$, where

$\sigma=$ Weibull distribution scale parameter

$\xi=$ Weibull distribution shape parameter.

This is equivalent to applying the three-parameter Weibull distribution \begin{equation} F_{\mathrm{W}_3}(y) = 1 - \exp\left[-\left\{\frac{(x-\mu)}{\sigma}\right\}^\xi\right] \end{equation} to the original values $x$ unadjusted for the threshold (i.e. $x=\mu+y$), where the location parameter $\mu$ is the user-specified threshold value and is not fitted by the maximum likelihood calculation.

Generalized Pareto distribution

The generalized Pareto distribution function, for the excesses $y$ above threshold, is \begin{equation} F_\textrm{GPD}(y) = \begin{cases} 1 - \left(1+\xi\frac{y}{\sigma}\right)_+^{-1/\xi} & \text{for $\xi\neq0$, where } \left(1+\xi\frac{y}{\sigma}\right)_+ = \max\{1+\xi\frac{y}{\sigma}, 0\} \\ 1 - \exp\left(\frac{-y}{\sigma}\right) & \text{for $\xi=0$} \end{cases} \end{equation} where

$\sigma=$ GPD scale parameter

$\xi=$ GPD shape parameter.