

\section{Introduction}

A common approach to practical extreme value analysis is to use the GEV form of the asymptotic extreme value distributions to fit to the observed extreme values. The typical data that are used in this process are the extreme values observed over specified periods of time, e.g. over one year periods. Such a procedure, extracting only extremes over blocks of data, would immediately appear to be wasteful, since potentially very useful data might be discarded. The Peaks-Over-Threshold method, or simply the POT method, represents an approach to extreme value analysis that tries to avoid this waste of data by considering all data that exceed a prescribed high threshold. It will be shown that also for the POT method there are three limiting forms of the exceedance probability distributions corresponding to The Extremal Types Theorem. The Generalized Pareto (GP) distribution will take the place of the GEV distribution.   



\section{The Peaks-Over-Threshold Method}\index{POT method}\index{peaks-over-threshold method}


The basic assumption of the POT method is that the observed time series $x_1, x_2, \ldots$ are the outcomes of a sequence of independent and identically distributed (iid) random variables $X_1, X_2, \ldots$ with common distribution function $F_X(x)$, which belongs to the domain of attraction of one of the extreme value distributions. Instead of extracting the extreme observations over blocks of data, the focus of the POT method is on all data that exceed a given high threshold. If the chosen threshold $u$ is high enough, it would be natural to consider also the data exceeding $u$ as extremes. Denoting an arbitrary term of the $X_i$ sequence by $X$, the exceedance probability essential to the POT method, is the following conditional probability:

\begin{equation}\label{eq:pot.1}
\Prob(X > u + y \mid X > u ) = \frac{1-F(u+y)}{1-F(u)}, \quad  y > 0,
\end{equation}
where $y$ denotes the size of the exceedance above the threshold $u$. 
 
If the parent distribution $F$ were fully known, the distribution of the threshold exceedances in Eq. \eqref{eq:pot.1} would also be known. However, in practical applications this is rarely the case. Typically, only estimated approximations to the parent distribution based on the observed data would be available. A consequence of this is that the estimates of $F(x)$ for $x \le u$ in general becomes very uncertain for large to extreme values of $u$. This makes the direct application of Eq. \eqref{eq:pot.1} useless for practical purposes. Hence, limiting forms that would parallel the GEV distributions for extremes are sought. 

The main result needed is expressed by the following theorem:

\bigskip
Let $X_1, X_2, \ldots$ be a sequence of iid random variables with common distribution function $F$, and let $M_n = \max\{X_1, \ldots, X_n\}$. Assume that the conditions of the Extremeal Types Theorem are satisfied, so that for large $n$,
\begin{equation}\label{eq:pot.2}
\Prob(M_n \le z ) \approx G(z), \nonumber
\end{equation}
where
\begin{equation}\label{eq:pot.3}
G(x) = \exp \left\{ - \left( 1 + \gamma \Big( \frac{x -
	\mu}{\sigma} \Big) \right)^{- 1/ \gamma} \right\} \,,
\end{equation}
defined on the set $\{x \,: \, 1 + \gamma \big( (x -
\mu) / \sigma \big) > 0 \}$, where the parameters satisfy $- \infty < \mu < \infty$, $\sigma >
0$, $- \infty < \gamma < \infty$.
Then, for large enough $u$,
\begin{equation}\label{eq:pot.4}
H(y) = \Prob(X \le u + y \mid X > u ) \approx 1  - \Big(1 + \gamma \frac{ y}{\tilde{\sigma}} \Big)^{- 1/\gamma},
\end{equation} 
defined on $\{y: y>0 \:\: \mbox{and} \:\: (1+ \gamma y / \tilde{\sigma}) > 0 \}$, where $\tilde{\sigma} = \sigma + \gamma (u - \mu)$.
    
Note that the special case $\gamma = 0$ has to be interpreted as a limit, in analogy with the GEV case $\gamma = 0$, viz.
\begin{equation}\label{eq:pot.5}
H(y) = 1  - \exp \left( - \frac{y}{\tilde{\sigma}} \right), \quad y > 0,
\end{equation} 
which corresponds to an exponential distribution with parameter $1 / \tilde{\sigma}$.   

