\chapter{Dependent extremes} \section{Introduction} The threshold--based approach to modelling extremes, as discussed in Chapter 3, has obvious practical advantages over the more traditional `block maxima' approach of Chapter 2. Since extremes are, by their very nature, scarce, a modelling procedure which allows the inclusion of more data in the analysis has got to be a good thing. Indeed, the whole point of a threshold--based analysis is that we include \textit{all} extremes in the analysis, extreme in the sense that the observations used have all exceeded some pre--determined threshold $u$. At a practical level, we might expect a threshold based analysis to result in estimates of return levels with much reduced standard errors and, correspondingly, narrower confidence intervals. \\ \\ However, using information on \textit{all} extremes, as natural and appealing as it may sound, brings with it it's own set of problems -- not least, the issue of \textit{temporal dependence}, sometimes referred to as \textit{serial correlation}. For example, if we are analysing extremes of temperature, and we have hourly records, then we might expect dependence from hour to the next. Such dependence is also likely to be present in the extremes of the process. For example, in a threshold analysis, we are likely to include all hourly observations from a heatwave, and these extreme temperatures are also very likely to depend on one another from one hour to the next. This, of course, seems obvious, and is common to the extremes of most environmental time series (e.g. hourly wind speeds during a storm, sea--surge measurements during a hurricane...). And this is what makes such datasets interesting! However, such temporal dependence violates one of the most basic modelling assumptions we have made in this course so far -- the assumption that our series of extremes are independent. In fact, this is the first ``I'' of ``IID'' -- we will look at the ``ID'' part in the next chapter! \newpage \section{Example: Hourly wind speeds at High Bradfield} \noindent Hourly maximum wind gusts (in knots) were collected at High Bradfield, in the Peak District, over a period of 10 years from January 1st 2003 to December 31st 2012, giving total of 87,672 observations (including any missing values; 81,835 after the missing values have been removed). Suppose these hourly observations are in the vector \texttt{brad} in \textbf{\sf R}; the series without the missing values has been stored in the vector \texttt{brad2}. A plot of the partial autocorrelation function can be useful in ascertaining the extent of any temporal dependence in a series. Similarly, a plot of each observation against it's neighbour can be produced to check for the presence of serial correlation. In \textbf{\sf R}: \begin{verbatim} > pacf(brad2) > plot(brad[1:87671],brad[2:87672] \end{verbatim} produces the plots shown in Figure 4.1. The partial autocorrelation function shows significant autocorrelation between observations one and two hours apart (lags 1 and 2); however, this plot doesn't focus on the extremes of the series. Neither does the plot of the time series against the version at lag 1 (right hand side) -- however, it is clear from this plot that there is dependence between consecutive extremes even above a high threshold (horizontal and vertical lines in the plot). This could cast doubt on the validity of the GPD when fitted to all exceedances of this threshold. \\ \\ Figure 4.2 delves a bit deeper; here, we see a plot of a portion of 150 observations from the wind speed series at Bradfield. Again, a possible threshold of 50 knots has been superimposed. It is clear that extreme wind gusts tend to occur in \textit{clusters}; this makes perfect sense, as a storm might last for several hours with sustained gusts. How can we proceed? Ignoring this dependence, and proceeding with a fit of the GPD to the set of all threshold exceedances is an option, and a reasonable fit to the data might still be achieved. However, there are now fewer independent observations in the data than the likelihood assumes, and there is a real chance that the standard errors attached to estimates of the GPD parameters, and associated return levels, will be under--estimated. \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{0}[][]{\sf \tiny 0} \psfrag{10}[][]{\sf \tiny 10} \psfrag{20}[][]{\sf \tiny 20} \psfrag{30}[][]{\sf \tiny 30} \psfrag{40}[][]{\sf \tiny 40} \psfrag{50}[][]{\sf \tiny 50} \psfrag{60}[][]{\sf \tiny 60} \psfrag{80}[][]{\sf \tiny 80} \psfrag{100}[][]{\sf \tiny 100} \psfrag{PACF}[][]{\sf \scriptsize Partial autocorrelations} \psfrag{Lag}[][]{\sf \scriptsize Lag} \psfrag{Xt}[][]{\sf \tiny $X_{t}$} \psfrag{Xt+1}[][]{\sf \tiny $X_{t+1}$} \includegraphics[width=16cm,height=7cm]{brad.ps} \caption{Partial autocorrelation function (left) and plot of the hourly wind gusts against the next value in the series (right).} \end{figure} \begin{figure} \centering \psfrag{20}[][]{\sf \tiny 20} \psfrag{30}[][]{\sf \tiny 30} \psfrag{40}[][]{\sf \tiny 40} \psfrag{50}[][]{\sf \tiny 50} \psfrag{60}[][]{\sf \tiny 60} \psfrag{70}[][]{\sf \tiny 70} \psfrag{80}[][]{\sf \tiny 80} \psfrag{0}[][]{\sf \tiny 0} \psfrag{50}[][]{\sf \tiny 50} \psfrag{100}[][]{\sf \tiny 100} \psfrag{150}[][]{\sf \tiny 150} \psfrag{w}[][]{\sf \scriptsize Wind speed (knots)} \psfrag{t}[][]{\sf \scriptsize Time} \includegraphics[width=16cm,height=12cm]{brad2.ps} \caption{Time series plot of a small section of the Bradfield wind speed data, with threshold.} \end{figure} \section{Maxima of dependent (but stationary) series} \noindent The book by Leadbetter \textit{et al.} (1983) considers, in great detail, properties of extremes of dependent processes. A key result often used is `Leadbetter's $D(u_{n})$ condition', which ensures that long--range dependence is sufficiently weak so as not to affect the asymptotics of an extreme value analysis. This condition is stated more formally in the Definition below. \newpage {\large\textbf{Definition}\textit{ (Leadbetter's $D(u_{n})$ condition)}} \vspace{0.3cm} \noindent A stationary series $\tilde{X}_{1}, \tilde{X}_{2}, \ldots$ is said to satisfy the $D(u_{n})$ condition if, for all $i_{1}<\ldots l$, \begin{eqnarray} \Big\arrowvert\textrm{Pr}\left\{\tilde{X}_{i_{1}}\leq u_{n}, \ldots, \tilde{X}_{i_{p}} \leq u_{n}, \tilde{X}_{j_{1}} \leq u_{n}, \ldots, \tilde{X}_{j_{q}}\leq u_{n}\right\} \nonumber \\ -\textrm{Pr}\left\{\tilde{X}_{i_{1}}\leq u_{n}, \ldots, \tilde{X}_{i_{p}}\leq u_{n}\right\}\textrm{Pr}\left\{\tilde{X}_{j_{1}}\leq u_{n}, \ldots,\tilde{X}_{j_{q}}\leq u_{n}\right\}\Big\arrowvert&\leq& \alpha(n,l), \end{eqnarray} where $\alpha(n,l) \to 0$ for some sequence $l_{n}$ such that $l_{n}/n \to 0$ as $n \to \infty$. \\ \\ \noindent For sequences of independent variables, the difference in probabilities in the above expression is exactly zero for \textit{any} sequence $u_{n}$. More generally, we will require that the $D(u_{n})$ condition holds only for a specific sequence of thresholds $u_{n}$ that increases with $n$. For such a sequence, the $D(u_{n})$ condition ensures that, for sets of variables that are far enough apart, the difference in probabilities expressed in (4.1), while not zero, is sufficiently close to zero to have no effect on the limit laws for extremes. \vspace{0.5cm} {\large \textbf{Theorem} \textit{(Extremes of dependent sequences)}} \vspace{0.3cm} \noindent Let $\tilde{X}_{1}, \tilde{X}_{2}, \ldots$ be a stationary series satisfying Leadbetter's $D(u_{n})$ condition, and let $\tilde{M}_{n}=\textrm{max}\{\tilde{X}_{1}, \ldots, \tilde{X}_{n}\}$. Now let $X_{1}, X_{2}, \ldots$ be an \textit{independent} series with $X$ having the same distribution as $\tilde{X}$, and let $M_{n}=\textrm{max}\{X_{1}, \ldots, X_{n}\}$. Then if $M_{n}$ has a non--degenerate limit law given by $\textrm{Pr}\left\{(M_{n}-b_{n})/a_{n} \leq x\right\} \to G(x)$, it follows that \begin{eqnarray} \textrm{Pr}\left\{(\tilde{M}_{n}-b_{n})/a_{n} \leq x\right\} \to G^{\theta}(x) \end{eqnarray} for some $0 \leq \theta \leq 1$. \noindent The parameter $\theta$ is known as the \textit{extremal index}, and quantifies the extent of extremal dependence: $\theta=1$ for a completely independent process, and $\theta \to 0$ with increasing levels of (extremal) dependence. Since $G$ in the above theorem is necessarily an extreme value distribution, and due to the \textit{max--stability} property (see Leadbetter \textit{et al.}, 1983), then the distribution of maxima in processes displaying short--range temporal dependence (characterised by the extremal index $\theta$) is also a GEV distribution; the powering of the limit distribution by $\theta$ only affects the location and scale parameters of this distribution. \\ \\ \noindent The above theorem implies that if maxima of a stationary series converge -- which, from Chapter 2, we know they will do -- then, provided an appropriate $D(u_{n})$ condition is satisfied, the limit distribution is related to the limit distribution of an independent series. The effect of dependence, as seen in expression (4.2), is just a replacement of $G$ as the limit distribution with $G^{\theta}$. In fact, if $G$ corresponds to the GEV distribution with parameters $(\mu, \sigma, \xi)$, then \begin{eqnarray*} G^{\theta}(x) &=& \textrm{exp}\left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]^{-1/\xi}\right\}^{\theta} \\ &=& \textrm{exp}\left\{-\left[1 + \xi\left(\frac{x-\mu^{*}}{\sigma^{*}}\right)\right]^{-1/\xi} \right\}, \end{eqnarray*} where $\mu^{*}=\mu - \frac{\sigma}{\xi}\left(1-\theta^{-\xi}\right)$ and $\sigma^{*}=\sigma \theta^{\xi}$. Thus, if the (approximate) distribution of $M_{n}$ is GEV with parameters $(\mu,\sigma,\xi)$, then the (approximate) distribution of $\tilde{M}_{n}$ is GEV with parameters $(\mu^{*}, \sigma^{*},\xi)$. %\subsubsection{Toy example} %Let $Y_{0}, Y_{1}, Y_{2}, \ldots$ be an independent sequence of random variables with distribution function \subsection{Modelling block maxima} \noindent Provided long--range dependence is weak, we can proceed to model block maxima from series with short--range extremal dependence as outlined in Chapter 2, since the distribution of block maxima falls within the same family of distributions as would be appropriate if the series were truly independent. This is fantastic news! Short--range temporal dependence is a much more plausible assumption than complete independence, and our modelling approach is still valid! However, the main difference -- excluding the change in parameters from $(\mu, \sigma, \xi)$ to $(\mu^{*}, \sigma^{*}, \xi)$ -- is that our implied $n$ (the number we are taking the maxima over) is now effectively reduced due to the dependence, so convergence of maxima to the limit distribution will be slower. And shouldn't we be using threshold methods anyway, which use information on \textit{all} extremes and not just those that are the maximum within their block? \subsection{Modelling threshold exceedances} \noindent Though the modelling procedure for fitting the GEV to a set of annual maxima is unchanged for series which display short--term temporal dependence, some revision is needed of the threshold exceedance approach. If all threshold exceedances are used in our analysis, and the GPD fitted to the set of threshold excesses, the likelihoods we use will be incorrect since they assume independence of sample observations. In practice, several techniques have been developed to circumvent this problem, including: \begin{itemize} \item [\textbf{1.}] filtering out an (approximately) independent set of threshold exceedances \item [\textbf{2.}] fitting the GPD to \textit{all} exceedances, ignoring dependence, but then appropriately adjusting the inference (usually an inflation of standard errors) to take into account the reduction in information \item [\textbf{3.}] Explicitly modelling the temporal dependence in the process \end{itemize} \noindent Though the first approach above is by far the most widely--used, research has recently focussed on the relative merits of the other two approaches. The third approach makes use of multivariate extreme value theory, and so we shall re--visit this idea in more detail in Chapter 6. For now, let us consider the first approach which is probably the most commonly--used. \subsubsection{Declustering} Since the mid--1990s, various methods for \textit{declustering} a series of extremes, to extract a set of independent extremes, have been discussed in the literature. The most natural, commonly--used method of declustering is that of \textit{runs declustering}. This is how it works: \begin{itemize} \item [1.] Choose an auxiliary `declustering parameter' (which we call $\kappa$) \item [2.] A cluster of threshold excesses is then deemed to have terminated as soon as at least $\kappa$ consecutive observations fall below the threshold \item [3.] Go through the entire series identifying clusters in this way \item [4.] The maximum (or `peak') observation from each cluster is then extracted, and the GPD fitted to the set of cluster peak excesses. \end{itemize} This approach is often referred to as the \textit{peaks over threshold} approach (POT, Davison and Smith, 1990) and is widely accepted as the main pragmatic approach for dealing with clustered extremes. Although this approach is quite easy to implement, there are issues surrounding the choice of $\kappa$; if \begin{itemize} \item $\kappa$ is too small, the cluster peaks will not be far enough apart to safely assume independence \item $\kappa$ is too large, there will be too few cluster exceedances on which to form our inference \end{itemize} \noindent It has also been shown (Fawcett and Walshaw, 2012\footnote{Estimating return levels from serially dependent extremes, \textit{Environmetrics} \textbf{23}(3), pp 272--283}) that parameter estimates can be sensitive to the choice of $\kappa$, and $\kappa$ is all too often chosen arbitrarily. \newpage \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} Let us return to the wind speeds at High Bradfield. Suppose we implement runs declustering with (a) $\kappa=2$; (b) $\kappa=4$ and (c) $\kappa=10$. Identify the clusters associated with (a), (b) and (c) on the snapshots below (top to bottom). \begin{center} \psfrag{20}[][]{\sf \tiny 20} \psfrag{30}[][]{\sf \tiny 30} \psfrag{40}[][]{\sf \tiny 40} \psfrag{50}[][]{\sf \tiny 50} \psfrag{60}[][]{\sf \tiny 60} \psfrag{70}[][]{\sf \tiny 70} \psfrag{80}[][]{\sf \tiny 80} \psfrag{0}[][]{\sf \tiny 0} \psfrag{50}[][]{\sf \tiny 50} \psfrag{100}[][]{\sf \tiny 100} \psfrag{150}[][]{\sf \tiny 150} \psfrag{w}[][]{\sf \scriptsize Wind speed (knots)} \psfrag{t}[][]{\sf \scriptsize Time} \includegraphics[width=16cm,height=6cm]{brad2.ps} \vspace{0.5cm} \includegraphics[width=16cm,height=6cm]{brad2.ps} \vspace{0.5cm} \includegraphics[width=16cm,height=6cm]{brad2.ps} \end{center} }} \newpage Suppose we use $\kappa=10$. Obviously, we wouldn't want to identify clusters by hand for the full Bradfield wind speed series (recall that we have 10 years of hourly observations!). Unfortunately, there is no function in \textsf{ismev} to perform runs declustering. I have written the following \textbf{\sf R} code to perform this declustering, for $\kappa=10$: \begin{verbatim} cluster10=function(dataset,threshold){ x=list() z=list() j=1 { for(i in (11):length(dataset)){ if(dataset[i-10]>threshold & dataset[i-9]<=threshold & dataset[i-8]<=threshold & dataset[i-7]<=threshold & dataset[i-6]<=threshold & dataset[i-5]<=threshold & dataset[i-4]<=threshold & dataset[i-3]<=threshold & dataset[i-2]<=threshold & dataset[i-1]<=threshold & dataset[i]<=threshold){ x=max(dataset[j:i]) ifelse(i !=length(dataset), j<-i+1, NA) z=c(z,x)}}} return(z)} \end{verbatim} Now, in \textbf{\sf R}, typing \begin{verbatim} > cluster.peaks=as.numeric(cluster10(brad2,50)) \end{verbatim} stores the set of cluster peaks in the vector \texttt{cluster.peaks}. \subsubsection{Fitting the model} Now we can fit the GPD to this set of cluster peak excesses using the function \texttt{gpd.fit} in \textsf{ismev}: \begin{verbatim} A=gpd.fit(cluster.peaks,50) $threshold [1] 50 $nexc [1] 299 $conv [1] 0 \end{verbatim} \newpage \begin{verbatim} $nllh [1] 992.3284 $mle [1] 11.9363320 -0.1607165 $rate [1] 1 $se [1] 0.8970203 0.0488930 \end{verbatim} Thus our GPD, fitted to the set of cluster peak excesses, has \[ \hat{\sigma}=11.936 (0.897) \quad \textrm{ and } \quad \hat{\xi}=-0.161 (0.049) \] (standard errors in parentheses). Notice, however, that the function has returned an estimate of the threshold exceedance rate (which we called $\lambda_{u}$ in the last chapter) of 1; this is because now, all of the data we are supplying \texttt{gpd.fit} with are above the threshold (recall that we are using the set of extracted cluster peaks!). Thus, we can put this `right' by typing \begin{verbatim} > A$rate = length(cluster.peaks)/length(brad2) \end{verbatim} giving \[ \hat{\lambda}_{u,\textrm{cp}} = 0.00365 (0.0002), \] again with the standard error in parentheses; the subscript notation on $\hat{\lambda}$ -- $u,\textrm{cp}$ -- is to remind us that this is the rate of cluster peak exceedances above the threshold, rather than the full rate of threshold exceedances. \subsubsection{Return level estimation} We can now produce estimates of return levels in the usual way; that is, for the GPD we would use equation (3.9) from Chapter 3, but now the threshold exceedance rate $\hat{\lambda}_{u}$ must be replaced with the \textit{cluster peak} exceedance rate $\hat{\lambda}_{u, \textrm{cp}}$ -- which in this example is 0.00365, or 0.365\%. \\ \\ So, for example, an estimate of the 50 year gust speed at High Bradfield, would be \[ 50+\frac{11.936}{-0.161}\left[(50\times \{365.25\times24\}\times 0.00365)^{-0.161}-1\right] = 101.533 \textrm{ knots.} \] We can also use profile likelihood to obtain a 95\% confidence interval for this estimate -- in the usual way, using the \textsf{ismev} function \texttt{gpd.prof}: \begin{verbatim} > gpd.prof(A,50,npy=365.25*24,xlow=92,xup=130) > abline(v=101.533,lty=2) \end{verbatim} This gives the plot shown in Figure 4.3; the corresponding 95\% confidence interval is approximately (94.7 knots, 117 knots). \begin{figure}[!h] \centering \psfrag{-993}[][]{\sf \tiny -993} \psfrag{-994}[][]{\sf \tiny -994} \psfrag{-995}[][]{\sf \tiny -995} \psfrag{-996}[][]{\sf \tiny -996} \psfrag{-997}[][]{\sf \tiny -997} \psfrag{-998}[][]{\sf \tiny -998} \psfrag{100}[][]{\sf \tiny 100} \psfrag{110}[][]{\sf \tiny 110} \psfrag{120}[][]{\sf \tiny 120} \psfrag{130}[][]{\sf \tiny 130} \psfrag{P}[][]{\sf \scriptsize Profile log--likelihood} \psfrag{Return}[][]{\sf \scriptsize Return level} \includegraphics[width=16cm,height=12cm]{bradprof.ps} \caption{Profile log--likelihood as a result of a cluster peaks analysis of the wind speed data at High Bradfield.} \end{figure} \newpage \section{Analyses of cluster peaks: words of warning} \vspace{0.5cm} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \vspace{16.5cm}}}