


\section{Introduction}

Classical extreme value statistics is concerned with the distributional properties of the maximum of a number of independent and identically distributed (iid) random variables when the number of variables becomes large. A partial result was obtained by \cite{frec27}, while \cite{fiti28} discovered that there are three types of possible limiting or asymptotic distributions, which are now contained in the Extremal Types Theorem, which is discussed in the next section. These three asymptotic distributions are typically referred to as the Gumbel, Fr\'{e}chet and Weibull distributions. It is also common practice to refer to them as Type I, Type II and Type III, in the same order. Important contributions to this theory were later made by \cite{gned43}, \cite{gumb58} and \cite{deha70}.


\section{The Asymptotic Limits of Extreme Value Distributions} \index{extreme value distribution!asymptotic limits}\label{sec:alevd}


The classical extreme value theory starts by looking at a sequence of independent and identically
distributed (iid) random variables $X_1, X_2, \ldots$ with common
distribution function $F_X(x)$. The extreme value of a finite
number $X_1, \ldots, X_n$ is then $M_n = \max\{X_1, \ldots,
X_n\}$. The distribution of $M_n$ can be easily derived as
\begin{align}\label{eq:8.33}
F_{M_n}(x) &= \Prob(M_n \le x) = \Prob(X_1 \le x, \ldots, X_n \le
x ) \nonumber \\
&= \Prob(X_1 \le x) \cdot \ldots \cdot \Prob(X_n \le x) = \big(
F_X(x)\big)^n.
\end{align}
This relation is not very helpful in practice, because in most cases the distribution function $F_X(x)$ is not known exactly. Therefore, it would have to be estimated from recorded data. However, small discrepancies in the estimates of $F_X(x)$ can lead to substantial
discrepancies, in a relative sense, in the values of $\big( F_X(x)\big)^n$ for large values of $n$. In classical
extreme value theory, one proceeds by studying the behaviour of
$\big( F_X(x)\big)^n$ as $n \to \infty$, but with a twist.
Obviously, for any $x$ such that $F_X(x) < 1$, $\big(
F_X(x)\big)^n \to 0$ as $n \to \infty$. This necessitates a
rescaling. Specifically, instead of studying $M_n$, one introduces
a renormalized version of $M_n$:
\begin{align}\label{eq:8.34}
M_n^* = \frac{M_n - b_n}{a_n}
\end{align}
for suitable sequences of constants $a_n > 0$ and $b_n$ that are
chosen to stabilize the location and scale of $M_n^*$ as $n \to
\infty$. It is then proven that there are, in fact, only three
types of limiting distributions for this renormalized $M_n^*$.
This is the famous Extremal Types Theorem \citep{leliro83}, \index{extremal types theorem} which can be expressed
as follows.