\bigskip
The family of distribution functions defined by Eq. \eqref{eq:pot.4} is called the Generalized Pareto (GP) distributions\index{generalized Pareto distribution}. The theorem above implies that if the block extremes have the limiting distribution $G$, then threshold excesses have a corresponding limiting distribution within the GP  
family of distributions. And it follows from this theorem that the parameters of the GP distribution\index{GP distribution} of threshold excesses are uniquely determined by those of the associated GEV distribution of block extremes. In particular, the shape parameter $\gamma$ of the GP distribution is identical to that of the corresponding GEV distribution. It is important to note, however, that this statement is valid only asymptotically. The practical significance of this is that, in general, the $\gamma$-parameter estimated on the basis of a finite set of data rarely equals that of the correct asymptotic extreme value distribution. Also relevant in this context, is the observation made by \cite{fiti28} that the best fit to a finite set of extreme value data from a stationary Gaussian process is provided by a Type III distribution and not by the asymptotically correct Type I distribution. 

A full proof of the theorem above will not be given here, but a sketch of the proof may serve to illustrate the main ideas. A more precise argument is provided by \cite{leliro83}.

For large $n$, and for suitable values of the argument $z$,
\begin{equation}\label{eq:pot.6}
\Prob(M_n \le z ) = F^n(z) \approx \exp \left\{ - \left( 1 + \gamma \Big( \frac{z -
	\mu}{\sigma} \Big) \right)^{- 1/ \gamma} \right\} \,, \nonumber
\end{equation}
Hence,
\begin{equation}\label{eq:pot.7}
n \log F(z) \approx - \left( 1 + \gamma \Big( \frac{z -
	\mu}{\sigma} \Big) \right)^{- 1/ \gamma}  \,, 
\end{equation}
By a Taylor expansion, it is found that $\log (1-\varepsilon) \approx - \varepsilon$ for small positive $\varepsilon$. Consequently, for large values of $z$, it follows that $\log F(z) \approx - (1-F(z))$. Substituting into Eq. \eqref{eq:pot.7}, followed by a rearrangement, leads to
\begin{equation}\label{eq:pot.8}
1 -  F(u) \approx  \frac{1}{n} \left( 1 + \gamma \Big( \frac{u -
	\mu}{\sigma} \Big) \right)^{- 1/ \gamma}  \,, \nonumber
\end{equation}
for large values of $u$. Similarly, for $y > 0$,
\begin{equation}\label{eq:pot.9}
1 -  F(u + y) \approx  \frac{1}{n} \left( 1 + \gamma \Big( \frac{u + y -
	\mu}{\sigma} \Big) \right)^{- 1/ \gamma}  \,, \nonumber
\end{equation}
Hence, it follows that,
\begin{align}\label{eq:pot.10}
\Prob \Big( X > u + y \mid X > u \Big) &\approx \left( \frac{1 + \gamma (u + y -
	\mu) /\sigma }{1 + \gamma (u - \mu) /\sigma}  \right)^{- 1/ \gamma} \nonumber \\
&= \left( 1 + \gamma \frac{y}{\tilde{\sigma}}  \right)^{- 1/ \gamma} 
\end{align}
where $\tilde{\sigma} = \sigma + \gamma (u - \mu)$, as required by Eq. \eqref{eq:pot.4}.





\section{Threshold Selection} \index{POT method!threshold selection}

Selection of an appropriate threshold to use when applying the POT method on measured data is an essential ingredient of this method. Its importance is reflected by the fact that the predicted extreme values will in many cases show a significant dependence on this choice. Unfortunately, in practice, there are no fully reliable methods to guide the selection of an appropriate threshold. In spite of this, a couple of methods have seen extensive use. These two methods will be discussed in this section.

