\chapter{Multivariate extremes} \section{Introduction} In this section we consider the problems we face if we wish to model the extremal behaviour of two or more (dependent) processes simultaneously. There are several reasons why we may wish to do this: \begin{itemize} \item to model the extreme behaviour of a particular variable over several nearby locations (e.g. rainfall over a network of sites -- simultaneous flooding at several locations could cause devastation); \item to model the joint extremes of two or more different variables at a particular location (e.g. wind and rain at a site -- the combined effects of wind and rain during a hurricane can result in extreme storm surge); \item to model the joint behaviour of extremes which occur as consecutive observations in a time--series (e.g. consecutive hourly maximum wind gusts during a storm). \end{itemize} All of these problems suggest fitting an appropriate limiting multivariate distribution to the relevant data. However, as we shall see, the derivation of such a multivariate distribution is not as easy as we might hope. The analogy with the Normal distribution as a model for means breaks down as we move into $n$ dimensions! It is not even clear what the `relevant data' should be! Most of the increased complexity is apparent in the move from $1$ to $2$ dimensions, so we will focus largely on bivariate problems. \section{Componentwise maxima models} \subsection{Example: Rainfall at 8 locations in Scotland} Suppose we want to study the joint extremes of daily rainfall accumulations at the network of 8 sites is southwest Scotland shown in Figure 6.1. \vspace{-1cm} \begin{figure}[!h] \begin{center} \includegraphics[height=6cm,width=12cm]{Daveloc.ps} \end{center} \caption{Eight rainfall recording stations in southwest Scotland.} \end{figure} \noindent Such issues are of great interest, especially currently -- given the severe flooding experienced in the U.K. Suppose we have sequences of daily total rainfall at each location. There is likely to be strong inter--site dependence in extremes, in the sense that days with heavy rain are liable to occur simultaneously across locations. The raw multivariate observations are $8$--dimensional vectors of the daily rainfall over the eight sites. Now suppose we wish to take a block--maxima approach, with `blocks' being years. For any given year, the 8--dimensional vector of annual maxima is unlikely to be one of the raw multivariate observations. Let's simplify to the bivariate case. Let $(X_1, Y_1), (X_2, Y_2), \ldots$ be i.i.d. vectors with distribution function $F(x,y)$. Now consider the componentwise block maxima \[ M_{x,n} = \max_{i=1,\ldots,n}\{X_i\} \qquad \mbox{and} \qquad M_{y,n} = \max_{i=1,\ldots,n}\{Y_i\}. \] We define the {\it vector of componentwise maxima} to be \[ {\bf M_n} = (M_{x,n}, M_{y,n}). \] ${\bf M_n}$ is not necessarily one of the original observations $(X_i, Y_i)$. Nevertheless, we are interested in the limiting behaviour of ${\bf M_n}$ as $n \rightarrow \infty$. The first point to note is that standard univariate extreme value results apply in each margin. When considering the dependence, this allows us to make a simplifying assumption. We assume that the $X_i$ and $Y_i$ variables have a known marginal distribution. It is convenient to assume a unit Fr\'{e}chet distribution (see Chapter 2), which has CDF \[ F(z) = \textrm{exp}(-1/z), \qquad z>0. \] This gives rise to a very simple normalisation of maxima: \[ \textrm{Pr} (X_i < x) = \textrm{Pr} (M_{x,n}/n < x) = \exp(-1/x), \qquad x>0, \] (and similarly for $Y_i$). So if we consider the re--scaled vector \[ {\bf M_n^*} = \left(\max_{i=1,\ldots,n}\{X_i\}/n, \max_{i=1,\ldots,n}\{Y_i\}/n\right), \] the margins are unit Fr\'{e}chet for all $n$, and hence we can characterise the limiting joint behaviour of ${\bf M_n^*}$ without having to worry about the margins. Unfortunately no limiting parametric family exists (for bivariate extremes, or multivariate extremes in general)! \subsection{Theorem: limiting distributions for bivariate extremes} Let ${\bf M_n^*} = (M_{x,n}^*, M_{y,n}^*)$ be the normalised maxima as above, where the $(X_i,Y_i)$ are i.i.d. with standard Fr\'{e}chet marginal distributions. Then if \[ \textrm{Pr}(M_{x,n}^*, M_{y,n}^*) \rightarrow G(x,y), \] where $G$ is non--degenerate, then $G$ has the form \begin{equation} G\left(x,y\right)=\exp\left\{-V\left(x,y\right)\right\};\quad x>0,\;\;y>0 \end{equation} where: \begin{equation} V\left(x,y\right)=2\int_{0}^{1}\max\left(\frac{\omega}{x},\frac{1-\omega}{y}\right)dH\left(\omega\right) \end{equation} and $H$ is a distribution function on $\left[0,1\right]$ satisfying the mean constraint: \begin{equation} \int_{0}^{1}\omega\;dH\left(\omega\right)=0.5. \end{equation} Hence the class of bivariate extreme value distributions is in one-to-one correspondence with distribution functions $H$ satisfying the constraint (6.3). If $H$ is differentiable with density $h$, then (6.1) becomes \begin{equation*} V\left(x,y\right)=2\int_{0}^{1}\max\left(\frac{\omega}{x},\frac{1-\omega}{y}\right)h(\omega)d\omega. \end{equation*} However some simple models arise when $H$ is not differentiable. For example, if $H$ places mass $0.5$ on each of $\omega=0$ and $\omega=1$, then we get \[ G(x,y) = \exp\{-(x^{-1}+y^{-1})\}, \qquad x>0, y>0, \] corresponding to independent $x$ and $y$. Since the GEV provides the complete class of marginal limit distributions, then the complete class of bivariate extreme value distributions is obtained as follows. If we suppose $X$ and $Y$ are GEV with parameters $(\mu_x, \sigma_x, \xi_x)$ and $(\mu_y, \sigma_y, \xi_y)$ respectively, then the transformations \[ \tilde{x}=\left[1+\xi_x\left(\frac{x-\mu_x}{\sigma_x}\right)\right]^{1/\xi_x} \quad \mbox{and} \quad \tilde{y}=\left[1+\xi_y\left(\frac{y-\mu_y}{\sigma_y}\right)\right]^{1/\xi_y} \] obtain unit Fr\'{e}chet margins. Hence \[ G(x,y) = \exp\{-V(\tilde{x}, \tilde{y})\} \] is a bivariate extreme value distribution with the appropriate margins for valid $V(.)$, provided $[1+\xi_x(x-\mu_x)/\sigma_x]>0]$ and $[1+\xi_y(x-\mu_y)/\sigma_y]>0]$. %\subsection{Modelling bivariate extremes in practice} %In practice, modelling usually involves identifying a parametric sub--family with appropriate flexibility to handle the structure inherent in the data. Models can be fitted, e.g. by maximum--likelihood estimation, either in two steps (marginal components followed by dependence function), or in a single sweep. All of these procedures, including the choice of models, are handled in a very similar way when dealing with threshold exceedances. We consider the details in the next section. \section{Bivariate threshold excess models} \subsection{Example: wave--surge data at Newlyn, Cornwall} Here, we choose a different type of example of dependence to the rainfall problem considered in Section 6.2.1. Specifically, we consider two variables recorded concurrently at the same site. A series of 3--hourly measurements on sea--surge were obtained from Newlyn, southwest England, giving -- at each time point -- measurements of the wave height (in metres) and the surge (also in metres). Figure 6.2 below shows these two variables plotted against each other; this plot suggests a tendency for extremes of one variable to coincide with extremes of the other. This dependence could be important -- the impact of an event that is simultaneously extreme in both variables is likely to be much greater than if extremes of either component occurred in isolation. \begin{figure}[!h] \begin{center} \includegraphics[height=12cm,width=8cm,angle=-90]{wavesurge.ps} \end{center} \caption{Concurrent wave and surge heights at Newlyn, Cornwall.} \end{figure} \subsection{From componentwise maxima to threshold excesses} Now we want to define our bivariate extremes as those observations which exceed a threshold in one or other margin. For our bivariate observation $(X,Y)$, let's focus on $X$. We have already seen that the distribution function for the exceedances of a threshold $u$ by a variable $X$, conditional on $X>u$ for large enough $u$, is given by: \begin{equation*} H(x)=1-\lambda_{u_{x}}\left\{ 1+\frac{\xi_{x}\left( x-u_{x} \right)}{\sigma_{x}} \right\}^{-1/\xi_{x}} \end{equation*} defined on $\left\{x-u_{x}: x-u_{x}>0 \text{ and }\left(1+\xi_{x}\left(x-u_{x}\right)/\sigma_{x}\right)>0\right\}$, where $\xi_{x}\neq0$, $\sigma_{x}>0$, and $\lambda_{u_{x}}=\textrm{Pr}\left(X>u_{x}\right)$. Now we can obtain a unit Fr\'{e}chet margin with the transformation: \begin{equation} \tilde{x}=-\left(\log\left\{1-\lambda_{u_{x}}\left[1+\frac{\xi_{x}\left(x-u_{x}\right)}{\sigma_{x}}\right]^{-1/\xi_{x}}\right\}\right)^{-1}. \end{equation} If we apply the analogous transformation to the $Y$ margin, we obtain \begin{equation*} G\left(x,y\right)=\exp\left\{-V\left(\tilde{x},\tilde{y}\right)\right\};\quad x>u_{x},\;\;y>u_{y}, \end{equation*} where $V(x,y)$ is as defined in Equation (6.1), again satisfying the mean constraint in Equation (6.3). \section{Modelling bivariate extremes in practice} In practice, modelling usually involves identifying a parametric sub--family with appropriate flexibility to handle the structure inherent in the data. Models can be fitted, e.g. by maximum--likelihood estimation, either in two steps (marginal components followed by dependence function), or in a single sweep. All of these procedures, including the choice of models, are handled in a very similar way when dealing with either bivariate componentwise maxima or bivariate threshold exceedances. \subsection{Modelling the dependence structure} The class of bivariate extreme value models contains many families of distributions which can be used to model the dependence structure in the data. The dependence structure must satisfy the conditions on $H\left(\omega\right)$. Possible choices are: \begin{itemize} \item Logistic Model --- symmetric\\[-8mm] \item Negative Logistic Model\\[-8mm] \item Bilogistic Model --- asymmetric\\[-8mm] \item Dirichlet Model\\[-6mm] \end{itemize} Here we will focus on the logistic model and the bilogistic model as two commonly used but contrasting choices. \vspace{0.5cm} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} \normalsize \vspace{3.5cm}}} \subsubsection{1. The Logistic model} Here, for $V(x,y$) in Equation (6.1), we have \begin{equation*} \left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}, \end{equation*} where $x>0$, $y>0$ and $\alpha\in\left(0,1\right)$. \begin{itemize} \item $\alpha\rightarrow1$ corresponds to independent variables. \item $\alpha\rightarrow0$ corresponds to perfectly dependent variables. \item This model is symmetric --- the variables are \textit{exchangeable}. \end{itemize} \subsubsection{2. The Bilogistic model} Now we have the following form for $V(x,y)$: \begin{equation*} -x\gamma^{1-\alpha}-y\left(1-\gamma\right)^{1-\beta}, \end{equation*} where $0<\alpha<1$, $0<\beta<1$ and $\gamma=\gamma\left(x,y;\alpha,\beta\right)$ is the solution of: \begin{equation*} \left(1-\alpha\right)x\left(1-\gamma\right)^{\beta}=\left(1-\beta\right)y\gamma^{\alpha} \end{equation*} \begin{itemize} \item Independence is obtained when $\alpha=\beta\rightarrow1$ and when one of $\alpha$ or $\beta$ is fixed and the other approaches 1. \item When $\alpha=\beta$ the model reduces to the logistic model. \item The value of $\alpha-\beta$ determines the extent of \textit{asymmetry} in the dependence structure. \end{itemize} \subsection{Likelihood calculations} After transformation to unit Fr\'{e}chet margins, we can obtain the probability density function of the chosen dependence model by differentiation of Equation (6.1) to give $g(x,y)$; from this, the likelihood can be formed (and maximised) in the usual way. However, inference for the bivariate threshold excess setup is complicated by the fact that a bivariate pair may exceed a specified threshold in just one of its components. For example, consider Figure 6.3 below, an updated version of Figure 6.2 for the wave--surge data now including marginal thresholds (obtained by using mean residual life plots for both components). We obtain contributions to the likelihood function for a pair of observations in the following way, where $\theta$ represents the parameter(s) in our dependence model (e.g. $\theta=\alpha$ if we use the logistic model; $\theta=(\alpha,\beta)$ in the bilogistic model): \begin{equation*} g\left(x,y;\theta \right)= \begin{cases} \left.\displaystyle \frac{\partial^{2} G}{\partial x \partial y}\right|_{\left(\tilde{x},\tilde{y}\right)}\quad&\text{if}\left(x,y\right)\in\text{Region 1}\\ \left.\displaystyle\frac{\partial G}{\partial x}\right|_{\left(\tilde{x},\tilde{u}_{y}\right)}\quad&\text{if}\left(x,y\right)\in\text{Region 2}\\ \left.\displaystyle\frac{\partial G}{\partial y}\right|_{\left(\tilde{u}_{x},\tilde{y}\right)}\quad&\text{if}\left(x,y\right)\in\text{Region 3}\\ G\left(\tilde{u}_{x},\tilde{u}_{y}\right)\quad&\text{if}\left(x,y\right)\in\text{Region 4} \end{cases} \end{equation*} Then, we have \[ L(\theta;x,y) = \prod_{i=1}^{n}g(\tilde{x}_{i},\tilde{y}_{i}). \] %\begin{center} % \pnode(-10em,-8ex){A}% % \pnode(-5em,-10ex){B}% % \pnode(9em,-8ex){C}% % \pnode(4em,-10ex){D}% % \pnode(-10em,-18ex){E}% % \pnode(-5em,-20ex){F}% % \pnode(9em,-18ex){G}% % \pnode(4em,-20ex){H}% % \makebox(0,0){\hspace{-25em}\vspace{-15ex}Region 3}% % \makebox(0,0){\ncline[linecolor=magenta,linestyle=solid,arrows=->]{A}{B}}% % \makebox(0,0){\hspace{23em}\vspace{-15ex}Region 1}% % \makebox(0,0){\ncline[linecolor=magenta,linestyle=solid,arrows=->]{C}{D}}% % \makebox(0,0){\hspace{-25em}\vspace{-35ex}Region 4}% % \makebox(0,0){\ncline[linecolor=magenta,linestyle=solid,arrows=->]{E}{F}}% % \makebox(0,0){\hspace{23em}\vspace{-35ex}Region 2}% % \makebox(0,0){\ncline[linecolor=magenta,linestyle=solid,arrows=->]{G}{H}}% %\end{center} %\begin{figure}[!h] % \begin{center} % \vspace{-0.4cm} % \rotatebox{-90}{\scalebox{0.25}{\includegraphics{wavethresh}}} % \end{center} % \end{figure} \begin{figure}[!h] \begin{center} \psfrag{Wave Height (m)}[c][c]{\small \sf Wave Height (m)} \psfrag{Surge (m)}[c][c]{\small Surge (m)} \vspace{-0.4cm} \scalebox{0.7}{\includegraphics{arrow}} \end{center} \caption{Threshold classification of bivariate data.} \end{figure} %\begin{figure}[!h] %\begin{center} %\includegraphics[height=12cm,width=8cm,angle=-90]{wavesurge.ps} %\end{center} %\caption{Concurrent wave and surge heights at Newlyn, Cornwall.} %\end{figure} %%\newpage \vspace{2cm} \section*{Example 6.1: Wave--surge analysis at Newlyn} Consider the wave--surge data from Newlyn, Cornwall, shown in Figures 6.2 and 6.3. Flood defences in Newlyn have been designed to withstand a sea swell resulting from, at most, a wave height of $x=9 \textrm{ metres}$ combined with a surge height of $y=0.7 \textrm{ metres}$. A threshold--based approach to modelling is to be used for the wave--surge data shown in Figures 6.2 and 6.3. Mean residual life plots suggest marginal thresholds of $u_{x}=6.1 \textrm{ metres}$ and $u_{y}=0.32 \textrm{ metres}$ (as shown in Figure 6.2) for identifying wave height and surge height as extreme. \newpage \begin{itemize} \item [(a)] Assuming extreme wave heights occur independently of extreme surge heights, find the probability that the flood defence system in Newlyn will be overwhelmed. \item [(b)] Now assume there is extremal dependence between wave height and surge height. \begin{itemize} \item [(i)] Assuming the logistic model for this dependence, obtain the likelihood contribution to $L(\alpha;x_{i},y_{i})$ if $x>u_{x}$ but $y library(ismev) > data(wavesurge) > wavesurge wave surge 1 1.50 -0.009 2 1.83 -0.053 3 2.44 -0.024 4 1.68 0.000 5 1.49 0.079 6 1.20 0.068 7 1.35 -0.009 . . . . . . . . . \end{verbatim} For convenience, we can ``pull out'' the wave and surge components: \begin{verbatim} > wave=wavesurge[,1] > surge=wavesurge[,2] \end{verbatim} \newpage \begin{itemize} \item [(a)] Thus, we can estimate $(\hat{\sigma}_{x},\hat{\xi}_{x})$ and $(\hat{\sigma}_{y}, \hat{\xi}_{y})$ using \texttt{gpd.fit}. For example: \begin{verbatim} > gpd.fit(wave,6.1) $threshold [1] 6.1 $nexc [1] 141 $conv [1] 0 $nllh [1] 155.1482 $mle [1] 1.3335727 -0.1875704 $rate [1] 0.04872149 $se [1] 0.14084399 0.06565004 \end{verbatim} Using the same command on \texttt{surge}, i.e. \texttt{gpd.fit(surge, 0.32)}, we have \[ \hat{\lambda}_{u_{x}} = 0.049 (0.004) \quad \hat{\sigma}_{x} = 1.334 (0.141) \quad \hat{\xi}_{x} = -0.188 (0.066) \] and \[ \hat{\lambda}_{u_{y}} = 0.051 (0.004) \quad \hat{\sigma}_{y} = 0.093 (0.011) \quad \hat{\xi}_{y} = -0.041 (0.080), \] with standard errors in parentheses. \end{itemize} \framebox[17cm]{ \parbox[t]{16cm}{ \huge \ding{46} \normalsize \textcolor{white}{Assuming independence, we have \begin{eqnarray*} \textrm{\sf Pr}(\textrm{\sf flood defence fails}) &=&\textrm{\sf Pr}(X>9) \times \textrm{\sf Pr}(Y>0.7) \\ &=& \hat{\lambda}_{u_{x}}\left[1+\hat{\xi}_{x}\left(\frac{9-u_{x}}{\hat{\sigma}_{x}}\right)\right]_{+}^{-1/\hat{\xi}_{x}} \times \hat{\lambda}_{u_{y}}\left[1+\hat{\xi}_{y}\left(\frac{0.7-u_{y}}{\hat{\sigma}_{y}}\right)\right]_{+}^{-1/\hat{\xi}_{y}}\\ &=&0.049\left[1-0.188\left(\frac{9-6.1}{1.334}\right)\right]_{+}^{1/0.188} \\ & & \times 0.051\left[1-0.041\left(\frac{0.7-0.32}{0.093}\right)\right]_{+}^{1/0.041} \\ &=&0.002995 \times 0.000583 \\ &=&0.00000175. \end{eqnarray*}} }} \begin{itemize} \item [(b)] \item [(i)] For the logistic model, we have \[ G(x,y) = \textrm{exp}\left\{-\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}\right\}. \] \end{itemize} \subsubsection*{Region 1:} If $x>u_{x}$ and $y>u_{y}$, then we are in Region 1 of Figure 6.3, and so \[ g(x,y; \alpha) = \displaystyle \frac{\partial^{2}G}{\partial x \partial y}. \] Differentiating w.r.t. $x$ gives \begin{eqnarray*} \frac{\partial G}{\partial x} &=& \textrm{exp}\left\{-\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}\right\} \\ & & \times -\alpha\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha-1} \times (-1/\alpha)x^{-(1/\alpha+1)} \\ & & \\ &=& \textrm{exp}\left\{-\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}\right\} \left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha-1}x^{-(1/\alpha-1)}. \end{eqnarray*} Differentiating the result w.r.t. $y$ then gives (after some algebra) \begin{eqnarray*} \frac{\partial^{2}G}{\partial x \partial y}&=&(xy)^{-(1/\alpha +1)}\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha-2}\left[\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}-\left(1-1/\alpha\right)\right] \\ & & \textrm{exp}\left\{-\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}\right\}\\ & & \\ &=&(\tilde{x}\tilde{y})^{-(1/\alpha +1)}\left(\tilde{x}^{-1/\alpha}+\tilde{y}^{-1/\alpha}\right)^{\alpha-2}\left[\left(\tilde{x}^{-1/\alpha}+\tilde{y}^{-1/\alpha}\right)^{\alpha}-\left(1-1/\alpha\right)\right]\\ & & \textrm{exp}\left\{-\left(\tilde{x}^{-1/\alpha}+\tilde{y}^{-1/\alpha}\right)^{\alpha}\right\}, \end{eqnarray*} evaluated at $x=\tilde{x}$ and $y=\tilde{y}$. \subsubsection*{Region 2:} If only $x>u_{x}$, then we are in Region 2 of Figure 6.3, and so \[ g(x,y;\alpha) = \frac{\partial G}{\partial x}. \] Thus, we have: \framebox[17cm]{ \parbox[t]{16cm}{ \huge \ding{46} \normalsize \textcolor{white}{ \begin{eqnarray*} \frac{\partial G}{\partial x} &=&\textrm{\sf exp}\left\{-\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}\right\} \times(-\alpha)\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha-1} \times (-1/\alpha)x^{-(1/\alpha+1)} \\ & & \\ &=& \textrm{\sf exp}\left\{-\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}\right\}\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha-1}x^{-(1/\alpha+1)} \\ & & \\ &=& \textrm{\sf exp}\left\{-\left(\tilde{x}^{-1/\alpha}+\tilde{u}_{y}^{-1/\alpha}\right)^{\alpha}\right\}\left(\tilde{x}^{-1/\alpha}+\tilde{u}_{y}^{-1/\alpha}\right)^{\alpha-1}\tilde{x}^{-(1/\alpha+1)} \end{eqnarray*} evaluated at $x=\tilde{x}$ and $y=\tilde{u}_{y}$. } \vspace{3cm} }} \subsubsection*{Region 3:} If only $y>u_{y}$, then we are in Region 3 of Figure 6.3, and so \[ g(x,y; \alpha) = \frac{\partial G}{\partial y}. \] Thus, we have: \framebox[17cm]{ \parbox[t]{16cm}{ \huge \ding{46} \normalsize \textcolor{white}{ \begin{eqnarray*} \frac{\partial G}{\partial y} &=&\textrm{\sf exp}\left\{-\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}\right\} \times(-\alpha)\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha-1} \times (-1/\alpha)y^{-(1/\alpha+1)} \\ & & \\ &=& \textrm{\sf exp}\left\{-\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}\right\}\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha-1}y^{-(1/\alpha+1)} \\ & & \\ &=& \textrm{\sf exp}\left\{-\left(\tilde{u}_{x}^{-1/\alpha}+\tilde{y}^{-1/\alpha}\right)^{\alpha}\right\}\left(\tilde{u}_{x}^{-1/\alpha}+\tilde{y}^{-1/\alpha}\right)^{\alpha-1}\tilde{u}_{x}^{-(1/\alpha+1)} \end{eqnarray*} evaluated at $x=\tilde{u}_{x}$ and $y=\tilde{y}$. } \vspace{3cm} }} \subsubsection*{Region 4:} If both $x$ and $y$ are sub-threshold, then we are in Region 4 of Figure 6.3. Thus, \[ g(x,y;\alpha) = G(\tilde{u}_{x},\tilde{u}_{y}) = \textrm{exp}\left\{-\left(\tilde{u}_{x}^{-1/\alpha}+\tilde{u}_{y}^{-1/\alpha}\right)^{\alpha}\right\}. \] \begin{itemize} \item [(ii)] For each pair of wave height/surge height measurements, we would now use the marginal GPD fits from part (a) to transform to unit Fr\'{e}chet margins (via Equation 6.4); then the full likelihood \[ L(\alpha; x,y) = \prod_{i=1}^{2894}g(\tilde{x}_{i},\tilde{y}_{i}) \] can be maximised in the usual way to estimate the dependence parameter $\alpha$. However, we can make use of the \texttt{fbvpot} command in \textsf{evd}: \begin{verbatim} > fbvpot(wavesurge, threshold=c(6.1,0.32), model='log') Call: fbvpot(x = wavesurge, threshold = c(6.1, 0.32), model = "log") Deviance: 2025.254 AIC: 2035.254 Dependence: 0.3023499 Threshold: 6.1 0.32 Marginal Number Above: 141 147 Marginal Proportion Above: 0.0487 0.0508 Number Above: 48 Proportion Above: 0.0166 Estimates scale1 shape1 scale2 shape2 dep 1.265706 -0.139014 0.091818 0.006741 0.763539 Standard Errors scale1 shape1 scale2 shape2 dep 0.13280 0.06886 0.01052 0.08404 0.02933 Optimization Information Convergence: successful Function Evaluations: 37 Gradient Evaluations: 8 \end{verbatim} Thus, we have \[ \hat{\lambda}_{u_{x}} = 0.049 (0.004) \quad \hat{\sigma}_{x} = 1.266 (0.133) \quad \hat{\xi}_{x} = -0.139 (0.069) \] \[ \hat{\lambda}_{u_{y}} = 0.051 (0.004) \quad \hat{\sigma}_{y} = 0.092 (0.011) \quad \hat{\xi}_{y} = 0.007 (0.084) \] and \[ \hat{\alpha} = 0.764 (0.029) \] Notice that there are slight differences in the estimates of the marginal GPD parameters compared to those we found in part (a) (although these differences are not significant). This is because the full log--likelihood $\ell(\lambda_{u_{x}},\lambda_{u_{y}},\sigma_{x},\sigma_{y},\xi_{x},\xi_{y},\alpha)$ has been optimised in a single sweep for all of the parameters \textit{simultaneously}. Although our estimated dependence parameter does not indicate the presence of strong extremal dependence (recall that as $\alpha \to 1$ we have decreasing dependence), at least we have correctly accounted for this dependence. \end{itemize} \vspace{0.3cm} \framebox[17cm]{ \parbox[t]{16cm}{ \huge \ding{46} \normalsize {\color{white} \begin{itemize} \item [(iii)] Assuming extremal dependence, we have \begin{eqnarray*} \textrm{\sf Pr}(\textrm{\sf Flood defence fails}) &=&1-G(\widetilde{9},\widetilde{0.7}) \\ \nonumber & & \\ &=&1-\textrm{\sf exp}\left\{-\left(\widetilde{9}^{-1/0.764}+\widetilde{0.7}^{-1/0.764}\right)\right\}. \end{eqnarray*} Using our estimated marginal GPD parameters from the output above, and the transformation given by Equation (6.4), we have \[ \widetilde{9} = 321.4721 \qquad \widetilde{0.7} = 1161.417; \] this gives, on substitution into the above expression, \[ \textrm{\sf Pr}(\textrm{\sf Flood defence fails})=1-\textrm{\sf exp}\left\{-\left(321.4721^{-1/0.764}+1161.417^{-1/0.764}\right)\right\} = 0.0006. \] Although this probability is small, notice that it is more than 350 times as large as the probability obtained assuming independence between wave height and surge height. Thus, assuming independence gives us a false sense of security, as it grossly underestimates the chance of the flood defence being overwhelmed. This is typical of what happens when we ignore extremal dependence! \vspace{6cm} \end{itemize}} }} \begin{itemize} \item [(c)] To check for asymmetry in the dependence structure, we can fit a bilogistic model in \textbf{\sf R} and then compare this to the fitted logistic (as the logistic is \textit{nested} within the bilogistic; when $\alpha=\beta$, the bilogistic reduces to the logistic). \begin{verbatim} > fbvpot(wavesurge,threshold=c(6.1,0.32),model="bilog") Call: fbvpot(x = wavesurge, threshold = c(6.1, 0.32), model = "bilog") Deviance: 2024.823 AIC: 2036.823 Dependence: 0.3032494 Threshold: 6.1 0.32 Marginal Number Above: 141 147 Marginal Proportion Above: 0.0487 0.0508 Number Above: 48 Proportion Above: 0.0166 Estimates scale1 shape1 scale2 shape2 alpha beta 1.28033 -0.14768 0.09074 0.01283 0.79655 0.71996 Standard Errors scale1 shape1 scale2 shape2 alpha beta 0.13597 0.06955 0.01047 0.08315 0.05209 0.07883 Optimization Information Convergence: successful Function Evaluations: 48 Gradient Evaluations: 10 \end{verbatim} Notice the similarity in the estimates of $\alpha$ and $\beta$; in fact, constructing confidence intervals for these two parameters gives: \begin{eqnarray*} \alpha &:&(0.694, 0.899) \\ \beta &:&(0.565, 0.874) \end{eqnarray*} Notice that these intervals overlap, suggesting that $\alpha \approx \beta$; thus, allowing for asymmetry in the dependence structure is probably not worthwhile. %\vspace{3.5cm}}} \end{itemize} \end{document}