\bigskip
If there exist sequences of constants $a_n > 0$ and $b_n$ such
that
\begin{align}\label{eq:8.35}
\Prob \Big( \frac{M_n - b_n}{a_n} \le x \Big) \to G(x) \,, \:\: n
\to \infty,
\end{align}
where $G(x)$ is a nondegenerate distribution function, then
$G(x)$ belongs to one of the following three families,
\begin{align}\label{eq:8.36}
\mbox{I} \qquad G(x) = \exp \left\{ - \exp \left[ - \Big( \frac{x
- b}{a} \Big) \right] \right\} \,, \:\: - \infty < x < \infty\,;
\end{align}
\begin{align}\label{eq:8.37}
\mbox{II} \qquad G(x) = \left\{
\begin{array}{l@{\quad , \quad}l} 0 & x \le b \,, \\
\exp \left\{ - \Big( \frac{x - b}{a} \Big)^{- c} \right\} & x
> b \,;
\end{array} \right.
\end{align}
\begin{align}\label{eq:8.38}
\mbox{III} \qquad G(x) = \left\{
\begin{array}{l@{\quad , \quad}l} \exp \left\{ - \Big( \frac{b - x}{a} \Big)^{ c} \right\} & x
< b\,, \\
1 & x \ge b \,;
\end{array} \right.
\end{align}
for parameters $a > 0$, $b$ and for families II and III, $c > 0$. 

\bigskip
These three types of extreme value distributions are also
commonly referred to as Gumbel, Fr\'{e}chet, and Weibull,
respectively.\index{Gumbel distribution}\index{Fr\'{e}chet distribution}\index{Weibull distribution}  Note that the Weibull distribution given here is not the same as the commonly known Weibull
distribution, which corresponds to the type III extreme value
distribution for minima. Also, carefully note that even if the
Weibull distribution is the only type of extreme value
distribution with a finite upper limit on its values, this does
not mean that extremes of limited data must follow this
distribution. For such data, it may very well happen that the
rescaling constants $a_n \to 0$ as $n$ increases. Hence, even the Gumbel distribution may be the appropriate asymptotic limit for the extreme values of bounded data. 

It may be verified that it is, in fact, possible to express all
three types of extreme value distributions in a common form, which
is known as the generalized extreme value (GEV)\index{GEV distribution} distribution\index{generalized extreme value distribution}. This is
achieved as follows:
\begin{align}\label{eq:8.39}
G(x) = G(x; \mu, \sigma, \gamma) = \exp \left\{ - \left[ 1 + \gamma \Big( \frac{x -
\mu}{\sigma} \Big) \right]^{- 1/ \gamma} \right\} \,,
\end{align}
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$. This
distribution has three parameters: a location parameter $\mu$, a
scale parameter $\sigma$, and a shape parameter $\gamma$. The type
II distributions correspond to $\gamma > 0$, while type III
corresponds to $\gamma < 0$. The case $\gamma = 0$ must be
interpreted as a limiting case when $\gamma \to 0$, which leads to
the Gumbel distribution:
\begin{align}\label{eq:8.40}
G(x) = \exp \left\{ - \exp \left[ -\Big( \frac{x - \mu}{\sigma}
\Big) \right] \right\} \,, \:\: - \infty < x < \infty.
\end{align}
The statistical moments of the GEV distributions can now be calculated based on the explicit formulas of Eqs.~\eqref{eq:8.39} and \eqref{eq:8.40}. Denoting the random variable determined by a GEV distribution  by $M$, its first two moments are,
\begin{align}\label{eq:8.40a}
\EE(M) = \mu + (e_1 - 1) \frac{\sigma}{\gamma}  \,, \:\: \gamma \ne 0, \gamma < 1,
\end{align}
and
\begin{align}\label{eq:8.40b}
\Var(M) = \left(e_{2}-e_{1}^{2}\right) \frac{\sigma^{2}}{\gamma^{2}}  \,, \:\: \gamma \ne 0, \gamma < 1/2,
\end{align}
where $e_k = \Gamma (1 - k \gamma)$, $k = 1, 2$, and $\Gamma(\cdot)$ is the gamma function.
For $\gamma \ge 1$, $\EE(M) = \infty$, while $\EE(M) = \mu + \lambda_E \, \sigma$ when $\gamma = 0$, that is, for the Gumbel case. Here, $\lambda_E = 0.5772..$ denotes Euler's constant. For $\gamma \ge 1/2$, $\Var(M) = \infty$, while $\Var(M) = \sigma^2 \pi^2 /6$, when $\gamma = 0$.


For statistical inference on experimental data, the unified form expressed by Eq.~\eqref{eq:8.39} has the advantage that the data themselves determine which type of
distribution is appropriate, thereby avoiding a prior subjective
judgment about any specific tail behaviour. The uncertainty in the
estimated value of $\gamma$ is also a reflection of the
uncertainty about the correct distribution for the data. Unfortunately, in practice, it may very well happen that the uncertainty in $\gamma$ may cover all three types of extreme value distribution, which would necessitate a more careful analysis of the data. Also note that the data used for estimation purposes are never truly asymptotic, thereby introducing additional uncertainty when trying to identify the correct asymptotic distribution. Since the results of extrapolation to determine long return period design values may depend very much on the asymptotic extreme value distribution used, identifying the correct one is clearly important in such cases. 


\section{The Block Maxima Method}\label{sec:bmm} \index{block maxima method}\index{GEV distribution!block maxima method}


In practical application of the GEV distributions to a long time series of observed data it is assumed that the maximum observation of a reasonably large chunk of the time series follows a GEV distribution. This is recognized by observing that from \eqref{eq:8.35} we would assume that for large $n$,
\begin{align}\label{eq:8.41}
\Prob \Big( \frac{M_n - b_n}{a_n} \le x \Big) \approx G(x) \,,
\end{align}
But this may be rewritten as ($y = a_n x + b_n$),
\begin{align}\label{eq:8.42}
\Prob \big( M_n \le y \big) \approx G \Big(\frac{y - b_n}{a_n} \Big) = G^*(y) \,,
\end{align}
where $G^*$ is also a member of the GEV family of distributions. Hence, if the main theorem applies, that is, by \eqref{eq:8.35}, $M_n^* = (M_n - b_n)/ a_n$ approximately follows a GEV distribution, then $M_n$ itself will approximately follow a GEV distribution, but with different parameters. Anyway, in practice, it is the parameters of $G^*$ that would be of most interest.

This leads to the following approach, which is often referred to as the block maxima method. Assume that a sequence of independent observations $x_1, x_2, \ldots$ from a stationary time series is long enough to allow segmenting it into blocks of data of length $n$, for some large value of $n$, generating a series of $m$ block maxima, $M_{n,1}, \ldots, M_{n,m}$, say, to which a GEV distribution is tentatively fitted. A typical application of the block maxima method would be to yearly extreme value observations of an environmental parameter, e.g. wind speed. In such a case, it is also often referred to as an annual maxima method. There is a practical argument behind extracting the maximum over the period of one year, because by choosing shorter periods, the assumption that the sampled maxima are outcomes of a common distribution would easily be violated due to seasonal variations. Still, of course, the underlying assumption that the block maxima are extracted from a set of iid random variables is clearly violated. Fortunately, by experience, this does not seem to pose a serious obstacle to the practical use of the block maxima method.  

A quantity of specific interest in applications is the return period level $x_p$, where $G(x_p) = 1 - p$. For the annual maxima method, $x_p$ has a return period of $1/p$ years. That is, $x_p$ would be exceeded on the average every $1/p$ years. Inverting \eqref{eq:8.39}, it is found that for $\gamma \ne 0$,
\begin{align}\label{eq:8.43}
x_p = \mu - (\sigma/\gamma) \big[ 1 - \big( - \log(1-p) \big)^{-\gamma} \big] \,, 
\end{align}
while for $\gamma = 0$,
\begin{align}\label{eq:8.44}
x_p = \mu - \sigma \log \big( - \log(1-p) \big) \,, 
\end{align}   

\cite{cole01} discusses how to estimate confidence intervals on $x_p$ using profile likelihood methods, which seem to provide reasonable accuracy. In this book the focus is on the bootstrap method, cf. Section~\ref{sec:bootstrapping!GEV}.
 



\section{Outline Proof of the Extremal Types Theorem} \index{extremal types theorem!outline proof}


The proof of the Extremal Types Theorem is not a very complicated proof, but it is rather lengthy and technical \citep{leliro83}. Since it is not central to the focus of this book, only a sketch will be given here to illustrate the main ingredients. The concept of max-stability is needed. It is defined as follows,

\bigskip
A distribution $G$ is called max-stable if, for every $m = 2, 3, \ldots$, there are constants $\alpha_m > 0$ and $\beta_m$ such that
\begin{align}\label{eq:8.45}
G^m (\alpha_m x + \beta_m) = G(x) \,, 
\end{align}  
\bigskip

$G^m$ is the distribution function of $M_m = \max\{Z_1, \ldots, Z_m \}$, where the $Z_i$ are iid random variables with distribution function $G$. Therefore, max-stability is a property satisfied by distributions that are invariant under the operation of taking sample maxima, except for a change of scale and location. The following result brings forward the connection between max-stability and extreme value distributions \citep{leliro83},

\bigskip

A distribution is max-stable if, and only if, it is a GEV distribution.

\bigskip

To check that a GEV distribution is max-stable, is a straightforward exercise in algebra. The converse is much harder. Anyway, this result can now be used to prove the Extremal Types Theorem. Consider first $M_{nk} = \max\{X_1, \ldots, X_{nk}  \} $ of a sample of $nk$ iid random variables $X_i$, for some large value of $n$. This large sample can be divided into $k$ subsamples of $n$ variables in each. Hence, there will be $k$ iid random variables like $M_{n} = \max\{X_1, \ldots, X_{n} \}$. $n$ is chosen large enough to claim that,     
\begin{align}\label{eq:8.46}
\Prob \Big( \frac{M_n - b_n}{a_n} \le x \Big) \approx G(x) \,, 
\end{align}
for suitable constants $a_n$ and $b_n$ and for the limiting distribution $G$. Hence, for any integer $k \ge 2$, since $nk > n$,
\begin{align}\label{eq:8.47}
\Prob \Big( \frac{M_{nk} - b_{nk}}{a_{nk}} \le x \Big) \approx G(x) \,, 
\end{align}
Eq.~\eqref{eq:8.46} leads to $\Prob \big(M_n  \le z \big) \approx G \big( (z-b_n)/a_n \big)$, while Eq.~\eqref{eq:8.47} gives $\Prob \big(M_{nk}  \le z \big) \approx G \big( (z-b_{nk})/a_{nk} \big)$. However, $M_{nk}$ is obviously the maximum of $k$ variables having the same distribution as $M_n$.  But, then, 
\begin{align}\label{eq:8.48}
\Prob \big( M_{nk}  \le z \big) = \Big[ \Prob \big( M_n  \le z \big) \Big]^k \,, 
\end{align}

From this, it is deduced that (in the limit),
\begin{align}\label{eq:8.49}
G \Big( \frac{z - b_{nk}}{a_{nk}} \Big) = G^k \Big( \frac{z - b_{n}}{a_{n}}  \Big) \,, 
\end{align}
From this it follows that $G$ and $G^k$ are identical apart from location and scale parameters. Hence, G is max-stable, and by the result above, it is a member of the GEV family of distributions. 



\section{Domains of Attraction for the Extreme Value Distributions}\label{sec:daevd} \index{extreme value distribution!domains of attraction}

In practice, the exact statistical distribution of the data being analyzed, is rarely known. However, in many cases there may be rather strong evidence as to what type of distribution to expect. For instance, average wind speeds over periods of 10 minutes in northern Europe have been found to follow a Weibull type distribution. Then it would be useful to know what kind of extreme value distribution to expect for such data. The answer to such questions is the subject of the theory of domains of attraction for extreme value distributions.  It is beyond the scope of our treatment of this topic here to go into much details, but some useful results seem worthwhile presenting. A more thorough discussion is given by \cite{gned43} and \cite{leliro83}. 

 A time series $X_1, X_2, \ldots$ of iid random variables with distribution function $F$, and with a density function $f$, is considered. $x_F$ is defined to be the right endpoint of $F$ by $x_F = \sup \{x ; F(x) < 1\} \quad (x_F \le \infty)$. Then the following sufficient conditions due to von Mises apply:

\bigskip

 Suppose that $F$ is absolutely continuous with density $f$. Then sufficient conditions for $F$ to belong to each of the three possible domains of attractions are:
  
 \bigskip
 {\it Type I}:   $f$ has a negative derivative $f'$ for all $x$ in some interval $(x_0, x_F)$, $(x_F \le \infty)$, and 
 \begin{equation}\label{eq:vm1}
 \lim_{x \nearrow x_F} \frac{f'(x) (1-F(x))}{f^2(x)} = -1. \nonumber
 \end{equation}
 
\bigskip
 {\it Type II}: $f(x) > 0$ for $x \ge x_0$ finite, and for some constant $\alpha > 0$,
 \begin{equation}\label{eq:vm2}
 \lim_{x \to \infty} \frac{x f(x)}{1-F(x)} = \alpha. \nonumber
 \end{equation}

\bigskip
{\it Type III}: $f(x) > 0$ for all $x$ in some finite interval $(x_0, x_F)$, $f(x) = 0$ for $x > x_F$, and for some constant $\alpha > 0$,
\begin{equation}\label{eq:vm3}
\lim_{x \nearrow x_F} \frac{(x_F - x) f(x)}{1-F(x)} = \alpha. \nonumber
\end{equation}

\bigskip
Using these results, it is straightforward to verify that the following list of distributions belong to the domain of attraction of the Type I (Gumbel) case, just to mention a few well-known cases: Normal, lognormal, exponential, Weibull, gamma, and, of course, the Gumbel distribution itself.  

Distributions belonging to the domain of attraction of Type II are e.g. the Pareto, the Generalized Pareto for positive shape parameter, and the Type II extreme value distribution itself. For Type III may be mentioned e.g. the uniform distributions, distributions truncated on the upper side (provided a smooth density function), and the Type III extreme value distribution itself.



\section{Parameter Estimation for the GEV Distributions}\label{sec:pegev} \index{GEV distribution!parameter estimation}\index{parameter estimation!GEV distribution}



The practical application of the block maxima method involves the need to decide on how to divide the observed data into blocks. Obviously, there will be two conflicting issues that have to be dealt with. The desire to have large blocks so that the distribution of the block maxima will approximate a GEV distribution, may easily lead to a sample of few block maxima. Statistical inference on small samples may entail large uncertainties. On the other hand, increasing the sample of block maxima by choosing smaller blocks, may violate the asymptotic approximation by assuming a GEV distribution for the block maxima. These issues may be further complicated by the issues of independence and stationarity, which were discussed in Section~\ref{sec:bmm}. While establishing general rules for the choice of block size relative to the amount of data available is hardly feasible, for some practical cases the accumulated experience has lead to what may be called a consensus. For example, in wind engineering, the choice of one year as a block size has become very close to a standard procedure. An important consideration for this choice is that the data may then reasonably be assumed to belong to the same population since seasonal effects have effectively been removed.     

When the sample of block maxima has been determined, the next step would be to estimate the parameters of the GEV model, or one of the three types, if that can be ascertained a priori. In this book the focus is on two rather popular estimation methods, which is the method of moments (primarily for the Gumbel model) and the maximum likelihood method. The probability weighted moments method has also been used to some extent. For this method, please cf. \cite{howawo85}. 

To simplify notation, the block maxima are denoted by $Z_1, \ldots, Z_k$, assuming $k$ blocks. These random variables are assumed to be iid with a common GEV distribution the parameters of which are to be estimated from the outcomes of the block maxima, that is, the observed data.   


\subsection{Estimation by the method of moments} \index{GEV distribution!method of moments}
  
The exposition of the method of moments for parameter estimation is limited to the Gumbel model. It is for this case that it seems to be most popular, maybe due to its simplicity in this case. The general Gumbel model has two parameters. Since the first two statistical moments  $m_1 = \EE(Z)$ and $m_2 = \EE(Z^2)$ of a Gumbel distributed variable $Z$ can be expressed in term of the these two parameters, for estimation the following two empirical moments are calculated:
\begin{align}\label{eq:8.50}
\hat{m}_j = (1/k) \sum_{i=1}^{k} z_i^j \,, j=1,2,
\end{align}  
where $z_1, \ldots, z_k$ are the observed data. 

Assuming that $Z$ has the general Gumbel distribution $G(z) = \exp \{-\exp[-(z-\mu)/\sigma]\}$, then $m_1 = \EE(Z) = \mu + 0.5772 \sigma$ and $m_2 = \EE(Z^2) = m_1^2 + \pi^2 \sigma^2 / 6$, cf. Section~\ref{sec:alevd}. Denote by $\mu_k$ and $\sigma_k$ the estimated values of the parameters based on the $k$ observations of block maxima. It is then obtained that,
\begin{equation}\label{eq:51}
\mu_k = \hat{m}_1 - 0.5772 \sigma_k
\end{equation}
and
\begin{equation}\label{eq:52}
\sigma_k = (\sqrt{6}/\pi) \sqrt{\hat{m}_2 - \hat{m}_1^2}
\end{equation}





\subsection{Maximum likelihood estimation} \index{GEV distribution!maximum likelihood estimation}

The maximum likelihood (ML) method is very popular and has widespread use in almost every branch of statistics. It turns out that the application of the ML methods for estimation on GEV models requires some caution. Fortunately, it seems that for the applications relevant for this book, the restrictions that need to be observed are rarely an issue. Specifically, for values of the shape parameter $\gamma > - 0.5$, the ML estimators behave regularly. The only thing to note is that there are some small sample issues related to the use of ML estimators also for GEV models, cf. \cite{codi99}.

Based on the assumption that $Z_1, \ldots, Z_k$ are iid random variables having a common GEV distribution, then the log-likelihood function for the GEV parameters when $\gamma \ne 0$ has the following expression:

\begin{equation}\label{eq:53}
\ell(\mu, \sigma, \gamma) = - k \log \sigma - (1 + 1/\gamma)\sum_{i=1}^k \log \left[1 + \gamma \Big(\frac{z_i - \mu}{\sigma} \Big) \right] - \sum_{i=1}^k \left[1 + \gamma \Big(\frac{z_i - \mu}{\sigma} \Big) \right]^{- 1/ \gamma},
\end{equation}
provided that 
\begin{equation}\label{eq:54}
1 + \gamma \Big(\frac{z_i - \mu}{\sigma} \Big) > 0, \quad \mbox{for} \quad i = 1, \ldots, k.
\end{equation}
If the last condition is violated, the likelihood becomes zero, and the log-likelihood therefore $- \infty$. The case $\gamma = 0$ needs to be considered separately, using the Gumbel model. In this case the log-likelihood becomes
\begin{equation}\label{eq:55}
\ell(\mu, \sigma) = - k \log \sigma - \sum_{i=1}^k  \Big(\frac{z_i - \mu}{\sigma} \Big)  - \sum_{i=1}^k \exp \left\{ - \Big(\frac{z_i - \mu}{\sigma} \Big) \right\}.
\end{equation}

To obtain the numerical maximum likelihood estimates from the observed data by using \eqref{eq:53} and \eqref{eq:55}, standard numerical optimization programs may be utilized. If \eqref{eq:53} is used, care must be exercised to avoid numerical problems in cases where the optimization algorithms tend to parameter estimates in the close vicinity of $\gamma = 0$. Then it is strongly advisable to use \eqref{eq:55}.

Confidence intervals on the estimated parameter values can be calculated exploiting that the approximate distribution of the estimators $(\hat{\mu}, \hat{\sigma}, \hat{\gamma})$ is multivariate normal with mean value $(\mu, \sigma, \gamma)$. This is discussed by \cite{cole01}. 




\section{Model Validation} \index{GEV distribution!model checking}

As is well known from basic courses in statistics, the use of probability (or PP) plots and quantile (or QQ) plots may reveal very useful information about the extent of agreement between an assumed or estimated probability distribution and the empirical distribution of the data. These are also highly useful tools for a visual check of fitted GEV models in particular cases. For a thorough discussion of the use of these plots, cf. \cite{begosete04}.   

A probability or PP plot is a direct comparison of the fitted distribution model to the empirical distribution. Assume that the sample of block maxima has been ordered by increasing value: $z_{(1)} \le z_{(2)} \le \ldots \le z_{(k)}$. The empirical distribution function, $\tilde{G}$ say, evaluated at $z_{(i)}$ is given by,
\begin{equation}
\tilde{G}(z_{(i)}) = i/ (k+1). 
\end{equation}
The proposed GEV model distribution is obtained by substituting the parameter estimates into \eqref{eq:8.39},
\begin{equation}\label{eq:gev.x}
\hat{G}(z_{(i)}) = \exp \left\{ - \left[ 1 + \hat{\gamma} \Big( \frac{z_{(i)} -
	\hat{\mu}}{\hat{\sigma}} \Big) \right]^{- 1/ \hat{\gamma}} \right\}, 
\end{equation}
provided $\hat{\gamma} \ne 0$. If $\hat{\gamma} = 0$, the plot is constructed using the Gumbel distribution. If the GEV model is a good approximation, then
\begin{equation}\label{eq:gev.y}
\hat{G}(z_{(i)}) \approx \tilde{G}(z_{(i)}) 
\end{equation}
for each index $i$, so that the PP plot consisting of the points
\begin{equation}\label{eq:gev.z}
\Big(\hat{G}(z_{(i)}), \tilde{G}(z_{(i)}) \Big) \:\: i = 1, \ldots, k,   
\end{equation}
should follow approximately the unit diagonal.  

For the case of extreme value distributions, a quantile or QQ plot is usually considered to be more informative than a PP plot because it shows more clearly the agreement at high values of the observed data, which is of primary concern when fitting extreme value models. Assuming again that $\hat{\gamma} \ne 0$, the QQ plot is traced out by the point graph, 
\begin{equation}
\big( \hat{G}^{-1}(i/ (k+1)) , z_{(i)} \big) , \:\: i = 1, \ldots, k,  
\end{equation}
where 
\begin{equation}
\hat{G}^{-1}(i/ (k+1)) = \hat{\mu} - \frac{\hat{\sigma}}{\hat{\gamma}} \left[ 1 - \left \{ - \log \big(i/ (k+1) \big)  \right \}^{- \hat{\gamma}} \right] .  
\end{equation}
This graph should also approximately follow a straight line. These procedures are discussed at greater length in Chapter~\ref{chap:mc}. 





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

The bootstrapping method is a statistical technique of fairly
recent origin that can be used for estimating confidence
intervals on quantities derived from a statistical distribution on the basis of a limited
sample generated by that same distribution \citep{efti93,dahi97}. It is based on resampling from a
distribution determined by the available sample of data. Despite
the fact that the name of the method alludes to lifting oneself up
by the bootstraps (Baron von Munchausen), the method appears to be
reasonably effective for the specific purpose of estimating
confidence bands. For convenience, a brief discussion of some
basic features of the bootstrapping method is provided here.

Assume that $\mathbf{z} = (z_1, z_2, \ldots, z_n)$ is a
sample or vector consisting of $n$ independent observations of a
random variable $Z$, and that this is the only empirical information available about $Z$. Confidence intervals for a statistical quantity require the estimation of quantiles from the distribution of a relevant estimator. There are in principle two available options for obtaining bootstrap estimates of such quantiles. One is the nonparametric
approach,\index{bootstrapping!nonparametric} where a purely empirical distribution function is
established for $Z$ on the basis of the observed data by
allocating a probability of $1/n$ to each of the observed data
points. The other is the parametric bootstrap,\index{bootstrapping!parametric} which is obtained
by assuming that $Z$ has a specified distribution function $F_Z(z;\theta)= \Prob(Z \le z)$, where $\theta$
denotes a vector of unknown parameters, which determine the distribution.
These parameters are then estimated from the observed data
$\mathbf{z}$, giving $\hat{\theta}$, and
$F_Z(z;\hat{\theta})$ is adopted as the distribution of $Z$.

In this section on the block maxima method using GEV models, only the parametric bootstrap
is used. The goal is to estimate some statistical quantity $V$, e.g. a high quantile like $100 (1 - \alpha)\%$ ($0 < \alpha << 1$), given by the unknown distribution. Let $\hat{V}$ denote the estimate of $V$ obtained from the fitted model distribution $F_Z(z;\hat{\theta})$, which is a GEV distribution. The parametric bootstrapping\index{bootstrapping!parametric} technique for estimating confidence intervals on $V$ is based on resampling from the GEV model obtained.

This is done as follows: Let $Z^*$ denote the random variable with distribution function $F_Z(z;\hat{\theta})$. $\ell$ bootstrap samples $\mathbf{z}^*_j$,
$j=1, \ldots, \ell$, with $n$ independent observations of $Z^*$ in
each sample are now generated. Each sample
$\mathbf{z}^*_j$ is used to fit a new GEV model from which an estimate $V^*_j$ of
$V$ is obtained.

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}\label{eq:confint}
(\,\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 method, usually 20-30
bootstrap samples are sufficient. However, to avoid making the assumption
that the bootstrap estimates $V^*_j$ are generated by a normal distribution, which is the basis for Eq.~\eqref{eq:confint}, 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$). Such estimates may be further improved as described by
\cite{dahi97}. However, such details will not be discussed here.





\section{The Asymptotic Extreme Value Distributions for Dependent Sequences} \index{extreme value theory!classical}
\index{extreme value distribution!dependent data}\label{sec:dss}

The assumed sequence of iid random variables underlying the classical approach to extreme value distributions is obviously not a very practical model for many physical phenomena where dependence effects are obvious. Fortunately, it turns out that the extremal types theorem still applies provided some conditions are satisfied. These conditions relate to the long range dependence structure of the sequence of random variables. The typical model adopted is one of a stationary time series $X_1, X_2, \ldots$. Stationarity means that the joint probability law of a group of random variables from the sequence is invariant with respect to time shifts. That is, e.g. ${X_1, X_2}$ has the same joint distribution as ${X_{51}, X_{52}}$.


The condition that has to be satisfied by the long range dependence of the stationary time series to allow for an extremal types theorem, can be formulated as follows:

\bigskip
A stationary time series $X_1, X_2, \ldots$ is said to satisfy the $D(u_n)$ condition if, for all $i_1< \ldots < i_p < j_1 < \ldots < j_q$ with $j_1 - i_p > l$,

\begin{align}\label{eq:8.61}
&\left| \Prob \Big( X_{i_1} \le u_n, \ldots, X_{i_p} \le u_n, X_{j_1} \le u_n, \ldots,  X_{j_q} \le u_n \Big) - \right. \nonumber \\
&\left. \Prob \Big( X_{i_1} \le u_n, \ldots, X_{i_p} \le u_n \Big) \Prob \Big(  X_{j_1} \le u_n, \ldots,  X_{j_q} \le u_n \Big) \right|    \le \alpha(n, l), 
\end{align}
where $\alpha(n, l_n) \to 0$ for some sequence $l_n$ such that $l_n/n \to 0$ as $n \to \infty$.
\bigskip

A scrutiny of this condition conveys the understanding that it is required that block maxima tend to become independent random variables if the blocks are sufficiently far apart. Provided that this condition is satisfied, the following theorem applies \citep{leliro83}:

\bigskip
Let $X_1, X_2, \ldots$ be a stationary time series, and define $M_n = \max\{X_1, \ldots,
X_n\}$. If there exist sequences of constants $a_n > 0$ and $b_n$ such
that
\begin{align}\label{eq:8.62}
\Prob \Big( \frac{M_n - b_n}{a_n} \le x \Big) \to G(x) \,, \:\: n
\to \infty,
\end{align}
where $G(x)$ is a nondegenerate distribution function, and the $D(u_n)$ condition is satisfied with $u_n = a_n z + b_n$ for every real $z$, and $G(z) > 0$, then $G(x)$ belongs to the class of generalized extreme value distributions.
\bigskip


The practical significance of this theorem is, in fact, substantial since very few time series met in practice would consist of independent data. Hence, without it, the application of GEV distributions to practical problems, would be very hard to justify. Its application, of course, presupposes that the $D(u_n)$ condition is satisfied, which may seem like a highly nontrivial criterion to check for a given stationary time series. Fortunately, it turns out that for most applications it may be routinely assumed to be satisfied. E.g. a stationary Gaussian time series $X_1, X_2, \ldots$ with an autocovariance function $\rho_n = \EE \left[(X_i - \mu)(X_{i+n} - \mu)\right]$, where $\mu = \EE{X_i}$, will satisfy $D(u_n)$ if $\rho_n \log n \to 0$ when $n \to \infty$. Time series of measured data in engineering applications where $\rho_n$ decays slower that $1/ \log n$ is actually very hard to imagine. 
However, it should be clearly understood that the convergence to the appropriate asymptotic limit may depend very much on the dependence structure of the time series. For instance, if there is strong dependence between consecutive data points of the time series, the convergence to the asymptotic limit may be very slow. 

From the original time series $X_1, X_2, \ldots$, a time series of independent variables $\tilde{X}_1, \tilde{X}_2, \ldots$ may be constructed. Let $\tilde{M}_n = \max\{\tilde{X}_1, \ldots,
\tilde{X}_n\}$. Then the following theorem has been proved by \cite{lead83},

\bigskip
 If there exist sequences of constants $a_n > 0$ and $b_n$, and a nondegenerate distribution function $\tilde{G}(x)$ such that,
\begin{align}\label{eq:8.63}
\Prob \Big( \frac{\tilde{M}_n - b_n}{a_n} \le x \Big) \to \tilde{G}(x) \,, \:\: n
\to \infty,
\end{align}
if the $D(u_n)$ condition is satisfied with $u_n = a_n z + b_n$ for every real $z$, and $G(z) > 0$, and if $\Prob \big( (M_n - b_n) /{a_n} \le x \big)$ converges for some $x$, then
\begin{align}\label{eq:8.64}
\Prob \Big( \frac{M_n - b_n}{a_n} \le x \Big) \to G(x) = \big(\tilde{G}(x) \big)^\theta \,, \:\: n
\to \infty,
\end{align}
for some constant $\theta \in [0, 1]$.
\bigskip

The constant $\theta$ is called the \textit{extremal index}, and unless it is equal to one, the limiting distributions for the independent and the original stationary sequences, are not the same. If $\theta > 0$, then $G(x)$ is an extreme value distribution, but with different parameters than $\tilde{G}(x)$. If $(\mu, \sigma, \gamma)$ are the parameters of $G(x)$, and $(\tilde{\mu}, \tilde{\sigma}, \tilde{\gamma})$ are the parameters of $\tilde{G}(x)$, then their relationship is 
\begin{align}\label{eq:8.65}
\gamma = \tilde{\gamma}, \quad \mu = \tilde{\mu} - \frac{\tilde{\sigma}}{\gamma}\, (1 - \theta^\gamma), \quad \sigma = \tilde{\sigma} \, \theta^\gamma,
\end{align}
or, if $\gamma = 0$, taking limits, it is obtained that $\mu = \tilde{\mu} + \sigma \log \theta$ and $\sigma = \tilde{\sigma}$. Note that the shape parameter $\gamma$ remains the same.


The results cited previously are derived under the assumption of
a sequence of random variables. For many applications in this
book, a sequence of random variables $X_1, X_2,
\ldots$ is initially not considered, but rather a stochastic process in continuous time
$X(t)$. To fit into the framework discussed in this section, one
could envisage sampling the considered process at discrete time
points and obtaining a sequence $X_j = X(t_j)$, $j = 1, 2,
\ldots$. In practice, this is often done by extracting local peak values from an observed realization of the process. As has been seen, if the obtained time series is stationary, then, under suitable conditions, the extremal types theorem still holds true. Unfortunately, no satisfactory general
theory of extremes is available for the case of continuous time
stochastic processes, but some results have been proven
\citep{leliro83}. Of course, in order for a realization of a stochastic process to be stored in computer memory, it has to be sampled. In that sense, what is available for further analysis is, in fact, a time series of outcomes of random variables. Hence, the sampling frequency relative to the characteristics of the stochastic process determines to what extent one may consider the stored time series to be a good replica of the realization of the stochastic process, which may then be used for further processing, like extraction of peak values. 