The practical application of the POT method would typically proceed as follows. The observed data $x_1, x_2, \ldots, x_n$ are assumed to be outcomes of independent and identically distributed random variables $X_1, X_2, \ldots, X_n$, which, of course, would require some degree of justification. The extreme events to be used are the observed exceedances above a high threshold $u$, that is, the data $x_i : x_i > u$. Denote these exceedances as $x_{(1)}, x_{(2)}, \ldots, x_{(k)}$, and define the threshold excesses by $y_j = x_{(j)} - u$ for $j = 1, \ldots, k$. According to our main result in the previous section, the excess data $y_j$ may be regarded as independent realizations of a random variable whose distribution function is approximately of the GP type. By fitting the excess data to the GP distribution, an approximate distribution of the excess variable is obtained. 

In the choice of threshold, there are clearly two conflicting issues that have to be dealt with. Choosing a high threshold is desirable from the point of view of not violating too much the asymptotic basis for the theory. However, if too few data are retained, there will be high uncertainty in every estimate. On the other hand, if the threshold is too low, the apparent uncertainty of estimates may be reduced, but the assumption that the excess values follow a GP distribution may be seriously in error. Hence, in practice, there is a need to balance the two issues. There seems to be two procedures that are used for this purpose. One is carried out prior to model estimation by investigating how a diagnostic follows an expected pattern. The other is an assessment of the stability of parameter estimates based on the fitting of models for a range of different thresholds.

The first of the two methods uses the expected value of a GP variable $Y$ with distribution function 
\begin{align}\label{eq:pot.11}
H(y) =  1  - \Big(1 + \gamma  y/{\sigma} \Big)^{- 1/\gamma},
\end{align}
The expected value is then,  
\begin{align}\label{eq:pot.12}
\EE[Y] = \frac{\sigma}{1 - \gamma},
\end{align}
provided $\gamma < 1$, otherwise the mean value is infinite. Assume that the GP distribution is a valid model for the excesses of a fixed threshold $u_0$. Let $X$ denote an arbitrary term among $X_1, X_2, \ldots, X_n$. Then, by \eqref{eq:pot.12},
\begin{align}\label{eq:pot.13}
\EE[X - u_0 \mid X > u_0 ] = \frac{\sigma_{u_0}}{1 - \gamma},
\end{align}
provided $\gamma < 1$, where the convention of using $\sigma_{u}$ to denote the scale parameter corresponding to excesses of the threshold $u$ has been adopted. If the GP model holds for the threshold $u_0$, by necessity it also holds for any $u > u_0$. Hence, for $u > u_0$, it is obtained that,
\begin{align}\label{eq:pot.14}
\EE[X - u \mid X > u ] = \frac{\sigma_{u}}{1 - \gamma} = \frac{\sigma_{u_0} + \gamma (u - u_0) }{1 - \gamma},
\end{align}
where the last equation follows from the main result of the previous section. This tells us that for $u > u_0$, $\EE[X - u \mid X > u ]$ is a linear function of $u$. This can then be verified by calculating the mean value of the excesses for a range of thresholds. Thus, a plot of the following empirical point estimates,
\begin{align}\label{eq:pot.15}
\left\{ \Big(u, \frac{1}{n_u} \sum_{i=1}^{n_u} (x_{(i)} - u) \Big) \: : \: u < x_{\max}\right\},
\end{align}
where $x_{(1)}, x_{(2)}, \ldots, x_{(n_u)}$ consist of the $n_u$ observations that exceed $u$, and $x_{\max}$ denotes the largest observation, should be an approximately straight line for a range of $u$-values where the GP model is applicable. This empirical graph is commonly called the {\it mean residual life plot}\index{mean residual life plot}. Confidence intervals can be added to the graph by using the approximate normality of the sample means. We shall have occasion to illustrate this diagnostic for the applicability of the POT method for real-life data, and it will be clear that it is not always of much help. 




\section{Return Periods} \index{return period} 
The return period $R$ of a given wind speed, in years, is defined
as the inverse of the probability that the specified wind speed
will be exceeded in any one year. If $\lambda$  denotes the mean
crossing rate of the threshold $u$ per year (i.e., the average
number of data points above the threshold $u$ per year), the
return period $R$ of the value of $X$ corresponding to the level
$x_R = u + y_R$ is given by the relation
\begin{equation}
R = \frac{1}{\lambda \, \Prob(X > x_R)} = \frac{1}{\lambda \,
	\Prob(Y > y_R)}.
