\chapter{Classical models for extremes} \section{Background and theoretical motivation} \noindent Suppose that $X_{1}, X_{2}, \ldots, X_{n}$ is a sequence of independent and identically distributed (IID) random variables with common distribution function $F$. One way of characterising extremes is by considering the distribution of the maximum order statistic \begin{eqnarray} M_{n}&=& \textrm{max}\left\{X_{1}, X_{2}, \ldots, X_{n}\right\}. \end{eqnarray} Think back to the example in Section 1.4 concerning sea surges at Wassaw Island. Sea surge measurements were taken every hour; for each year (1955--2004) you were presented with the \textit{annual maximum} sea surge. Thus, here $n=365\times 24 = 8760$ (for non--leap years, anyway), and we might use the notation: \[ M_{8760,i}, \qquad i=1, \ldots, 50, \] to denote, generally, the 50 annual maxima given in Table 1.1 (of course, for leap years we would write $M_{8784,i}$). The assumption of IID might be reasonable here. For example, each annual maximum is likely to occur during the hurricane season (which is usually at its peak in September/October), and so it seems likely that maximum hourly observations from one year to the next will be far enough apart to be independent; issues of non--stationarity, however, might arise in long--range datasets owing to the effects of climate change, for example, and we will come back to this in Chapter 4. \noindent Assuming, for now, that our maxima are IID, how can we obtain the distribution of $M_{n}$? \framebox[17cm]{ \parbox[t]{16cm}{ \huge \ding{46} \normalsize \vspace{4cm}}} \subsection{A brief history of Extreme Value Theory} \noindent One of the earliest books on the statistics of extreme values is \textcolor{gray}{\bf E.J. Gumbel} (1958, see Figure 2.1). Gumbel traces the origins back to 1709, when \textcolor{gray}{\bf N. Bernoulli} considers the problem of estimating he age of the longest survivor in a group of people. Research into extreme values as a subject in it's own right began much later, between 1920 and 1940, when work by \textcolor{gray}{\bf E.L. Dodd}, \textcolor{gray}{\bf M. Fr\'{e}chet}, \textcolor{gray}{\bf E.J. Gumbel}, \textcolor{gray}{\bf R. von Mises} and \textcolor{gray}{\bf L.H.C. Tippett} investigated the asymptotic distribution of the \textit{largest order statistic}. This led to the main theoretical result: the \textit{Extremal Types Theorem} (see Section 2.1.2), which was developed in stages by Fisher, Tippett and von Mises, and eventually proved in general by \textcolor{gray}{\bf B. Gnedenko} in 1943. Until 1950, development was largely theoretical. In 1958, Gumbel started applying theory to problems in engineering. In the 1970s, \textcolor{gray}{\bf L. de Haan} and \textcolor{gray}{\bf J. Pickands} generalised the theoretical results, giving a better basis for statistical models. Since the 1980s, methods for the application of Extreme Value Theory have become much more widespread. Current researchers who have played a significant role in developing applications and methodology include \textcolor{gray}{\bf Richard Smith} (University of North Carolina), \textcolor{gray}{\bf Anthony Davison} (Ecole Polytechnique Federale de Lausanne) \textcolor{gray}{\bf Jonathan Tawn} (Lancaster University) and \textcolor{gray}{\bf Stuart Coles} (University of Padua and \textsc{SmartOdds}). There are still gaps between the theory and the models, and also between the models and common practice in applications -- this is where our work fits in (\textcolor{gray}{\bf Fawcett} and \textcolor{gray}{\bf Walshaw}, see Figure 2.1), and this is the area of research we will focus on in the latter part of this course. \begin{figure}[!h] \centering \includegraphics[height=7.1cm]{gumbel.ps} \hspace{0.7cm} \includegraphics[height=7.1cm]{gumbelbook.ps} \hspace{0.7cm} \includegraphics[height=7.1cm]{LFbook.ps} \caption{Emil Julius Gumbel (left), and his book ``Statistics of Extremes'' (centre); the latest book to hit the shops in August 2016 (right, pre--order on Amazon now!!).} \end{figure} \newpage \subsection{The Extremal Types Theorem} \noindent The obvious question now, is, what possible distributions might be considered candidates for $M_{n}$ as $n \to \infty$? It might also be convenient to find out if we can formulate this set of candidate distributions into a single class -- say $G$ -- which is independent of $F$; if so, we can estimate the distribution of $M_{n}$ using $G$, without any reference to $F$. \noindent Clearly, the limiting distribution of $M_{n}$ is \textit{degenerate}, that is, the distribution converges to a single point on the real line with probability 1, this single point being the upper endpoint of $F$ (which itself will be $\infty$ in some applications). This is analogous to the sample mean $\bar{X}$ converging to the population mean $\mu$ with certainty in the Central Limit Theorem; here, the degenerate limit is prevented by allowing a linear rescaling, so that \[ \frac{\bar{X}-b_{n}}{a_{n}}\xrightarrow{\hspace{0.25cm} D \hspace{0.25cm} } N(0,1) \] where $b_{n} = \mu$ and $a_{n}=\sigma/\sqrt{n}$, where $\sigma$ and $n$ are the population standard deviation and sample size, respectively. Can we apply a similar linear rescaling to $M_{n}$ to avoid convergence of the distribution to a single point? The answer, of course, is ``yes'', and is provided by the main result in classical extreme value theory -- the \textit{Extremal Types Theorem} -- a result for the maximum $M_{n}$ which is analogous to the Central Limit Theorem for the mean $\mu$. \vspace{0.8cm} \noindent {\large \textbf{Theorem} \textit{(The Extremal Types Theorem)}} If there exist sequences of constants $a_{n}>0$ and $b_{n}$ such that, as $n \to \infty$, \begin{eqnarray} \textrm{Pr}\left\{(M_{n}-b_{n})/a_{n}\leq x\right\} \to G(x) \end{eqnarray} for some non--degenerate distribution $G$, then $G$ is of the same type as one of the following distributions: \begin{eqnarray} \mbox{\emph{I}}: G(x) &=& \textrm{exp}\left\{-\textrm{exp}(-x)\right\} \qquad -\infty0, \alpha>0; \end{array} \right.\\ \mbox{\emph{III}}:G(x) &=& \left\{ \begin{array}{ll} \textrm{exp}\left\{-(-x)^{\alpha}\right\} &x<0, \alpha>0 \\ 1 &x \geq 0. \end{array} \right. \end{eqnarray} The three types of distribution in the Extremal Types Theorem ($\mbox{\emph{I}}$, $\mbox{\emph{II}}$ and $\mbox{\emph{III}}$) have become known as the Gumbel, Fr\'{e}chet and Weibull types (respectively), and are known collectively as the \textit{extreme value distributions}. For both the Gumbel and Fr\'{e}chet distributions the limiting distribution $G$ is unbounded; that is, the upper--endpoint tends to $\infty$. The Weibull distribution has a finite upper bound. It should be noted that the Extremal Types Theorem does not \textit{ensure} the existence of a non--degenerate limit for $M_{n}$; nor does it specify \textit{which} of types $\mbox{\emph{I}}$, $\mbox{\emph{II}}$ or $\mbox{\emph{III}}$ is applicable \textit{if} a limit distribution \textit{does} exist (i.e. in which \textit{domain of attraction} the distribution of $G$ lies). However, when such a distribution does exist, we find that, by analogy with the Central Limit Theorem, the limiting distribution of sample maxima follows one of the distributions given by the Extremal Types Theorem, \textit{no matter what the parent distribution $F$.} \noindent So we know that \[ \frac{M_{n}-b_{n}}{a_{n}} \xrightarrow{\hspace{0.25cm}D\hspace{0.25cm}} G, \] where -- if it exists -- $G$ is given by one of the extreme value distributions (either the Gumbel, Fr\'{e}chet or Weibull distributions). But how do we know which one of these distributions to use? \subsection{The Generalised Extreme Value (GEV) distribution} \noindent In practice, working with, and having to choose between, three distributions is inconvenient. However, there exists a parameterisation which encompasses all three types of extreme value distribution. Von Mises (1954) and Jenkinson (1955) independently derived the \textit{generalised extreme value distribution} (GEV), often denoted $\mathcal{G}(\mu, \sigma, \xi)$, whose cumulative distribution function is given by \begin{eqnarray} G(x;\mu, \sigma, \xi) &=& \textrm{exp}\left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}\right\}, \end{eqnarray} where $a_{+} = \textrm{max}(0,a)$. The situation where $\xi=0$ is not defined in (2.6), but is taken as the limit as $\xi \to 0$, given by \begin{eqnarray} G(x;\mu, \sigma)&=&\textrm{exp}\left\{-\textrm{exp}\left(\frac{x-\mu}{\sigma}\right)\right\}. \end{eqnarray} The parameters $\mu$ ($-\infty<\mu<\infty$), $\sigma$ ($>0$) and $\xi$ ($-\infty<\xi<\infty$) are location, scale and shape parameters respectively. The value of the shape parameter $\xi$ differentiates between the three types of extreme value distribution in Theorem 2.1: $\xi=0$, leading to (2.7), corresponds to the Gumbel distribution (type $\mbox{\emph{I}}$); $\xi>0$ and $\xi<0$ correspond to the Fr\'{e}chet (type $\mbox{\emph{II}}$) and Weibull (type $\mbox{\emph{III}}$) distributions respectively. In practice, when we estimate the shape parameter $\xi$, the standard error for $\xi$ accounts for our uncertainty in choosing between the three models given by (2.3)---(2.5). But what about the constants $a_{n}$ and $b_{n}$? We know that \[ \frac{M_{n}-b_{n}}{a_{n}} \xrightarrow{\hspace{0.25cm} D \hspace{0.25cm}} \mathcal{G}(\mu,\sigma,\xi), \qquad \textrm{ as } n\to \infty. \] After some algebra, it turns out that \[ M_{n} \xrightarrow{\hspace{0.25cm} D \hspace{0.25cm}} \mathcal{G}(\mu^{*}, \sigma^{*}, \xi), \qquad \textrm{ as }n \to \infty, \] with $a_{n}$ and $b_{n}$ being absorbed into $\mu^{*}$ and $\sigma^{*}$. Since the GEV parameters need to be estimated anyway, in practice we just ignore the normalisation constants and fit the GEV directly to our set of maxima $M_{n,i}$. However, before we consider applications of the GEV to real data, let us first consider some theoretical examples which demonstrate that, with careful choices of $a_{n}$ and $b_{n}$, one of the three extreme value distributions is always achieved when the parent distribution $F$ is known. \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \subsection*{Example 2.1} \noindent Suppose $X_{1}, X_{2}, \ldots, X_{n}$ is a sequence of independent $\textit{Exp}(1)$ variables, that is \[ F(x) = 1-e^{-x}, \qquad x>0. \] By letting $a_{n}=1$ and $b_{n}=\textrm{log }n$, show that the limit distribution of $(M_{n}-b_{n})/a_{n}$ is of extreme value type, and identify the distribution. \vspace{0.6cm} {\large \textbf{Solution}} \huge \ding{46} \normalsize \vspace{15cm} }} \vspace{0.5cm} \framebox[17cm]{ \parbox[t]{16cm}{ \subsection*{Example 2.2} Suppose $X_{1}, X_{2}, \ldots, X_{n}$ is a sequence of independent $\textit{Fr\'{e}chet}(1)$ variables, that is \[ F(x) = e^{-1/x}, \qquad x>0. \] By letting $a_{n}=n$ and $b_{n}=0$, show that the limit distribution of $(M_{n}-b_{n})/a_{n}$ is of extreme value type, and identify the distribution. \vspace{0.6cm} {\large \textbf{Solution}} \huge \ding{46} \normalsize \vspace{15cm} }} \newpage \subsection{Typical application} \subsubsection{1. Data pre--processing} Unless you already have a filtered set of block/annual maxima, the first step in any application of the GEV is to obtain your $M_{n,i}$. This requires choosing a \textit{block length} $n$ (usually the calendar year), and then discarding all but the maximum value in each block. Be careful -- if your block length is too short, the limiting arguments supporting the GEV might not be valid, and your maxima may be too close together to assume independence; if $n$ is too large, you won't have enough data to work with! \subsubsection{2. Estimation of the GEV parameters} A typical application of the generalised extreme value distribution is to fit (2.6) to a series of block (often annual) maxima. Numerous ways of fitting this to a set of maxima have been considered, though numerical maximum likelihood estimation has become generally accepted as the most robust procedure. For example, differentiation of (2.6) gives the probability density function of the GEV; this can be found to be \begin{eqnarray} g(x; \mu, \sigma, \xi) &=& \frac{1}{\sigma}\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]_{+}^{-(1/\xi+1)}\textrm{exp}\left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}\right\}. \end{eqnarray} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} Use the GEV probability density function in (2.8) to form the likelihood function $L(\mu,\sigma,\xi;\boldsymbol{x})$. Also, obtain the GEV log--likelihood function $\ell(\mu,\sigma,\xi;\boldsymbol{x})$ and the corresponding (log) likelihood equations. \huge \ding{46} \normalsize \vspace{9.25cm} }} \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} How would you use the (log) likelihood equations in order to obtain maximum likelihood estimates of $\mu$, $\sigma$ and $\xi$? Why can't we obtain closed form solutions for $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\xi}$? How can we get around this? \huge \ding{46} \normalsize \vspace{5cm} }} The usual techniques can be used to obtain standard errors for the estimated GEV parameters, i.e. via inversion of the information matrix: these can be used to construct confidence intervals and perform hypothesis tests. \subsubsection{3. Checking model adequacy} As with all statistical models, there are various goodness--of--fit properties that should be considered to check the overall adequacy of the fitted GEV. These include probability plots, quantile--quantile ($Q$-$Q$ plots) and simply plotting a histogram of the data against the fitted density. Again, these will be reviewed shortly via a real--life data demonstration in \textbf{\sf R}. \subsubsection{4. Return level estimation} \noindent As discussed in Section 1.3, interest usually lies not in estimates of the GEV parameters themselves, but in how we can use the fitted model to estimate other quantities, such as \begin{itemize} \item The height of a sea wall to protect against the once in a hundred year sea--surge; \item The ``fifty year wind speed'' to provide new structures enough protection against wind damage. \end{itemize} Such quantities, in extreme value terminology, are usually referred to as \textit{return levels}. Provided we have faith in our fitted model being suitable beyond the range of our observed data, we can estimate the $r$--year return level $z_{r}$ for any period by setting the GEV distribution function (Equation 2.6) equal to $1-1/r$ and solving for $x=\hat{z}_{r}$, provided we have annual maxima. \noindent For example, suppose we fit the GEV to the set of annual maxima given in Table 1.1 and obtain estimates of the location, scale and shape as $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\xi}$ (respectively). Suppose further that the authorities require an estimate of $z_{100}$, the sea surge we might expect to be exceeded once in a hundred years (this could be used to help estimate the corresponding height of a new sea wall). Then we can write down the following probability statement: \begin{eqnarray*} \textrm{Pr}(\textrm{annual maximum}>z_{100}) &=& \frac{1}{100}, \end{eqnarray*} i.e. \begin{eqnarray} 1 - \textrm{Pr}(\textrm{annual maximum}\leq z_{100}) &=&\frac{1}{100}. \end{eqnarray} Now the left--hand--side of (2.9), in terms of our fitted GEV, is \[ 1-G(\hat{z}_{100};\hat{\mu},\hat{\sigma},\hat{\xi}), \] giving \begin{eqnarray*} 1-\textrm{exp}\left\{-\left[1+\hat{\xi}\left(\frac{\hat{z}_{100}-\hat{\mu}}{\hat{\sigma}}\right)\right]^{-1/\hat{\xi}}\right\} &=& 0.01 \qquad \textrm{i.e.} \\ & & \\ \textrm{exp}\left\{ -\left[ 1+\hat{\xi}\left(\frac{\hat{z}_{100}-\hat{\mu}}{\hat{\sigma}}\right) \right]^{-1/\hat{\xi}} \right\} &=&0.99. \end{eqnarray*} Solving for $\hat{z}_{100}$ gives an estimate of the 100--year return level as \[ \hat{z}_{100} = \hat{\mu}+\frac{\hat{\sigma}}{\hat{\xi}}\left[\left(-\textrm{log}(0.99)\right)^{-\hat{\xi}}-1\right]; \] more generally, estimates of the $r$--year return level $z_{r}$ are given by \begin{eqnarray} \hat{z}_{r} &=&\hat{\mu}+\frac{\hat{\sigma}}{\hat{\xi}}\left[\left(-\textrm{log}\left(1-r^{-1}\right)\right)^{-\hat{\xi}}-1\right]. \end{eqnarray} Care needs to be taken when $\xi=0$; in this case, as with forming the likelihood in the first place, we should work with the limiting form (as $\xi \to 0$) of the distribution function (i.e. Equation (2.7)). Of course, a point estimate of the return level alone is not enough; any estimate requires a corresponding standard error. As the return level is a function of the GEV parameters, which we have standard errors for, we can use the delta method to obtain standard errors for return level estimates. Further, confidence intervals formed in the usual way are often extremely misleading for return level estimates, owing to the (often severe) asymmetry encountered in the surface of the likelihood for these quantities; rather, the method of \textit{profile likelihood} is preferred to construct confidence intervals here (see later). \newpage \section{Worked example: sea surge at Wassaw Island} \noindent Recall the example from Section 1.4. The southeastern coast of the USA is an area often hit by hurricanes. The historic city of Savannah, Georgia, has suffered direct hits from 22 Hurricanes since 1871, and can expect to `brushed' or directly hit once every other year. Predictions show that the city is three years overdue it's next direct hit (keep an eye on the news this summer). The figures in Table 2.1 below show the annual maximum sea surges observed at Wassaw Island, a barrier island just off the coast of Savannah, between 1955 and 2004 (inclusive). \begin{table}[!h] \centering \begin{tabular}{|c c c c c c c c c c|}\hline 8.5 &8.9 &9.1 &8.9 &8.4 &9.7 &9.1 &9.6 &8.7 &9.3 \\ 9.6 &9.3 &8.7 &9.0 &8.8 &8.9 &8.9 &12.2 &7.8 &7.7 \\ 8.3 &8.1 &7.3 &6.8 &6.7 &7.3 &7.6 &8.2 &8.6 &9.8 \\ 9.5 &7.4 &7.3 &10.2 &10.3 &10.4 &8.8 &9.7 &10.0 &10.8 \\ 11.1 &12.7 &11.5 &11.8 &12.6 &13.0 &10.5 &10.5 &10.0 &9.4 \\ \hline \end{tabular} \caption{Annual maximum sea--surges at Wassaw Island, 1955--2004 (inclusive). Units are feet.} \end{table} Recall that, in Chapter 1, we tried to estimated quantities beyond the range of our observed data. For example, using simple empirical arguments, \[ \textrm{Pr}(\textrm{Sea surge exceeds 14 feet}) = 0, \] simply because, over the time-frame we have data for, we have not observed such an extreme event. Does this really mean this event is \textit{impossible}? Probably not, it's just that we're trying to estimate the probability of an event more extreme than we've ever observed before. We will now proceed with a detailed extreme value analysis of these annual maxima. \subsection{Data pre--processing} \noindent In order to fit the GEV distribution, we need to pre--process our data so we have a filtered set of extremes -- in particular, block, or annual, maxima. In this example no pre--processing is necessary, as we are given a set of annual maxima. In the first computer practical session, we will consider how to pre--process a dataset to obtain the set of block maxima. Figure 2.2 shows a time series plot and histogram of the 50 annual maxima. Notice there is no obvious trend in our dataset, although there could be an issue of dependence, with neighbouring years having similar values of maximum sea surge. These plots have been produced using \textsf{R}, and, in your own time, you should make sure you can reproduce them. The data are available to download/copy--and--paste from the MAS8304 webpage, in the file ``\texttt{wassaw.txt}''. To load them into \textsf{R}, save the file to your own workspace and then type: \begin{verbatim} > wassaw=scan('wassaw.txt') \end{verbatim} Alternatively, you can type them directly into an \textsf{R} terminal: \begin{verbatim} > wassaw=c(8.5, 8.9, ..., 9.4) \end{verbatim} The following code sets up a year counter (from 1955 to 2004): \begin{verbatim} > year=seq(1955,2004,1) \end{verbatim} Then, to produce the plots shown in Figure 2.2, you can type: \begin{verbatim} > par(mfrow=c(1,2)) \end{verbatim} Generically, \texttt{par} partitions the plotting window, and \texttt{mfrow=c(a,b)} tells \textsf{R} to split the window into \texttt{a} rows and \texttt{b} columns. Then: \begin{verbatim} > plot(wassaw~year,type='b') > hist(wassaw) \end{verbatim} Obviously, these are the basic commands; further arguments can be defined to choose different titles/axis labels etc. \begin{figure}[!h] \centering \psfrag{0}[][]{\sf \tiny 0} \psfrag{2}[][]{\sf \tiny 2} \psfrag{4}[][]{\sf \tiny 4} \psfrag{6}[][]{\sf \tiny 6} \psfrag{1960}[][]{\sf \tiny 1960} \psfrag{1970}[][]{\sf \tiny 1970} \psfrag{1980}[][]{\sf \tiny 1980} \psfrag{1990}[][]{\sf \tiny 1990} \psfrag{2000}[][]{\sf \tiny 2000} \psfrag{T}[][]{\sf \scriptsize Time series plot} \psfrag{S}[][]{\sf \tiny Sea surge (feet)} \psfrag{F}[][]{\sf \tiny Frequency} \psfrag{H}[][]{\sf \scriptsize Histogram} \psfrag{7}[][]{\sf \tiny 7} \psfrag{8}[][]{\sf \tiny 8} \psfrag{9}[][]{\sf \tiny 9} \psfrag{10}[][]{\sf \tiny 10} \psfrag{11}[][]{\sf \tiny 11} \psfrag{12}[][]{\sf \tiny 12} \psfrag{13}[][]{\sf \tiny 13} \psfrag{year}[][]{\sf \tiny Year} \includegraphics[width=17cm,height=6cm]{was1.ps} \caption{Time series plot (left) and histogram (right) of annual maximum sea--surges observed at Wassaw Island.} \end{figure} \newpage \subsection{Fitting the GEV} \noindent As discussed in Section 2.1.4, the most commonly--used technique for fitting the GEV is maximum likelihood, which requires maximisation of the GEV log--likelihood function $\ell(\mu,\sigma,\xi;x)$. In particular, for our data we have \begin{eqnarray} \ell(\mu,\sigma,\xi; \boldsymbol{x}) &=& -50\textrm{log}\sigma - (1+1/\xi)\sum_{i=1}^{50}\textrm{log}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+} \nonumber \\ & & \hspace{0.75cm} - \sum_{i=1}^{50}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}, \end{eqnarray} As previously discussed, there is no analytical solution for obtaining the maximum likelihood estimates of the GEV parameters. Instead, we need to maximise the log--likelihood in (2.11) numerically; we do this using \textsf{R}. \subsubsection*{Using \textsf{R} to estimate the GEV parameters} \noindent You should, at some point, try the following code for yourself. Having already stored our data in \texttt{wassaw}, the next step is to define the GEV parameter vector $\theta=(\mu, \sigma, \xi)$, and to provide \textsf{R} with ``starting values'' for the elements of $\theta$. Since $\mu$ and $\sigma$ are the GEV location and scale parameters, it is common to use the mean and standard deviation (respectively) as starting values for these parameters; the shape parameter $\xi$ can be given any starting value, but a common choice is 0.1. Thus, in \textsf{R}: \begin{verbatim} > theta=c(mean(wassaw), sd(wassaw), 0.1) \end{verbatim} Now, the fiddly bit is to define the GEV log--likelihood function. Read through the following code, and corresponding notes, carefully -- you will be expected to do this sort of thing in Assignment 2! \begin{verbatim} > gev.loglik=function(theta){ mu=theta[1] sigma=theta[2] xi=theta[3] m=min((1+(xi*(dataset-mu)/sigma))) #See note 1 if(m<0.00001)return(as.double(1000000)) #See note 2 if(sigma<0.00001)return(as.double(1000000)) #See note 3 if(xi==0){ loglik=-length(dataset)*log(sigma)-sum((dataset-mu)/sigma) + -sum(exp(-((dataset-mu)/sigma)))} #See note 4 else{ loglik=-length(dataset)*log(sigma) + -(1/xi+1)*sum(log(1+(xi*(dataset-mu)/sigma))) + -sum((1+(xi*(dataset-mu)/sigma))**(-1/xi))} #See note 5 return(-loglik)} #See note 6 \end{verbatim} \noindent \textit{Notes on the above code} \begin{enumerate} \item This line stores the following component in the scalar \texttt{m}: \begin{eqnarray*} \left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+}. \end{eqnarray*} \item Recall that, if the component \texttt{m} is less than zero, we take it to be equal to zero, as $a_{+}=\textrm{max}(0,a)$. Thus, if for a particular combination of $\theta=(\mu,\sigma,\xi)$ we get \texttt{m<0.00001}, we return a huge value for the negative log--likelihood (1e+06) ensuring that this combination will almost surely never be the combination that minimises the negative log--likelihood (and so will almost surely never maximise the log--likelihood -- see note 6). \item This line performs a similar task for \texttt{sigma} as the code referred to in note 2 does for \texttt{m}. $\sigma$ is the scale parameter and so can never be negative; thus, any attempt to minimise the negative log--likelihood using $\texttt{sigma<0.00001}$ returns a huge negative log--likelihood (1e+06), again ensuring that this choice of $\sigma$ can almost surely never minimise the negative log--likelihood (again, see note 6). \item This line defines the GEV log--likelihood when $\xi=0$, taken as the limit when $\xi \to 0$ (see Equation (2.7)). \item This line defines the GEV log--likelihood when $\xi \neq 0$ (see Equation (2.6)). \item The \textsf{R} routine \texttt{nlm} is a minimisation routine that minimises a supplied function (here \texttt{gev.loglik}) using a Newton--Raphson type algorithm. Since the aim is to find estimates of the GEV parameters that \textit{maximise} the log--likelihood, this will be equivalent to finding estimates of the GEV parameters that \textit{minimise} the \textit{negative} log--likelihood, which is why we ask the function \texttt{gev.loglik} to return \texttt{-loglik}. \end{enumerate} We implement the \texttt{nlm} routine in the following way, giving the output shown underneath: \begin{verbatim} > dataset=wassaw > nlm(gev.loglik,theta) $minimum [1] 89.52412 $estimate [1] 8.7112735 1.3114836 -0.1084451 $gradient [1] 3.350727e-06 2.316675e-05 2.145839e-06 $code [1] 1 $iterations [1] 10 \end{verbatim} \newpage \noindent Thus, our estimates of the GEV parameters are \begin{displaymath} \hat{\mu}=8.711 \qquad \hat{\sigma}=1.311 \qquad \hat{\xi}=-0.108 \end{displaymath} The \textit{minimised} \textit{negative} log--likelihood is $89.524$, giving a maximised log--likelihood at the above values for the GEV parameters of $-89.524$. \texttt{gradient} gives the gradient at the estimated minimum of the function \texttt{gev.loglik}, \texttt{iterations} gives the number of iterations performed (in the Newton--Raphson algorithm) before the solution to the minimisation problem was found, and \texttt{code} is an integer which indicates why the optimisation process terminated (\texttt{1} indicates the relative gradient is close to zero and so the current iterate is very close to the solution). \subsubsection*{Using \textsf{R} to estimate standard errors of the GEV parameters} \noindent Using \texttt{nlm} to minimise the negative log--likelihood function is fine, but our estimates of the GEV parameters do not have any associated estimate of variability attached. However, using the argument \texttt{hessian} in the \texttt{nlm} routine returns the matrix of second--order partial derivatives (the \textit{Hessian matrix}) of the function we are trying to minimise. In \textsf{R}, we type: \begin{verbatim} > A<-nlm(gev.loglik,theta,hessian=TRUE) \end{verbatim} \noindent This will give the same output as before, but now the Hessian matrix will also be returned. Note we have also stored the output in \texttt{A}. To look at the output, we type: \begin{verbatim} > A $minimum [1] 89.52412 $estimate [1] 8.7112735 1.3114836 -0.1084451 $gradient [1] 3.350727e-06 2.316675e-05 2.145839e-06 $hessian [,1] [,2] [,3] [1,] 27.375570 -4.837481 18.06140 [2,] -4.837481 56.151740 29.89004 [3,] 18.061397 29.890039 118.18086 $code [1] 1 $iterations [1] 10 \end{verbatim} \noindent Recalling theory from MAS2302, our observed information matrix will be \[ I_{O} = \left(\begin{array}{c c c} -\frac{\partial^{2}\ell}{\partial \mu^{2}} & & \\ & & \\ -\frac{\partial^{2}\ell}{\partial \sigma \partial \mu} &-\frac{\partial^{2}\ell}{\partial \sigma^{2}} & \\ && \\ -\frac{\partial^{2}\ell}{\partial \xi \partial \mu} &-\frac{\partial^{2}\ell}{\partial \xi \partial \sigma} &-\frac{\partial^{2}\ell}{\partial \xi^{2}} \end{array} \right), \] evaluated at $\mu=\hat{\mu}$, $\sigma=\hat{\sigma}$ and $\xi=\hat{\xi}$. Notice that this is just the negative of the hessian of the log--likelihood function; however, the hessian returned by \textsf{R} in the code above is for the \textit{negative} log--likelihood (recall that \texttt{gev.loglik} returns \texttt{-loglik}), and so -- in this example -- the hessian is \textit{exactly} the observed information matrix. Recall also from MAS2302/MAS3320 that the variance--covariance matrix is obtained via inversion of the observed information matrix $I_{O}$. One way of doing this is to use the \textsf{R} command \texttt{solve}. This command solves, for $X$, a system of equations of the form \begin{eqnarray*} \mathbf{a}X = \mathbf{b}, \end{eqnarray*} where both $\mathbf{a}$ and $\mathbf{b}$ can be numeric or complex matrices, $\mathbf{a}$ corresponding to the coefficients of the linear system and $\mathbf{b}$ giving the right--hand--side of this system. In our example, $\mathbf{a}$ is the observed information matrix. If $\mathbf{b}$ is missing in the execution of \texttt{solve}, then the identity matrix is assumed and \textsf{R} will thus return the inverse of $\mathbf{a}$ -- exactly what we want here! In \textsf{R}, this gives: \begin{verbatim} > solve(A$hessian) [,1] [,2] [,3] [1,] 0.043869792 0.008491490 -0.008852201 [2,] 0.008491490 0.022223135 -0.006918367 [3,] -0.008852201 -0.006918367 0.011564254 \end{verbatim} \noindent Thus, the square roots of the diagonal elements of the above matrix will give the expected standard errors for the GEV parameters ($\mu$, $\sigma$ and $\xi$ respectively). We could obtain these standard errors in the following way: \begin{verbatim} > varcovar=solve(A$hessian) > sqrt(diag(varcovar)) [1] 0.2094512 0.1490743 0.1075372 \end{verbatim} \noindent The first line of code stores the matrix shown above in the matrix \texttt{varcovar}; the second line of code finds the square roots of the diagonal elements of \texttt{varcovar}. The resulting expected standard errors are shown underneath. Thus, we now have the following inference for our annual maximum sea surges, in terms of the GEV distribution: \begin{displaymath} \hat{\mu}=8.711 \textrm{ (0.209)}\qquad \hat{\sigma}=1.311 \textrm{ (0.149)} \qquad \hat{\xi}=-0.108 \textrm{ (0.108)} \end{displaymath} \noindent From this, we can construct confidence intervals in the usual way (``Wald'' confidence intervals). For example, standard symmetric 95\% confidence intervals are found as: \begin{displaymath} (8.301, 9.121) \qquad (1.019, 1.603)\qquad \textrm{and} \qquad (-0.320, 0.104) \end{displaymath} for $\mu$, $\sigma$ and $\xi$ (respectively). \subsection{Model adequacy} \noindent We usually fit a statistical model to data to draw conclusions about some aspect of the population from which the data were obtained. The more accurate the fitted model is, the more reliable these conclusions are likely to be. Since inferences are sensitive to the accuracy of the fitted model, it is important that we check that the model fits well. Ideally, we would like to check that our model describes well the variations in the wider population, but this is usually not possible unless there are other sources of data for us to judge the model,against. Thus, it is common to assess the goodness--of--fit of a particular model using the data that were used to estimate it in the first place. \noindent A \textit{probability plot} is a plot of the points \begin{displaymath} \left\{\left(\hat{F}\left(x_{(i)}\right),\frac{i}{n+1}\right): i=1, \ldots, n\right\}, \end{displaymath} where $x_{(i)}, i=1, \ldots, n$ is an ordered sample of independent observations and $\hat{F}$ is a candidate model for the true probability distribution $F$. The quantity $i/n+1$ corresponds to the \textit{empirical distribution function} evaluated at $x_{(i)}$. If $\hat{F}$ is a reasonable model for the true distribution, the points in the probability plot will lie close to the unit diagonal. \noindent A \textit{quantile--quantile plot} ($q$--$q$ plot, or just quantile plot) is a plot of the points \begin{displaymath} \left\{\left(\hat{F}^{-1}\left(\frac{i}{n+1}\right),x_{(i)}\right): i=1, \ldots, n\right\}. \end{displaymath} \noindent The quantity $\hat{F}^{-1}(i/(n+1))$ gives a model--based estimate of the $i/(n+1)$ quantile provided by the candidate distribution $\hat{F}$ whilst $x_{(i)}$ itself provides an empirical estimate of this quantile. Again, a well--fitting model would provide points on this plot lying close to the unit diagonal. \subsubsection*{Producing a probability plot in \textsf{R}} Recall that in the last section, the wassaw dataset was stored in \texttt{dataset}. In \textsf{R}, the code \begin{verbatim} > ordered=sort(dataset) \end{verbatim} then stores the ordered annual maximum sea surges observed at Wassaw Island, $x_{(i)}, i=1, \ldots, 50$, in the vector \texttt{ordered}. For each point in our ordered sample, the following code then stores the empirical distribution function defined by $i/(n+1)$ in the vector \texttt{empirical}: \begin{verbatim} > empirical=vector('numeric',length(ordered)) > for(i in 1:length(empirical)){ empirical[i]=i/(length(dataset)+1)} \end{verbatim} \noindent The function \texttt{GEV.DF} defines the distribution function for the GEV, as provided by equations (2.6) and (2.7): \begin{verbatim} > GEV.DF=function(data,mu,sigma,xi){ if(xi==0){ GEV=exp(-exp(-((data-mu)/sigma)))} else{ GEV=exp(-(1+xi*((data-mu)/sigma))^(-1/xi))} return(GEV)} \end{verbatim} \noindent Then the following code stores a model--based estimate of the distribution function, evaluated at each point in the ordered sample \texttt{ordered}, in the vector \texttt{model} (recall that maximum likelihood estimates of the GEV parameters were stored in the object \texttt{A} previously): \begin{verbatim} > model=vector('numeric',length(dataset)) > for(i in 1:length(model)){ model[i]=GEV.DF(ordered[i],A$est[1],A$est[2],A$est[3])} \end{verbatim} Plotting \texttt{model} against \texttt{empirical} produces the corresponding probability plot for the set of sea surge annual maxima; this can be done using the code: \begin{verbatim} > plot(model~empirical,main='Probability plot') > abline(0,1) \end{verbatim} where the code \texttt{abline(0,1)} superimposes the line with intercept 0 and gradient 1 (i.e. the unit diagonal). The resulting plot is shown in Figure 2.3 (left-hand-side). Notice that the model--based probability estimates match up very closely with what we see in our data, suggesting that the GEV is plausible. \subsubsection*{Producing a quantile--quantile plot in \textsf{R}} \noindent In a similar fashion, we can produce a quantile plot using the following code in \textsf{R}: \begin{verbatim} > model.quantile=vector("numeric",length(dataset)) > GEV.INV=function(data,mu,sigma,xi){ if(xi==0){ INV=mu-sigma*log(-log(1-data))} else{ INV=mu+(sigma/xi)*(((-log(data))^(-xi))-1)} return(INV)} > for(i in 1:length(model.quantile)){ model.quantile[i]=GEV.INV(empirical[i],A$est[1],A$est[2],A$est[3])} > plot(model.quantile~ordered,main='Quantile plot') > abline(0,1) \end{verbatim} \noindent The function \texttt{GEV.INV} computes the inverse of the GEV distribution at each of the points in \texttt{data}; we evaluate this for each value in \texttt{empirical} to find $\hat{F}(i/(n+1))$ and store the results in the vector \texttt{model.quantile}. We then plot \texttt{model.quantile} against \texttt{ordered} to obtain the quantile plot, again superimposing the line of equality using \texttt{abline}. The resulting plot is shown in Figure 2.3 (right-hand-side). Again, this plot shows a reasonable fit of the GEV to the annual maximum sea surges at Wassaw. \begin{figure}[] \centering \psfrag{0.0}[][]{\sf \tiny 0.0} \psfrag{0.2}[][]{\sf \tiny 0.2} \psfrag{0.4}[][]{\sf \tiny 0.4} \psfrag{0.6}[][]{\sf \tiny 0.6} \psfrag{0.8}[][]{\sf \tiny 0.8} \psfrag{1.0}[][]{\sf \tiny 1.0} \psfrag{7}[][]{\sf \tiny 7} \psfrag{8}[][]{\sf \tiny 8} \psfrag{9}[][]{\sf \tiny 9} \psfrag{10}[][]{\sf \tiny 10} \psfrag{11}[][]{\sf \tiny 11} \psfrag{12}[][]{\sf \tiny 12} \psfrag{13}[][]{\sf \tiny 13} \psfrag{P}[][]{\sf \scriptsize Probability plot} \psfrag{MQ}[][]{\sf \tiny model.quantile} \psfrag{O}[][]{\sf \tiny ordered} \psfrag{QQ}[][]{\sf \scriptsize Quantile plot} \psfrag{E}[][]{\sf \tiny Empirical} \psfrag{M}[][]{\sf \tiny Model} \psfrag{9}[][]{\sf \tiny 9} \psfrag{10}[][]{\sf \tiny 10} \psfrag{11}[][]{\sf \tiny 11} \psfrag{12}[][]{\sf \tiny 12} \psfrag{13}[][]{\sf \tiny 13} \psfrag{year}[][]{\sf \tiny Year} \includegraphics[width=17cm,height=8cm]{was2.ps} \caption{Probability plot (left) and quantile--quantile plot (right) for the GEV fitted to annual maximum sea--surges observed at Wassaw Island.} \end{figure} \subsection{Return level estimation} In question 2 of Section 1.4, you were asked to provide an estimate of the height of a new sea wall to protect the city of Savannah against the storm surge we would expect to see \begin{itemize} \item [(i)] once in ten years; \item [(ii)] once in a hundred years. \end{itemize} Using the data, and the data alone, we could not obtain a meaningful estimate of (ii) because we only have 50 years of data. However, we can now use our fitted GEV to extrapolate beyond the range of our data to estimate such return levels. Using Equation (2.10), we find that \[ \hat{z}_{10} = 8.711-\frac{1.311}{0.108}\left[\left(-\textrm{log}\left(1-10^{-1}\right)\right)^{0.108}-1\right] = 11.33 \textrm{ feet}; \] similarly, \[ \hat{z}_{10} = 8.711-\frac{1.311}{0.108}\left[\left(-\textrm{log}\left(1-100^{-1}\right)\right)^{0.108}-1\right] = 13.46 \textrm{ feet}. \] Thus, an estimate of the height of the sea--wall might be about 11.5 feet or 13.5 feet to protect against the once in ten year, or once in a hundred year, storm surges (respectively). In fact, due to the invariance property of maximum likelihood estimators, our estimates of $z_{10}$ and $z_{100}$ are also the maximum likelihood estimates of these quantities. \noindent As with our inference for the GEV parameters, it is preferable to quote estimates of return levels with their estimated standard errors. Since $z_{r}$ is a function of the GEV parameters $\mu$, $\sigma$ and $\xi$, we can use the \textit{delta method} to find the approximate variance of $\hat{z}_{r}$. Specifically, \begin{eqnarray*} \textrm{Var}(\hat{z}_{r}) &\approx& \nabla z_{r}^{T}V\nabla z_{r}, \end{eqnarray*} where $V$ is the variance--covariance matrix of $(\hat{\mu},\hat{\sigma},\hat{\xi})^{T}$ and \begin{eqnarray*} \nabla z_{r}^{T} &=& \left[\frac{\partial z_{r}}{\partial \mu},\frac{\partial z_{r}}{\partial \sigma},\frac{\partial z_{r}}{\partial \xi} \right] \\ &=&\left[1, -\xi^{-1}(1-y_{r}^{-\xi}), \sigma \xi^{-2}(1-y_{r}^{-\xi})-\sigma \xi^{-1}y_{r}^{-\xi}\textrm{log }y_{r}\right], \end{eqnarray*} where $y_{r}=-\textrm{log}(1-r^{-1})$, evaluated at $(\hat{\mu},\hat{\sigma},\hat{\xi})$. \subsubsection*{Using \textsf{R} to obtain standard errors for return levels} \noindent Recall that $V$ is stored in the matrix \texttt{varcovar} in \textsf{R}. Also, recall that our estimates of the GEV parameters are stored in \texttt{A\$est}. We can define $\nabla z_{10}$, for example, in \textsf{R} as \begin{verbatim} > y10=-log(1-(1/10)) > del=matrix(ncol=1,nrow=3) > del[1,1]=1 > del[2,1]=-((A$est[3])^(-1))*(1-(y10^(-A$est[3]))) > del[3,1]=((A$est[2])*((A$est[3])^(-2))*(1-((y10)^(-A$est[3])))) -((A$est[2])*((A$est[3])^(-1))*((y10)^(-(A$est[3])))*log(y10)) > del.transpose=t(del) \end{verbatim} \noindent Then the \textsf{R} command for matrix multiplication -- \texttt{\%*\%} -- can be used to obtain an estimate of the standard error for $\hat{z}_{10}$ in the following way: \begin{verbatim} > sqrt(del.transpose%*%varcovar%*%del) [,1] [1,] 0.3614568 \end{verbatim} \noindent Estimated standard errors for other return levels can be obtained in a similar way. For example, for the standard error for the 100--year return level, we would replace \begin{verbatim} > y10=-log(1-(1/10)) \end{verbatim} with \begin{verbatim} > y100=-log(1-(1/100)) \end{verbatim} and then \texttt{y10} would be replaced with \texttt{y100} throughout. Table 2.2 shows a range of estimated return levels, with associated standard errors in parentheses. Of course, we could use these standard errors to construct confidence intervals for our return level estimates; however, as we shall shortly discuss, such confidence intervals can be misleading. \begin{table}[!h] \centering \begin{tabular}{|c| c c c c|} \hline &\multicolumn{4}{|c|}{Return level} \\ & $z_{10}$ &$z_{100}$ &$z_{200}$ &$z_{1000}$ \\ \hline MLE (st. err.) & 11.33 (0.361)&13.46 (0.938) &13.99 (1.182)&15.09 (1.821) \\ \hline \end{tabular} \caption{Estimated return levels for sea surges at Wassaw Island, with standard errors (in parentheses).} \end{table} \subsection{Using the \textsf{ismev} package in \textsf{R}} \noindent \noindent The \textsf{ismev} package provides an add--on suite of functions for \textsf{R} users to support the examples and computations in \textit{An Introduction to Statistical Modeling of Extreme Values} (Coles, 2001); the acronym of this book gives the package its name. Once installed, \textsf{ismev} gives the user access to the datasets used in Coles (2001), and the functions supplied by this package allow the user to reproduce the examples used in the book as well as perform rudimentary extreme value analyses on other datasets. All routines in the \textsf{ismev} package assume the data (where necessary) have been pre--processed so that, for example, the user has a set of block maxima ready to analyse. We shall give some thought to the practicalities of obtaining a set of block maxima and dealing with incomplete datasets in future computer practical sessions; for now, we assume any necessary data pre--processing has taken place. \noindent Installing the \textsf{ismev} package requires a version of \textsf{R} no older than \textsf{R} 1.50. Provided you have \textsf{R} 1.50, or something more recent, installed on your computer, you can load \textsf{ismev} using the command: \begin{verbatim} > library(ismev) \end{verbatim} Doing so will give you access to functions which can be used to perform a basic analysis of annual maxima in \textsf{R}: for example, \texttt{gev.fit} and \texttt{gev.diag}. The command \texttt{gev.fit} provides a maximum likelihood fitting routine for the GEV that also allows generalised linear modelling of each parameter (we will return to this in Chapter 5). For example, with our set of annual maximum sea surges stored in \texttt{wassaw}, typing \begin{verbatim} > gev.fit(wassaw) \end{verbatim} gives the following output: \newpage \begin{verbatim} $conv [1] 0 $nllh [1] 89.52412 $mle [1] 8.7114765 1.3114046 -0.1083845 $se [1] 0.2094549 0.1489799 0.1074819 \end{verbatim} If the convergence code, \texttt{conv}, takes the value zero, then successful convergence has been achieved, as is the case in our example above. \texttt{nllh} gives the value of the minimised negative log--likelihood, which has been found to be 89.524 -- exactly the value we obtained when doing this ``from first principles'' using the \texttt{nlm} routine in Section 2.2.2. The output also shows the maximum likelihood estimates, with their standard errors, giving (to three decimal places): \begin{displaymath} \hat{\mu} = 8.711 \textrm{ (0.209)} \qquad \hat{\sigma} = 1.311 \textrm{ (0.149)} \qquad \hat{\xi} = -0.108 \textrm{ (0.107)} \end{displaymath} Again, this is (almost!) exactly what we obtained when fitting the GEV `from first principles' in Section 2.2.2; any differences can be attributed to the numerical procedure used to optimise the problem. \noindent To check the goodness--of--fit of the GEV we can use the \textsf{ismev} function \texttt{gev.diag}. First of all, we need to store the output from \texttt{gev.fit} in \texttt{B}, for example: \begin{verbatim} > B=gev.fit(wassaw) \end{verbatim} Then, typing \begin{verbatim} > gev.diag(B) \end{verbatim} gives the plots shown in Figure 2.4. You should notice that the plots in the top row (probability plot and quantile plot) are exactly the same as those we constructed, from "first principles", in Figure 2.3. However, \texttt{gev.diag} also gives a "return level plot" (bottom left) and a plot of the fitted GEV density superimposed onto a histogram of the actual data. The return level plot is really just a quantile plot, but the $x$--axis is drawn on a logarithmic scale so as to compress the tail of the fitted GEV distribution and thus focus attention on long--period return levels, extrapolated far beyond the range of our observed data. Attached to these return level estimates are 95\% confidence bounds using $\textrm{estimate}\pm 1.98 \times \textrm{s.e.}$, where the standard errors have been obtained via the delta method (as in Section 2.2.4). Superimposed are empirical quantiles for the data we are analysing, for comparison. Both this plot, and the plot of the data with the fitted density, seem to indicate the adequacy of our fitted GEV. \begin{figure} \centering \psfrag{0.0}[][]{\tiny \sf 0.0} \psfrag{0.2}[][]{\tiny \sf 0.2} \psfrag{0.4}[][]{\tiny \sf 0.4} \psfrag{0.6}[][]{\tiny \sf 0.6} \psfrag{0.8}[][]{\tiny \sf 0.8} \psfrag{1.0}[][]{\tiny \sf 1.0} \psfrag{7}[][]{\tiny \sf 7} \psfrag{8}[][]{\tiny \sf 8} \psfrag{9}[][]{\tiny \sf 9} \psfrag{10}[][]{\tiny \sf 10} \psfrag{11}[][]{\tiny \sf 11} \psfrag{12}[][]{\tiny \sf 12} \psfrag{13}[][]{\tiny \sf 13} \psfrag{14}[][]{\tiny \sf 14} \psfrag{0.00}[][]{\tiny \sf 0.00} \psfrag{0.05}[][]{\tiny \sf 0.05} \psfrag{0.10}[][]{\tiny \sf 0.10} \psfrag{0.15}[][]{\tiny \sf 0.15} \psfrag{0.20}[][]{\tiny \sf 0.20} \psfrag{0.25}[][]{\tiny \sf 0.25} \psfrag{1e-01}[][]{\tiny \sf 0.1} \psfrag{1e+00}[][]{\tiny \sf 1} \psfrag{1e+01}[][]{\tiny \sf 10} \psfrag{1e+02}[][]{\tiny \sf 100} \psfrag{1e+03}[][]{\tiny \sf 1000} \psfrag{Model}[][]{\tiny \sf Model} \psfrag{Empirical}[][]{\tiny \sf Empirical} \psfrag{z}[][]{\tiny \sf z} \psfrag{f(z)}[][]{\tiny \sf f(z)} \psfrag{Return Level}[][]{\tiny \sf Return level} \psfrag{Return Period}[][]{\tiny \sf Return period} \psfrag{Probability Plot}[][]{\scriptsize \sf Probability plot} \psfrag{Quantile Plot}[][]{\scriptsize \sf Quantile plot} \psfrag{Return Level Plot}[][]{\scriptsize \sf Return level plot} \psfrag{Density Plot}[][]{\scriptsize \sf Density plot} \includegraphics[width=14.5cm,height=6.5cm]{wasplot3.ps} \caption{Diagnostic plots indicating the goodness--of--fit of the GEV to the Wassaw sea surge extremes.} \end{figure} \section{Further examples} \subsection*{Example 2.3: Rainfall in the Lake District} The data shown in Table 2.3 are annual maximum rainfall accumulations, obtained from daily records, for a period of 21 years (1991--2011 inclusive) at Eskdale in the Lake District. \begin{table}[!h] \centering \begin{tabular}{|c c c c c c c c c c c|} \hline 333 &213 &790 &343 &351 &521 &307 &305 &352 &277 &319\\ 319 &339 &262 &285 &297 &327 &620 &350 &545 &258 &\\ \hline \end{tabular} \caption{Annual maximum daily rainfall totals (mm) at Eskdale (1991--2011).} \end{table} \noindent Shown below are the maximum likelihood estimates of the GEV parameters, obtained using \textsf{R}; also shown is the observed information matrix $I_{O}$. \[ \hat{\mu} = 304.242; \qquad \hat{\sigma} = 68.977; \qquad \hat{\xi} = 0.249; \] \[ I_{O} = \left(\begin{array}{c c c} 0.0062 &-0.0046 &0.1962 \\ -0.0046 &0.0091 &-0.2114 \\ 0.1962 &-0.2114 &48.4122 \end{array} \right). \] \begin{itemize} \item [(a)] Estimate the standard errors for $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\xi}$, and use these to construct 95\% confidence intervals for the GEV parameters. Comment. \item [(b)] Find the correlation between $\hat{\sigma}$ and $\hat{\xi}$? \item [(c)] Estimate the 100 and 1000 year return levels for daily rainfall totals at Eskdale. \item [(d)] Obtain standard errors for your estimates in (b), and use these to construct 95\% confidence intervals in the usual way. Comment. \end{itemize} \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \textit{\bf \large Solution to Example 2.3} \huge \ding{46} \normalsize \vspace{22cm} }} \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \textit{\bf \large Solution to Example 2.3} \huge \ding{46} \normalsize \vspace{22cm} }} \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \textit{\bf \large Solution to Example 2.3} \huge \ding{46} \normalsize \vspace{22cm} }} \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \textit{\bf \large Solution to Example 2.3} \huge \ding{46} \normalsize \vspace{22cm} }} \newpage \subsection*{Example 2.4: Sunshine hours in Santiago de Compostela} Table 2.4 shows the annual maximum number of sunshine hours (per day) for Santiago de Compostela, Spain, for the years 1983--2007. \begin{table}[!h] \centering \begin{tabular}{|c c c c c c c c c|}\hline 12.8 &10.7 &11.9 &13.1 &13.3 &10.7 &11.1 &11.9 &12.1\\ 11.3 &12.5 &12.1 &12.0 &12.4 &11.8 &10.5 &10.9 &10.1 \\ 11.7 &12.7 &10.6 &11.6 &9.9 &9.6 &12.0 & & \\ \hline \end{tabular} \caption{Annual maximum sunshine hours (per day) observed at Santiago de Compostela, Spain, 1983--2007 (inclusive).} \end{table} \begin{itemize} \item [(a)] Use the \textsf{ismev} package in \textsf{R} to fit the GEV distribution to these data. Report estimates of the location, scale and shape parameters, with associated standard errors. \item [(b)] Explain why the output in part (a) shows that the Weibull (Type \textit{III} extreme value distribution) is most appropriate here, and give a practical justification for this. \end{itemize} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \textit{\bf \large Solution to Example 2.4} \huge \ding{46} \normalsize \vspace{11cm} }} \newpage \subsection*{Example 2.5: Degassing bursts at the Kilauea Volcano, Hawaii} The magnitudes of the vertical forces produced by seismic \textit{degassing bursts} are known to be associated with destructive volcanoes. Table 2.5 shows the maximum force produced by degassing bursts, every quarter in the years 2005--2011 (inclusive), for the Kilauea Volcano in Hawaii. \begin{table}[!h] \centering \begin{tabular}{|c|c c c c|} \hline &Jan--Mar &Apr--Jun &Jul--Sep &Oct--Dec \\ \hline 2005 &99983 &100067 &99905 &100367 \\ 2006 &99980 &99970 &100086 &99988 \\ 2007 &99912 &100084 &100432 &99921 \\ 2008 &100123 &99913 &100240 &100448 \\ 2009 &100162 &100296 &100015 &100122 \\ 2010 &99930 &99997 &100710 &100136 \\ 2011 &99946 &99918 &99962 &100052 \\ \hline \end{tabular} \caption{Maximum vertical forces produced (kilograms) by Seismic degassing bursts at the Kilauea Volcano in Hawaii.} \end{table} These data were stored the in the vector \texttt{bursts} in \textsf{R}, and the following output obtained: \begin{verbatim} gev.fit(haw) $conv [1] 0 $nllh [1] 178.2502 $mle [1] 9.998028e+04 8.710824e+01 5.921945e-01 $se [1] 21.6185212 21.3598264 0.3027163 \end{verbatim} \begin{itemize} \item [(a)] Use this information to estimate the degassing bursts we would expect to see (i) once every year, (ii) once every fifty years, and (iii) once every 100 years at this volcano (do not attempt to obtain standard errors for these estimates). How does your estimate in (i) compare to it's empirical counterpart? \item [(b)] Use Figure 2.5 (overleaf) to assess the goodness--of--fit of the GEV to these data. \end{itemize} \newpage \begin{figure} \centering \psfrag{Probability Plot}[][]{\sf \tiny Probability Plot} \psfrag{Empirical}[][]{\sf \tiny Empirical} \psfrag{Model}[][]{\sf \tiny Model} \psfrag{0.0}[][]{\sf \tiny 0.0} \psfrag{0.2}[][]{\sf \tiny 0.2} \psfrag{0.4}[][]{\sf \tiny 0.4} \psfrag{0.6}[][]{\sf \tiny 0.6} \psfrag{0.8}[][]{\sf \tiny 0.8} \psfrag{1.0}[][]{\sf \tiny 1.0} \psfrag{Quantile Plot}[][]{\sf \tiny Quantile Plot} \psfrag{100000}[][]{\sf \tiny 100000} \psfrag{100200}[][]{\sf \tiny 100200} \psfrag{100400}[][]{\sf \tiny 100400} \psfrag{100600}[][]{\sf \tiny 100600} \psfrag{100800}[][]{\sf \tiny 100800} \psfrag{Return Level Plot}[][]{\sf \tiny Return Level Plot} \psfrag{1e-01}[][]{\sf \tiny 1e-01} \psfrag{1e+00}[][]{\sf \tiny 1e+00} \psfrag{1e+01}[][]{\sf \tiny 1e+01} \psfrag{1e+02}[][]{\sf \tiny 1e+02} \psfrag{1e+03}[][]{\sf \tiny 1e+03} \psfrag{Return Period}[][]{\sf \tiny Return Period} \psfrag{Return Level}[][]{\sf \tiny Return Level} \psfrag{z}[][]{\sf \tiny z} \psfrag{f(z)}[][]{\sf \tiny f(z)} \psfrag{0.000}[][]{\sf \tiny 0.000} \psfrag{0.001}[][]{\sf \tiny 0.001} \psfrag{0.002}[][]{\sf \tiny 0.002} \psfrag{0.003}[][]{\sf \tiny 0.003} \psfrag{0.004}[][]{\sf \tiny 0.004} \psfrag{0.005}[][]{\sf \tiny 0.005} \psfrag{104000}[][]{\sf \tiny 104000} \psfrag{108000}[][]{\sf \tiny 108000} \includegraphics[width=15cm,height=7cm]{haw.ps} \caption{Diagnostic plots to assess the fit of the GEV to the degassing burst dataset in Table 2.5.} \end{figure} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \textit{\bf \large Solution to Example 2.5} \huge \ding{46} \normalsize \vspace{14cm} }} \newpage \section{Profile likelihood} \noindent As we saw in our solution to Example 2.3, standard errors for long range return levels, obtained via the delta method, can often be so large that confidence intervals become difficult to work with, or even meaningless. Actually, it turns out that constructing confidence intervals in the standard way ($\textrm{estimate}\pm1.96 \times \textrm{s.e.}$ for a 95\% confidence interval, for example) is not advisable for return levels. This is because of the severe asymmetry often observed in the likelihood surface for the return level, suggesting that the assumption of normality may not be valid for $z_{r}$. \noindent An alternative, and often more accurate, method for making inferences on a particular parameter can be found using the \textit{profile log--likelihood}. Formally, the log--likelihood for a parameter vector $\boldsymbol{\theta}$ can be written as $\ell(\boldsymbol{\theta}_{j},\boldsymbol{\theta}_{-j})$, where $\boldsymbol{\theta}_{-j}$ corresponds to all components of $\boldsymbol{\theta}$ excluding $\boldsymbol{\theta}_{j}$. The \textit{profile} log--likelihood for $\boldsymbol{\theta}_{j}$ is defined as \begin{eqnarray*} \ell_{p}(\boldsymbol{\theta}_{j}) &=& \textrm{max}_{\boldsymbol{\theta}_{-j}} \ell(\boldsymbol{\theta}_{j},\boldsymbol{\theta}_{-j}). \end{eqnarray*} Thus, for each value of $\boldsymbol{\theta}_{j}$ the profile log--likelihood is the maximised log--likelihood with respect to all other components of $\boldsymbol{\theta}$. For return levels using the GEV model, there are three main steps for the implementation of profile likelihood: \begin{itemize} \item [(1)] Re--parameterise the GEV model so that $z_{r}$ becomes one of the model parameters. For example, re--arrange equation (2.10) to write $\mu$ in terms of $\sigma$, $\xi$ and $z_{r}$, i.e.: \begin{eqnarray} \mu &=& z_{r}-\frac{\sigma}{\xi}\left[\left(-\textrm{log}\left(1-r^{-1}\right)\right)^{-\xi}-1\right]; \end{eqnarray} then obtain an expression for the log--likelihood $\ell(\sigma,\xi,z_{r})$ by substitution of (2.12) into \[ -m\textrm{log}\sigma - (1+1/\xi)\sum_{i=1}^{m}\textrm{log}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+} - \sum_{i=1}^{m}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}, \] where $m$ is the number of block maxima in the analysis. \item [(2)] For some fixed lower value of $z_{r}=z_{r,\textrm{low}}$, maximise the GEV log--likelihood \\ $\ell(\sigma,\xi,z_{r}=z_{r,\textrm{low}})$ with respect to the two remaining parameters ($\sigma$ and $\xi$), to obtain $\ell_{p}(z_{r})$ at $z_{r,\textrm{low}}$. Here, in terms of the more general notation above, $\boldsymbol{\theta}_{j} = z_{r}$ and $\boldsymbol{\theta}_{-j}=(\sigma,\xi)$. \item [(3)] Repeat step (2) for a range of values $z_{r}$ for $z_{r,\textrm{low}}\leq z_{r}\leq z_{r,\textrm{up}}$, and plot $\ell_{p}(z_{r})$ against $z_{r}$ to show the profile log--likelihood curve for $z_{r}$. \end{itemize} Figure 2.6 (left) shows a plot of the profile log likelihood for the 100 year return level $\ell_{p}(z_{100})$ for the Lake District rainfall data in Example 2.3, using $z_{100,\textrm{low}}=550$ and $z_{100,\textrm{up}} = 2900$. The plot to the right--hand--side is for the 1000--year return level. Both plots reveal strong asymmetry in the (profile) log--likelihood for the return levels, and it should be clear from these plots that constructing return levels in the usual way (as we did in Example 2.3, part (c)) will be misleading. So how can we use these plots? \noindent In this example, we are partitioning the GEV parameter vector $\boldsymbol{\theta}=(z_{r},\sigma,\xi)$ into two components $(\boldsymbol{\theta}^{(1)},\boldsymbol{\theta}^{(2)})$, where $\boldsymbol{\theta}^{(1)}=z_{r}$ and $\boldsymbol{\theta}^{(2)}=(\sigma,\xi)$, and the profile log--likelihood is now defined as \begin{eqnarray*} \ell_{p}(\boldsymbol{\theta}^{(1)}) &=& \textrm{max}_{\boldsymbol{\theta}^{(2)}}(\boldsymbol{\theta}^{(1)},\boldsymbol{\theta}^{(2)}). \end{eqnarray*} The following result leads to a procedure for making inferences on the maximum likelihood estimator of $\boldsymbol{\theta}^{(1)}$. \subsection*{Result} \noindent Let $x_{1}, \ldots, x_{n}$ be independent realisations from a distribution within a parametric family $\mathcal{F}$, and let $\hat{\boldsymbol{\theta}}_{0}$ be the maximum likelihood estimator of the $d$--dimensional model parameter $\boldsymbol{\theta}_{0}=(\boldsymbol{\theta}^{(1)},\boldsymbol{\theta}^{(2)})$, where $\boldsymbol{\theta}^{(1)}$ is a $k$--dimensional subset of $\boldsymbol{\theta}_{0}$. Then, under suitable regularity conditions, for large $n$ \begin{eqnarray*} D_{p}(\boldsymbol{\theta}^{(1)}) &=& 2 \left\{\ell(\hat{\boldsymbol{\theta}}_{0}) - \ell_{p}(\boldsymbol{\theta}^{(1)})\right\} \sim \chi_{k}^{2}. \end{eqnarray*} Thus, for our single component $\boldsymbol{\theta}^{(1)}=z_{r}$, the set of values $C_{\alpha}$ for which $\left\{z_{r}: D_{p}(z_{r})\leq c_{\alpha}\right\}$ provides a $(1-\alpha)$ confidence interval for $z_{r}$, where $c_{\alpha}$ is the $(1-\alpha)$ quantile of the $\chi_{1}^{2}$ distribution. \noindent We apply this result to the profile log--likelihood for $z_{100}$ and $z_{1000}$; a cut--off point equal to $\frac{1}{2}\times \chi_{1}^{2}(0.05)=1.921$ is shown in Figure 2.6 by the horizontal broken line. The points of intersection between this line and the profile log--likelihood for $z_{r}$ define the set of values $C_{0.05}$ giving a more appropriate 95\% confidence interval for $z_{r}$ than obtained using the usual approach. These points of intersection are also shown on Figure 2.6 (dotted vertical lines), and give confidence interval for $z_{100}$ and $z_{1000}$ based on the profile likelihood approach of $(602.5, 2765)$mm and $(800,103000)$mm respectively. These intervals are shown in the table below, along with the intervals we obtained the standard way in Example 2.3. \vspace{0.5cm} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \begin{tabular}{|c|c|c|} \hline &$\hat{z}_{100}$ &$\hat{z}_{1000}$ \\ \hline Standard 95\% CI &\color{white}{Lee Fawcett is no good, ch}&\textcolor{white}{Lee Fawcett is no good, ch} \\ & & \\ \hline Profiled 95\% CI &\color{white}{Lee Fawcett is no good, ch}&\textcolor{white}{Lee Fawcett is no good, ch} \\ & & \\ \hline \end{tabular} %\caption{Lake district rainfall data: Comparison of standard and profiled 95\% confidence intervals for the 100 and 1000 year return levels.} %\end{table} \vspace{0.4cm} {\large \textbf{Comments}} \vspace{2.65cm} }} \newpage \begin{figure} \centering \psfrag{-129}[][]{\sf \tiny --129} \psfrag{-128}[][]{\sf \tiny --128} \psfrag{-127}[][]{\sf \tiny --127} \psfrag{-126}[][]{\sf \tiny --126} \psfrag{-127.5}[][]{\sf \tiny --127.5} \psfrag{-127.0}[][]{\sf \tiny --127.0} \psfrag{-126.5}[][]{\sf \tiny --126.5} \psfrag{-126.0}[][]{\sf \tiny --126.0} \psfrag{-125.5}[][]{\sf \tiny --125.5} \psfrag{P}[][]{\sf \tiny Profile log--likelihood} \psfrag{500}[][]{\sf \tiny 500} \psfrag{1000}[][]{\sf \tiny 1000} \psfrag{1500}[][]{\sf \tiny 1500} \psfrag{2000}[][]{\sf \tiny 2000} \psfrag{2500}[][]{\sf \tiny 2500} \psfrag{4000}[][]{\sf \tiny 4000} \psfrag{6000}[][]{\sf \tiny 6000} \psfrag{8000}[][]{\sf \tiny 8000} \psfrag{10000}[][]{\sf \tiny 10000} \psfrag{12000}[][]{\sf \tiny 12000} \psfrag{Z100}[][]{\sf \scriptsize $z_{100}$} \psfrag{Z1000}[][]{\sf \scriptsize $z_{1000}$} \includegraphics[width=17.5cm,height=9cm]{profilelakes.ps} \caption{Profile log--likelihood curves for the 100 (left) and 1000 (right) year return levels for daily rainfall accumulations in Eskdale, the Lake District.} \end{figure} The plots shown above in Figure 2.6 can be produced in \textsf{R} using the \texttt{gev.prof} function within the \textsf{ismev} package. For example, suppose the results of fitting the GEV to the Lake District rainfall data (obtained using \texttt{gev.fit}) are stored in \texttt{A}. Then the code: \begin{verbatim} > gev.prof(A,xlow=550,xup=2900,100) \end{verbatim} should produce something like the plot for the 100 year return level shown in Figure 2.6. We will see more of this function in a forthcoming computer practical session. \section{Models for minima} Let $\bar{M}_{n}=\textrm{min}\left\{X_{1}, \ldots, X_{n}\right\}$. If we can assume the $X_{i}$ are independent and identically distributed, we can apply similar arguments to $\bar{M}_{n}$ as we applied to $M_{n}$ in Section 2.1.2. \noindent If there exist sequences of constants $a_{n}>0$ and $b_{n}$ such that, as $n \to \infty$, \begin{eqnarray*} \textrm{Pr}\left\{(\bar{M}_{n}-b_{n})/a_{n}\leq x\right\} &\to& \bar{G}(x) \end{eqnarray*} for some non--degenerate distribution $\bar{G}$, then $\bar{G}$ is a member of the GEV family of distributions for minima: \begin{eqnarray*} \bar{G}(x; \bar{\mu}, \sigma,\xi) &=& 1-\textrm{exp}\left\{-\left[1-\xi\left(\frac{x-\bar{\mu}}{\sigma}\right)\right]_{+}^{-1/\xi}\right\}, \end{eqnarray*} where $a_{+}=\textrm{max}(0,a)$, $-\infty<\bar{\mu}<\infty$, $\sigma>0$ and $-\infty<\xi<\infty$. This result can be useful where we are interested in modelling extremely small, rather than extremely large, observations (e.g. annual minimum air temperatures). Alternatively, we could negate our set of block minima and then model the corresponding set of maxima, giving identical maximum likelihood estimates of the GEV parameters but for the sign correction $\hat{\bar{\mu}}=-\hat{\mu}$. \section{The GEV: Words of warning} \vspace{0.5cm} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \vspace{15cm}}}