\end{equation}
Hence, it follows that
\begin{equation}
\Prob(Y \le y_R) =  1 - 1/(\lambda R).
\end{equation}
Invoking Eq.~(\ref{eq:pot.11}) for $\gamma \not= 0$ leads to the result
\begin{equation}
x_R  = u  + \sigma [ (\lambda R)^\gamma - 1]/\gamma.
\end{equation}
Similarly, for $\gamma=0$, it is found that
\begin{equation}
x_R  = u +  \sigma \ln (\lambda R),
\end{equation}
where $u$ is the threshold used in the estimation of $\gamma$ and $\sigma$. 
A discussion of how to estimate confidence intervals on $x_R$ using profile likelihood methods, is provided by \cite{cole01}. In this book the focus is on the bootstrap method, cf. Section~\ref{sec:bootstrapping!POT}.


\section{Parameter Estimation for the GP Distributions}\label{sec:pegp} \index{parameter estimation!GP distribution}\index{POT method!parameter estimation}

There is a range of possible estimation methods available for the parameters of the GP distribution. In our own work, we have mostly used three methods that serve the purpose reasonably well when the shape parameter does not deviate too much from zero ($\gamma < 0.5$): The Dekkers-Einmahl-de Haan estimators, the moment estimators and the maximum likelihood estimators. These three methods will be discussed below. Alternative methods include the Hill estimators and the probability weighted moments (PWM) estimators, cf. \cite{begosete04}. PWM estimators has been popular in flood frequency analysis \citep{howawo85}. A further development are provided by the L-moment estimators \citep{howa97}. It has been shown that L-moments have some issues related to their lack of sensitivity to the tail behaviour of the underlying statistical distribution. This is clearly of importance for properly representing extreme value distributions, cf. \cite{wima13}.


\subsection{Dekkers-Einmahl-de Haan estimators}\index{estimator!de Haan}\index{POT method!de Haan estimator}

Let $n$ denote the total number of data points, while the number
of observations above the threshold value $u$ is denoted by $k$.
The threshold $u$ then represents the $(k+1)$th highest data
point(s). An estimate for $\lambda$ is $\hat{\lambda} =
k/n_{yrs}$,  where $n_{yrs}$ denotes the length of the record in
years. The highest, second highest, ... , $k$th highest,
$(k+1)$th highest variates are denoted by $X_{n}^*$, $X_{n-1}^*$,
...,$ X_{n-k+1}^*$, $X_{n-k}^*=u$, respectively.

The parameter estimators proposed by \cite{deeiha89} and \cite{deha94} are based on
the following two quantities:
\begin{equation}
H_{k,n}  =  \frac{1}{k} \sum _{i=0}^{k-1} \{\ln(X_{n-i}^*) -
\ln(X_{n-k}^*)\}
\end{equation}
and
\begin{equation}
H_{k,n}^{(2)}  =  \frac{1}{k} \sum _{i=0}^{k-1} \{\ln(X_{n-i}^*) -
\ln(X_{n-k}^*)\}^2\,.
\end{equation}

Estimators for $\sigma$ and $\gamma$ are then given by the relations
\begin{equation}
\hat{\sigma} = \rho \, X_{n-k}^* \, H_{k,n} = \rho \, u \, H_{k,n}
\end{equation}
and
\begin{equation}
\hat{\gamma} = H_{k,n} + 1 - \frac{1}{2}\left\{ 1 - \frac{(H_{k,n})^2}{
	H_{k,n}^{(2)}}\right\}^{-1}
\end{equation}
where $\rho = 1$ if $\hat{\gamma} \ge 0$, while   $\rho = 1 -\hat{\gamma} $
if $\hat{\gamma} < 0$. These two estimators will be referred to as the Dekkers-Einmahl-de Haan estimators.

Subject to general conditions on the underlying probability law,
\cite{deeiha89} and \cite{deha94} showed that $\hat{\gamma} \to \gamma$ and $\hat{\sigma} \to \sigma$ as
$n \to \infty$ (in probability).

Closely related to the Dekkers-Einmahl-de Haan estimators are the Hill estimators.
Their application to the problem of estimating extreme wind speeds
was investigated by \cite{nacl99}. Their conclusion was that the
Hill estimators lead to results that are quite similar to those
provided by the Dekkers-Einmahl-de Haan estimators. Because the Hill estimators
require considerably higher numerical efforts than the Dekkers-Einmahl-de Haan
estimators, and rarely provide significantly better results, the Hill estimators were excluded from the present discussion. The interested reader is referred to \cite{nacl99,begosete04} for details.




\subsection{Moment estimators}\index{estimator!moment}\index{POT method!moment estimator}

In terms of the mean value $\EE(Y)$ and the standard deviation
$s(Y)$ of the exceedance variate $Y$, it can be shown that
\citep{howa87}
\begin{equation}
\sigma = \frac{1}{2} \EE(Y) \{ 1 + [\EE(Y) / s(Y) ]^2 \}
\end{equation}
and
\begin{equation}
\gamma = \frac{1}{2} \{ 1 - [\EE(Y) / s(Y) ]^2 \}\,.
\end{equation}

Hence, empirical estimates of the first two moments of $Y$ provide
estimates of $\sigma$ and $\gamma$. The resulting
estimators are referred to as the moment estimators for $\sigma$ and $\gamma$.




\subsection{Maximum likelihood estimators}\index{estimator!maximum likelihood}\index{POT method!maximum likelihood estimator}

The maximum likelihood estimators (MLEs) are often preferred due to their asymptotic efficiency. 
Let $y_1, \ldots, y_k$ denote the observed sample of exceedances above the threshold $u$ from an observed sample of peak values $x_1, \ldots, x_n$. 
The log-likelihood function $\ell(\sigma, \gamma | y_1, \ldots, y_k)$ for the sample 
$y_1, \ldots, y_k$ is given by
\begin{equation}
\ell(\sigma, \gamma | y_1, \ldots, y_k) = - k \ln \sigma - \left( \frac{1}{\sigma} +1 \right) \sum_{i=1}^k \ln \Big( 1 + \gamma \frac{y_i}{\sigma} \Big)\,, 
\end{equation}
provided $( 1 + \gamma y_i/ \sigma) > 0$ for $i = 1, \ldots, k$. If $\gamma=0$, the log-likelihood assumes the form
\begin{equation}
\ell(\sigma,0 | y_1, \ldots, y_k) = - k \ln \sigma - \frac{1}{\sigma} \sum_{i=1}^k y_i\,. 
\end{equation}
The MLE $\hat{\sigma}$ and $\hat{\gamma}$ are obtained by maximizing $\ell(\sigma, \gamma | y_1, \ldots, y_k)$ with respect to $\sigma$ and $\gamma$. These values are found by numerical methods, except for the special case $\gamma=0$, for which a simple, closed-form solution exists for $\hat{\sigma}$. It is given as $\hat{\sigma} = \sum_{i=1}^k y_i/k$.




\section{Model Validation}\index{POT method!model validation}


Probability (PP) and quantile (QQ) plots may be useful tools to get a grip on the suitability of the GP model for the excesses beyond a chosen threshold $u$, cf. \cite{begosete04}. Denote the threshold excesses by $y_1, y_2, \ldots, y_k$ and the estimated GP model by $\hat{H}$. Ordering the threshold excesses by increasing magnitude: $y_{(1)} \le y_{(2)} \le \ldots \le y_{(k)}$, the PP plot consists of the point graph,
\begin{equation}
\big( i/ (k+1) ,\hat{H}( y_{(i)}) \big) , \: i = 1, \ldots, k,
\end{equation}
where
\begin{align}\label{eq:pot.xx}
\hat{H}(y) =  1  - \Big(1 + \hat{\gamma}  y/{\hat{\sigma}} \Big)^{- 1/\hat{\gamma}},
\end{align}
provided $\hat{\gamma} \ne 0$. If $\hat{\gamma} = 0$, the plot is constructed using the exponential distribution.

Assuming again that $\hat{\gamma} \ne 0$, the QQ plot is traced out by the point graph, 
\begin{equation}
\big( \hat{H}^{-1}(i/ (k+1)) , y_{(i)} \big) , \: i = 1, \ldots, k,
\end{equation}
where 
\begin{equation}
\hat{H}^{-1}(y) = u + \frac{\hat{\sigma}}{\hat{\gamma}} \Big(y^{- \hat{\gamma}} - 1 \Big),
\end{equation}
 
If the GP model is a reasonable model for the distribution of the excesses of $u$, both the PP and the QQ plot should follow approximately a straight line. These procedures are discussed more thoroughly in Chapter~\ref{chap:mc}. 





\section{Estimating Confidence Intervals by Bootstrapping}\label{sec:bootstrapping!POT}\index{POT method!confidence interval by bootstrapping} \index{bootstrapping}\index{confidence interval}

The principle of the bootstrapping method was briefly explained in the previous chapter, cf. Section~\ref{sec:bootstrapping!GEV}. 
In this section on the POT method, only the nonparametric bootstrap
is used. The goal of this method is to estimate some statistical quantity $V$ given by the unknown distribution function on the basis of an observed sample $\mathbf{y} = (y_1, y_2, \ldots, y_n)$, which is a
sample or vector consisting of $n$ independent observations of a
random variable $Y$. Let $\hat{V}$ denote the estimate of $V$ based on the given sample.
The nonparametric bootstrapping\index{bootstrapping!nonparametric} technique for estimating
confidence intervals on $V$ is based on
resampling (with replacement) from the empirical distribution function (EDF) provided by the
observed sample $\mathbf{y}$, cf. Section~\ref{sec:bootstrapping!GEV}.

This is done as follows: the EDF gives rise to an empirical random
variable $Y^*$. $\ell$ bootstrap samples $\mathbf{y}^*_j$,
$j=1, \ldots, \ell$, with $n$ independent observations of $Y^*$ in
each sample are now generated. Each sample
$\mathbf{y}^*_j$ gives rise to an estimate $V^*_j$ of
$V$.

Simple estimates for confidence intervals on $V$ are
derived by calculating the sample standard deviation $s_V^*$:
\begin{equation}
s_V^*= \sqrt{\frac{1}{\ell-1} \sum_{j=1}^\ell  ( V^*_j - \bar{V}^*)^2},
\end{equation}
where $\bar{V}^* = (1/\ell)\,\sum_{j=1}^\ell V_j^*$. An approximate
confidence interval at level $1-q$ is then obtained as
\begin{equation}
(\,\hat{V} - w_{q/2}\,s_V^*\,,\, \hat{V} + w_{q/2}\,s_V^*\,),
\end{equation}
where $w_{q/2}$ denotes the $100(1-q/2)$\% standard normal
fractile. To get stable results for this procedure, usually 20-30
bootstrap samples are sufficient. To avoid making the assumption
that the bootstrap estimates $V^*_j$ are normally distributed, the true
distribution may be approximated by generating a large number of
bootstrap samples, usually several thousand are needed, especially
for small values of $q$. If $\ell$ samples were generated, the
$V_j^*$ are rearranged in increasing order $V^*_{(1)} \le V^*_{(2)} \le \ldots \le V^*_{(\ell)}$. A $100(1-q)$\%
confidence interval for $V$ is then,
\begin{equation}
(\,V_{(L)}^*\,,\, V_{(M)}^*\,)\,,
\end{equation}
where $(L) =[q\ell/2]$ and $(M) =[(1-q/2)\ell]$ ($[a]$ means the integer part
of $a$). \cite{dahi97} describe possible improvements of such estimates. However, they are not discussed here as they are considered to be of less practical interest.


