\documentclass[handout]{beamer} \usepackage[latin1]{inputenc} \usepackage{amsmath} \usepackage{color} \usepackage{xcolor} \usepackage{graphicx} \usepackage{pifont} \usepackage{psfrag} \usepackage{skull} \usepackage{tikz} \usepackage{bm} \usepackage{listings} %\usepackage{enumitem}% http://ctan.org/pkg/enumitem %\usetheme{Warsaw} \usecolortheme{ncl} \usefonttheme[onlymath]{serif} \newcommand{\penaltymark}[1]{\colorbox{nclRed}{\color{white}(-#1)}} \newcommand{\vsp}{\vspace{0.3cm}} \newcommand{\vspl}{\vspace{0.15cm}} \newcommand{\firstsect}{\setbeamercolor{structure}{fg=green!60!black}} \newcommand{\sol}{\setbeamercolor{structure}{fg=purple!80!black}} \newcommand{\normal}{\setbeamercolor{structure}{fg=blue!60!green}} \newcommand{\secondsect}{\setbeamercolor{structure}{fg=black!50!purple}} \newcommand{\thirdsect}{\setbeamercolor{structure}{fg=orange!90!yellow}} \newcommand{\fourthsect}{\setbeamercolor{structure}{fg=green!50!blue}} \newcommand{\fifthsect}{\setbeamercolor{structure}{fg=red!70!blue}} \newcommand{\cover}{\setbeamercolor{structure}{fg=blue!30!black}} \title[]{MAS3918/8918\\ Topics in Statistical Modelling: Extremes} \author{Dr. Lee Fawcett, Dr. Vianey Palacios-Ramirez} \date{Semester 2, 2024/25} \definecolor{dark-gray}{gray}{0.425} \definecolor{light-gray}{gray}{0.125} \setbeamercovered{transparent=10} \begin{document} \lstset{ language=R, % the language of the code %backgroundcolor=\color{light-gray}, % choose the background color; you must add \usepackage{color} or \usepackage{xcolor}; should come as last argument basicstyle=\footnotesize, % the size of the fonts that are used for the code breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace breaklines=true, % sets automatic line breaking captionpos=b, % sets the caption-position to bottom commentstyle=\color{mygreen}, % comment style deletekeywords={...}, % if you want to delete keywords from the given language escapeinside={\%*}{*)}, % if you want to add LaTeX within your code extendedchars=false, % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8 frame=single, % adds a frame around the code keepspaces=true, % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible) keywordstyle=\color{blue}, % keyword style morekeywords={*,...}, % if you want to add more keywords to the set numbers=left, % where to put the line-numbers; possible values are (none, left, right) numbersep=5pt, % how far the line-numbers are from the code numberstyle=\tiny\color{mygray}, % the style that is used for the line-numbers rulecolor=\color{black}, % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here)) showspaces=false, % show spaces everywhere adding particular underscores; it overrides 'showstringspaces' showstringspaces=false, % underline spaces within strings only showtabs=false, % show tabs within strings adding particular underscores stepnumber=1, % the step between two line-numbers. If it's 1, each line will be numbered stringstyle=\color{mymauve}, % string literal style tabsize=2, % sets default tabsize to 2 spaces title=\lstname % show the filename of files included with \lstinputlisting; also try caption instead of title } \begin{frame} \titlepage \end{frame} \frame{ \frametitle{Contact details} \begin{center} \begin{tabular}{l l} \textbf{My name:} & Lee Fawcett \\ & \\ \textbf{My job:} & Reader in Applied Statistics, Associate Dean (Global) \\ & \\ \textbf{Office:} & Room 2.07 Herschel Building \\ & \\ \textbf{Phone:} & 0191 208 7228 \\ & \\ \textbf{Email:} & \texttt{lee.fawcett@ncl.ac.uk} \end{tabular} \end{center} } \begin{frame} \frametitle{Timetable and Administrative arrangements} \begin{itemize} \item \textbf{Lectures} for the Extremes topic are on \textcolor{purple}{\bf Mondays} at \textcolor{purple}{\bf 4pm} (KGVI LT6) and again on \textcolor{purple}{\bf Thursdays} at \textcolor{purple}{\bf 2pm} (KGVI LT1). \pause \vspace{0.2cm} \item Every fortnight, starting this week, there will also be a session on \textcolor{purple}{\bf Fridays} at \textcolor{purple}{\bf 12pm} (Herschel LT2); these will usualy be used as \textbf{problems classes}. \vspace{0.2cm} \item \textbf{Office hours}: \textcolor{purple}{\bf Wednesdays 2pm--4pm}. \end{itemize} \end{frame} \frame{ \frametitle{Assessment and deadlines} {\scriptsize %\begin{center} \begin{tabular}{|c|l|l|} \hline Assessment &Weighting &Date/Deadline \\ \hline \hline Homework 1 &6\% (3\% from each topic) &4pm, Friday 21st February \\ Homework 2 &7\% (3.5\% from each topic) &4pm, Friday 7th March \\ Homework 3 &7\% (3.5\% from each topic) &4pm, Friday 28th March \\ \hline 2.5-hour exam &80\% (40\% from questions from each topic) &May/June 2025 (TBC) \\ \hline \end{tabular} %\end{center} } \vsp For this topic, there \underline{may} be the option to switch out the three assignments for a single project -- stay tuned and I'll let you know soon! } \frame{ \frametitle{Other stuff} All lecture notes will be available in advance on the MAS3918/MAS8918 Canvas page. \vsp Notes will be available in two formats: \begin{itemize} \item HTML screen-reader-friendly, \textcolor{purple}{\bf accessible lecture notes} \item \textcolor{purple}{\bf Traditional PDF lecture notes} for printing: \begin{itemize} \item [--] With gaps, for students who like to write stuff down during lectures \item [--] With the gaps filled in (\textcolor{red}{in red font}) \end{itemize} \end{itemize} } \frame{ \frametitle{Other stuff} \begin{itemize} \item A summarised verion of the lecture notes will be used as \textcolor{purple}{\bf slides} in the lectures \item These slides will be made available via Canvas after the lecture, along with the associated \textcolor{purple}{\bf ReCap recordings} \item \textcolor{purple}{\bf Solutions} to problem-solving exercises will be made available via Canvas after the deadline for each assessment has passed \item One last thing... There is an excellent \textcolor{purple}{\bf course textbook} available to buy, or to loan from the University library: \begin{center} \textit{An Introduction to Statistical Modeling of Extreme Values}, \end{center} by Stuart Coles. \end{itemize} } \frame{ \frametitle{History of this course} \begin{itemize} \item MAS369: \textcolor{purple}{\bf 1995--2000} -- Dr. Mark Dixon -- now the multi--millionaire founder of sports betting company \textsc{atass} \pause \vspl \item MAS369: resurrected in \textcolor{purple}{\bf 2002} for one year only -- Dave Walshaw (see OHP slides...)\pause \vspl \item Course: \textcolor{purple}{\bf Reading module} for Stage 3 MMath/MMathStat students (and select BSc students!) \pause \vspl \item New Stage 4 module in \textcolor{purple}{\bf 2011}: material in MAS369 out--of--date with current research -- a more advanced ``topic'' in MAS8306 \vspl \item Now we want to bring some of this material to stage 3 students -- very timely, given the current climate emergency \end{itemize} } \frame{ \frametitle{Aims of this course} \begin{itemize} \item To expose you to an area of \textcolor{purple}{\bf active research} in the School \vspl \item To expose you to a very \textcolor{purple}{\bf hands--on}, applied area of Statistics \vspl \item To give you a \textcolor{purple}{\bf firm grounding} in the \textcolor{purple}{\bf basics} of techniques for analysing environmental data... \pause \vspl \item ... but also to bring you bang up--to--date with \textcolor{purple}{\bf current research ideas} \pause \vspl \item Quite a n\'{i}che area of research, but along the way we will have to get to grips with some more generic statistical tools \end{itemize} } \frame{ \frametitle{A word on pre--requisite knowledge} \begin{itemize} \item Statistical ``Theory'': I will not assume any knowledge past 2nd year Statistics courses (\textcolor{purple}{\bf MAS2901} and \textcolor{purple}{\bf MAS2906} in particular) \vspl \item Statistical ``Practice'': Some basic \textcolor{purple}{\bf Statistical Modelling} techniques will be used, but will be fully demonstrated if they go beyond MAS2901/MAS2902 \pause \vspl \item We will make use of some techniques from these courses (e.g. \textcolor{purple}{\bf maximum likelihood}) -- these will be quickly reviewed, as--and--when they make an appearance, but you may want to consult your old notes for further details \end{itemize} } \frame{ \frametitle{} \Large \begin{center} \textit{Let's start with some (non--mathematical) background and motivation...} \end{center} } \normal \frame{ \vspace{1cm} \huge \textbf{Chapter 1} \vspace{1cm} \Large \textbf{Introduction} } \frame{ \frametitle{Practical Motivation} Finally, there is almost a \textcolor{purple}{\bf global consensus} amongst scientists that our planet's climate is changing. \pause \vsp Evidence for climatic change has been collected from a variety of sources, some of which can be used to reconstruct the earth's changing climates over \textcolor{purple}{\bf tens of thousands of years}. \pause \vsp Reasonable complete global records of the earth's surface temperature since the mid--1800's indicate a \textcolor{purple}{\bf positive trend} in the average annual temperature, and maximum annual temperature, most noticeable at the earth's poles. } \frame{ \frametitle{Practical Motivation} \begin{center} \includegraphics[width=11cm,height=7.5cm]{globwarm2.ps} \end{center} } \frame{ \frametitle{Practical Motivation} Glaciers are considered the \textcolor{purple}{\bf most sensitive indicators} of climate change. \pause \vsp As the earth warms, glaciers \textcolor{purple}{\bf retreat} and \textcolor{purple}{\bf ice sheets melt}. \pause \vsp Over the last 30 years or so this has resulted in a gradual \textcolor{purple}{\bf increase in sea and ocean levels}. } { \usebackgroundtemplate{\includegraphics[height=\paperheight,width=\paperwidth]{glacer2.ps}} \setbeamertemplate{navigation symbols}{} \begin{frame}[plain] \end{frame} } \frame{ \frametitle{Practical Motivation} \begin{center} \includegraphics[width=11cm,height=7.5cm]{glacier1.ps} \end{center} } \frame{ \frametitle{Practical Motivation} \textbf{Consequences} \vspl \begin{itemize} \item Irreversible changes to \textcolor{purple}{\bf sea} and \textcolor{purple}{\bf ocean ecosystems} \pause \vspl \item Threat of \textcolor{purple}{\bf flooding} to low--lying inhabited areas of land \pause \vspl \item Effect on the earth's \textcolor{purple}{\bf weather systems} -- for example: \pause \vspl \begin{itemize} \item [--] Larger amount of warmer water in North Atlantic Ocean \pause \\ $\longrightarrow$ Stronger, more frequent tropical storms/hurricanes (e.g. Katrina, Helene, Milton) \pause \vspl \item [--] More frequent tropical cyclones in the Indian Ocean/Pacific Ocean (e.g. Bangladesh, India, Malaysia) \pause \vspl \item [--] Flooding in the U.K. -- 18-month period from October 2022 wettest on record \vspl \item [--] Named storms in the U.K. -- e.g. Storm Darragh and \textcolor{purple}{\bf Storm \'{E}owyn} \end{itemize} \end{itemize} } { \usebackgroundtemplate{\includegraphics[height=\paperheight,width=1.2\paperwidth]{../../SCHOOLS/PHIL/cc4_2022.eps}} \setbeamertemplate{navigation symbols}{} \begin{frame}[plain] \end{frame} } { \usebackgroundtemplate{\includegraphics[height=\paperheight,width=1.2\paperwidth]{../../SCHOOLS/PHIL/cc1_2022.eps}} \setbeamertemplate{navigation symbols}{} \begin{frame}[plain] \end{frame} } \frame{ \frametitle{``Toon Monsoon'', June 2012} \begin{center} \includegraphics[width=11cm,height=6.5cm]{lightning.ps} \end{center} } \frame{ \frametitle{``Toon Monsoon'', June 2012} \begin{center} \includegraphics[width=11cm,height=6.5cm]{Heworth-6.ps} \end{center} } \frame{ \frametitle{``Toon Monsoon'', June 2012} \begin{center} \includegraphics[width=11cm,height=6.5cm]{castle-stairs1.ps} \end{center} } \frame{ \frametitle{``Toon Monsoon'', June 2012} \begin{center} \includegraphics[width=4.6cm,height=3.8cm]{boatstreet.ps} \hspace{0.4cm} \includegraphics[width=4.6cm,height=3.8cm]{car.ps} \end{center} } \frame{ \frametitle{The Great North Sea Flood, 1953} \begin{center} \includegraphics[width=11cm,height=6.5cm]{northsea.ps} \end{center} } \frame{ \frametitle{The Great North Sea Flood, 2025?} \begin{center} \includegraphics[width=11cm,height=6.5cm]{londonflood.ps} \end{center} } \frame{ \frametitle{The Great Storm of 1987} \begin{itemize} \item 15/16 October 1987 \pause \vspl \item Killed 22 people in England/France \pause \vspl \item Met Office criticism: \pause \vspl \begin{itemize} \item [--] Recent cutbacks -- no weather ship in the Southwest Approaches \pause \vspl \item [--] Weather presenter Michael Fish told the British public not to worry, and that there is no ``hurricane on the way'' \end{itemize} \end{itemize} } \frame{ \frametitle{Tornadoes in Birmingham, 2005} \begin{center} \includegraphics[width=11cm,height=6.5cm]{tornado.ps} \end{center} } \frame{ \frametitle{Tornadoes in Birmingham, 2005} \begin{center} \includegraphics[width=11cm,height=6.5cm]{tornado2.ps} \end{center} } \frame{ \frametitle{Practical Motivation} The world's changing climate is also likely to give rise to an increase in the frequency of periods of intense \textcolor{purple}{\bf heat} and \textcolor{purple}{\bf drought}. \pause \begin{itemize} \item Some parts of the world are familiar with the effects of extreme heat... \pause \item ... but we are beginning to see an increase in severe hot spells in more temperate latitudes. For example, the \textcolor{purple}{\bf 2003 and 2022 European heatwaves} which resulted in thousands of deaths \pause \item Increase in devastating forest fires (\textcolor{purple}{\bf New South Wales} and \textcolor{purple}{\bf Victoria} in Australia, and of course more recently in \textcolor{purple}{\bf LA}) \pause \item Increase in famine \end{itemize} } { \usebackgroundtemplate{\includegraphics[height=\paperheight,width=1.2\paperwidth]{../../SCHOOLS/PHIL/cc3_2022.eps}} \setbeamertemplate{navigation symbols}{} \begin{frame}[plain] \end{frame} } { \usebackgroundtemplate{\includegraphics[height=\paperheight,width=1.2\paperwidth]{../../SCHOOLS/PHIL/cc5_2022.eps}} \setbeamertemplate{navigation symbols}{} \begin{frame}[plain] \end{frame} } \frame{ \frametitle{Hurricane Katrina} \begin{center} \includegraphics[width=11cm,height=6.5cm]{katrina1.ps} \end{center} } \frame{ \frametitle{Hurricane Katrina} \textbf{Some facts} \pause \vspl \begin{itemize} \item \textcolor{purple}{\bf Category 5} Hurricane (Category 3 when it made landfall) \pause \vspl \item Killed \textcolor{purple}{\bf 1833} people \pause \vspl \item Caused \textcolor{purple}{\bf \$108 bn} worth of damage \pause \vspl \item Most damage/loss of life caused by \textcolor{purple}{\bf storm surge} \pause \vspl \item Storm surge reached nearly \textcolor{purple}{\bf 14 feet} above sea level \pause \vspl \item Lowest air pressure 902 mb \pause \vspl \item Political controversy \pause \vspl \item Billed as the ``\textcolor{purple}{\bf storm of the century}'' \end{itemize} } \frame{ \frametitle{Hurricane Katrina} \begin{center} \includegraphics[width=11cm,height=6.5cm]{evacuate2.ps} \end{center} } \frame{ \frametitle{Hurricane Katrina} \begin{center} \includegraphics[width=11cm,height=6.5cm]{katrina2.ps} \end{center} } \frame{ \frametitle{Hurricane Katrina} \begin{center} \includegraphics[width=11cm,height=6.5cm]{house2.ps} \end{center} } \frame{ \frametitle{So why statistical modelling?} Statistical modelling of extreme weather has a very practical motivation: \textcolor{purple}{\bf reliability} ---\pause \vsp --- anything we build needs to have a good chance of surviving the weather/environment for the whole of its working life. \pause \vsp This has obvious implications for \textcolor{purple}{\bf civil engineers} and \textcolor{purple}{\bf planners}. } \frame{ \frametitle{So why statistical modelling?} For example, in the context of some of the examples we have looked at so far, they need to know: \pause \vspl \begin{itemize} \item How tall to build sea walls/flood defences -- for example -- \pause \begin{itemize} \item [--] the levees surrounding New Orleans or Cedar Key, Florida; \pause \item [--] flood defences in New York that were breached during Superstorm Sandy; \pause \item [--] river and sea defences in the U.K.; \pause \vspl \end{itemize} \item How strong to make buildings -- for example, in the U.K. the \textit{British Standards Institute} specifies that new structures need to be strong enough to withstand the ``\textcolor{purple}{\bf Fifty year wind speed}''; \pause \vspl \item How tall to build reservoir dams; \pause \vspl \item How much fuel to stockpile. \end{itemize} } \frame{ \frametitle{So why statistical modelling?} This motivates the need to estimate what the: \pause \begin{itemize} \item Highest tide; \item Heaviest rainfall; \item Strongest wind; \item most severe cold spell; \end{itemize} \pause will be over some fixed period of future time. } \frame{ \frametitle{So why statistical modelling?} The only sensible way to do this is to use \textcolor{purple}{\bf data} on the variable of interest (wind, rain etc.) and fit an appropriate \textcolor{purple}{\bf statistical model}.\pause \vsp The models themselves are motivated by \textcolor{purple}{\bf asymptotic theory}, and this will be our starting point in Chapter 2. \pause \vsp Before this, let's start at a very basic level, and think about why we need statistical models at all. } \frame{ \frametitle{Simple example: sea surges at Wassaw Island} Wassaw Island is a barrier island off the southeast coast of Georgia, in the USA. \pause \vsp Any hurricanes about to make landfall on the mainland often hit such barrier islands first, substantially subduing the impending storm. \pause \vsp Sea surge data is collected at Wassaw Island every hour; the data shown below are the extracted \textcolor{purple}{\bf annual maxima} for the years 1975--2024. \pause \vsp {\small %\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} } } \frame{ \frametitle{Simple example: sea surges at Wassaw Island} \begin{itemize} \item [1.] Civil engineers are interested in various \textcolor{purple}\bf {exceedance probabilities}. \vsp For example, suppose they need to know the probability that sea--surge will exceed 8.75, 11.25 and 14 feet. Use the data to estimate these probabilities. \end{itemize} } \sol \frame{ \frametitle{Solution} Estimating these probabilities using relative frequencies gives: \begin{eqnarray*} \textsf{Pr}(\textsf{sea--surge exceeds 8.75 feet}) &=& \frac{33}{50} = 0.66 \pause \\ & & \\ \textsf{Pr}(\textsf{sea--surge exceeds 11.25 feet}) &=& \frac{6}{50} = 0.12 \pause \\ & & \\ \textsf{Pr}(\textsf{sea--surge exceeds 14 feet}) &=& \frac{0}{50} = 0. \end{eqnarray*} \pause Do you \textit{really} think it's \textit{impossible} for sea--surge to exceed 14 feet? } \normal \frame{ \frametitle{Simple example: sea surges at Wassaw Island} \begin{enumerate} \item [2.] How high should state authorities build flood defences at the nearby city of Savannah if they want to protect the city against the sea--surge they would expect to see: \begin{itemize} \item [(i)] Once every 10 years; \item [(ii)] Once every hundred years? \end{itemize} \end{enumerate} } \sol \frame{ \frametitle{Solution (1/2)} In ascending numerical order, the data are: \begin{center} {\small \begin{tabular}{|c c c c c c c c c c|}\hline 6.7 &6.8 &7.3 &7.3 &7.3 &7.4 &7.6 &7.7 &7.8 &8.1 \\ 8.2 &8.3 &8.4 &8.5 &8.6 &8.7 &8.7 &8.8 &8.8 &8.9 \\ 8.9 &8.9 &8.9 &9.0 &9.1 &9.1 &9.3 &9.3 &9.4 &9.5 \\ 9.6 &9.6 &9.7 &9.7 &9.8 &10.0 &10.0 &10.2 &10.3 &10.4 \\ 10.5 &10.5 &10.8 &11.1 &11.5 &11.8 &12.2 &12.6 &12.7 &13.0 \\ \hline \end{tabular}} \pause \end{center} Since we have annual data, an intuitive argument would be to find the value which has 10\% of the data above it. \pause \vsp We have 50 observations, and 10\% of these observations exceed $11.7999 \approx 11.8$ feet. } \frame{ \frametitle{Solution (2/2)} For the level which is exceeded once every \textcolor{purple}{\bf hundred} years, we would need the value which has 1\% of the data above it. \pause \vsp We don't have enough data -- we only have 50 years of data, and we want the value that is exceeded once every \textit{hundred} years! 1\% of 50 = 1/2 an observation! } \normal \frame{ \frametitle{Further comments} We have no chance of producing any sort of estimate that goes beyond the timeframe of our collected data -- however, it is \textcolor{purple}{\bf vitally important} for engineers that we \underline{do} estimate such quantities!\pause \vsp We are trying to estimate events more extreme than have ever occurred before -- an impossible task? \pause \vsp Not \textit{impossible}... To extrapolate beyond the range of our observed data, we need a \textcolor{purple}{\bf statistical model} that we can assume is valid for our observed extremes \textit{and beyond}. \pause \vsp This is where \textcolor{purple}{\bf extreme value theory} comes in! } \normal \frame{ \vspace{1cm} \huge \textbf{Chapter 2} \vspace{1cm} \Large \textbf{Classical Models for Extremes} } \frame{ \frametitle{Background and theoretical motivation} 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$. \pause \vsp One way of characterising extremes is by considering the distribution of the \textcolor{purple}{\bf maximum order statistic} \[ M_{n}= \textrm{\sf max}\left\{X_{1}, X_{2}, \ldots, X_{n}\right\}. \tag{1} \] } \frame{ \frametitle{Background and theoretical motivation} Think back to the example concerning sea surges at Wassaw Island. \pause \vsp Sea surge measurements were taken every hour; for each year (1975--2024) you were presented with the \textcolor{purple}{\bf annual maximum} sea surge. \pause \vsp 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 (of course, for leap years the notation would be $M_{8784,i}$). } \frame{ \frametitle{Background and theoretical motivation} The assumption of IID might be reasonable here. \pause \vsp 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. \pause \vsp 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. } \frame{ \frametitle{Background and theoretical motivation} \noindent Assuming, for now, that our maxima are IID, how can we obtain the distribution of $M_{n}$? } \sol \frame{ \frametitle{Background and theoretical motivation} This is trivial (in principle), since \pause \begin{eqnarray} \textsf{Pr}\left\{M_{n} \leq x\right\}&=&\textrm{\sf Pr}\left\{X_{1}\leq x, X_{2}\leq x, \ldots, X_{n} \leq x \right\} \nonumber \\ \pause &=& \textrm{\sf Pr}\left\{X_{1} \leq x\right\} \times \textrm{\sf Pr}\left\{X_{2} \leq x \right\} \times \cdots \times \textrm{\sf Pr}\left\{X_{n} \leq x\right\} \nonumber \\ \pause &=&\left\{F(x)\right\}^{n}. \nonumber \end{eqnarray} However, in practice the distribution function $F$ is unknown. This leads to an approach based on asymptotic arguments -- specifically, we look for limiting distributions for $\{F(x)\}^{n}$ as $n \to \infty$ --- \pause this is where the field of \textcolor{purple}{\bf Extreme Value Theory} was born. } \normal \frame{ \frametitle{A brief history of Extreme Value Theory} \noindent One of the earliest books on the statistics of extreme values is \textcolor{purple}{\bf E.J. Gumbel} (1958, see Figure 2.1). \pause \vsp Gumbel traces the origins back to 1709, when \textcolor{purple}{\bf N. Bernoulli} considers the problem of estimating the age of the longest survivor in a group of people. } \frame{ \frametitle{A brief history of Extreme Value Theory} Research into extreme values as a subject in it's own right began much later, between 1920 and 1940, when work by \textcolor{purple}{\bf E.L. Dodd}, \textcolor{purple}{\bf M. Fr\'{e}chet}, \textcolor{purple}{\bf E.J. Gumbel}, \textcolor{purple}{\bf R. von Mises} and \textcolor{purple}{\bf L.H.C. Tippett} investigated the asymptotic distribution of the \textit{largest order statistic}. \pause \vsp This led to the main theoretical result: the \textcolor{purple}{\bf 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{purple}{\bf B. Gnedenko} in 1943. } \frame{ \frametitle{A brief history of Extreme Value Theory} Until 1950, development was largely theoretical. \pause \vsp In 1958, Gumbel started applying theory to problems in engineering. \pause \vsp In the 1970s, \textcolor{purple}{\bf L. de Haan} and \textcolor{purple}{\bf J. Pickands} generalised the theoretical results, giving a better basis for statistical models. } \frame{ \frametitle{A brief history of Extreme Value Theory} Since the 1980s, methods for the application of Extreme Value Theory have become much more widespread. \pause \vsp Current researchers who have played a significant role in developing applications and methodology include \textcolor{purple}{\bf Richard Smith}, \textcolor{purple}{\bf Anthony Davison}, \textcolor{purple}{\bf Jonathan Tawn} and \textcolor{purple}{\bf Stuart Coles}. } \frame{ \frametitle{A brief history of Extreme Value Theory} 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{purple}{\bf Fawcett} and \textcolor{purple}{\bf Walshaw}). %\begin{figure}[!h] \begin{center} \includegraphics[height=3.8cm]{gumbel.ps} \hspace{0.25cm} \includegraphics[height=3.8cm]{gumbelbook.ps} \hspace{0.25cm} \includegraphics[height=3.8cm]{LFbook.ps} %\caption{Emil Julius Gumbel (left), and his book ``Statistics of Extremes'' (centre); the latest book to hit the shops in August 2013 (right, pre--order on Amazon now!!).} \end{center} } \frame{ \frametitle{The Extremal Types Theorem} \noindent The obvious questions now are: \pause \vspl \begin{itemize} \item What possible distributions might be considered candidates for the distribution for $M_{n} = \{F(x)\}^{n}$ as $n \to \infty$? \pause \vspl \item Can we formulate this set of candidate distributions into a single class -- say $G$ -- which is independent of $F$? \pause \vspl \item Can we estimate the distribution of $M_{n}$ using $G$, without any reference to $F$? \end{itemize} } \frame{ \frametitle{The Extremal Types Theorem} \noindent Clearly, the limiting distribution of $M_{n}$ is \textcolor{purple}{\bf degenerate}: \pause \vspl \begin{itemize} \item The distribution converges to a single point on the real line with probability 1 \pause \vspl \item In this case, this single point is the \textcolor{purple}{\bf upper endpoint} of $F$\pause \vspl \item In some applications, this will be $\infty$ \end{itemize} } \frame{ \frametitle{The Extremal Types Theorem} This is analogous to the sample mean $\bar{X}$ converging to the population mean $\mu$ with certainty in the \textcolor{purple}{\bf Central Limit Theorem}. \pause \vsp Here, the degenerate limit is prevented by allowing a linear rescaling, so that \pause \[ \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. \pause \vsp Can we apply a similar linear rescaling to $M_{n}$ to avoid convergence of the distribution to a single point? } \frame{ \frametitle{The Extremal Types Theorem} The answer, of course, is ``yes'', and is provided by the main result in classical extreme value theory -- the \textcolor{purple}{\bf Extremal Types Theorem} -- a result for the maximum $M_{n}$ which is analogous to the Central Limit Theorem for the mean $\mu$. } \frame{ \frametitle{The Extremal Types Theorem} \noindent {\large \textbf{Theorem} \textcolor{brown}{\textit{(The Extremal Types Theorem)}}} \pause \vsp If there exist sequences of constants $a_{n}>0$ and $b_{n}$ such that, as $n \to \infty$, \[ \textrm{\sf Pr}\left\{(M_{n}-b_{n})/a_{n}\leq x\right\} \to G(x) \tag{2} \] for some non--degenerate distribution $G$, then $G$ is of the same type as one of the following distributions: \[ \mbox{\textrm{\it I}}: G(x) = \textrm{\sf exp}\left\{-\textrm{\sf exp}(-x)\right\} \qquad -\infty0, \alpha>0; \end{array} \right.\tag{4} \] \pause \[ \mbox{\textrm{\it III}}:G(x) = \left\{ \begin{array}{ll} \textrm{\sf exp}\left\{-(-x)^{\alpha}\right\} &x<0, \alpha>0 \\ 1 &x \geq 0. \end{array} \right. \tag{5} \] } \frame{ \frametitle{The Extremal Types Theorem} The three types of distribution -- $\mbox{\textrm{\it I}}$, $\mbox{\textrm{\it II}}$ and $\mbox{\textrm{\it III}}$ -- have become known as the \textcolor{purple}{\bf Gumbel}, \textcolor{purple}{\bf Fr\'{e}chet} and \textcolor{purple}{\bf Weibull} types (respectively), and are known collectively as the \textcolor{purple}{\bf extreme value distributions}.\pause \vsp For both the Gumbel and Fr\'{e}chet distributions the limiting distribution $G$ is \textcolor{purple}{\bf unbounded}; that is, the upper--endpoint tends to $\infty$. Of the two, the Fr\'{e}chet distribution gives \textcolor{purple}{\bf heavier tails}. \pause \vsp For the Weibull distribution, the limiting distribution is \textcolor{purple}{\bf bounded}. } \frame{ \frametitle{The Extremal Types Theorem} It should be noted that this Theorem does not \textit{ensure} the existence of a non--degenerate limit for $M_{n}$. \pause \vsp Nor does it specify \textit{which} of types $\mbox{\textrm{\it I}}$, $\mbox{\textrm{\it II}}$ or $\mbox{\textrm{\it III}}$ is applicable \textit{if} a limit distribution \textit{does} exist (i.e. in which \textcolor{purple}{\bf domain of attraction} the distribution of $G$ lies). \pause \vsp 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$.} } \frame{ \frametitle{The Extremal Types Theorem} \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. \pause \vsp But how do we know which one of these distributions to use? } \frame{ \frametitle{The Generalised Extreme Value distribution} \noindent In practice, working with --- and having to choose between --- three distributions is inconvenient. \pause \vsp However, there exists a parameterisation which encompasses all three types of extreme value distribution. } \frame{ \frametitle{The Generalised Extreme Value distribution} \textcolor{brown}{\bf Von Mises (1954)} and \textcolor{brown}{\bf Jenkinson (1955)} independently derived the \textcolor{purple}{\bf generalised extreme value distribution} (GEV), often denoted $\mathcal{G}(\mu, \sigma, \xi)$, with CDF: \pause \[ G(x;\mu, \sigma, \xi) = \textrm{\sf exp}\left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}\right\}, \tag{6} \] \pause where $a_{+} = \textrm{\sf max}(0,a)$. The situation where $\xi=0$ is not defined in (6), but is taken as the limit as $\xi \to 0$, given by \pause \[ G(x;\mu, \sigma)=\textrm{\sf exp}\left\{-\textrm{\sf exp}\left(-{\color{magenta}\left[}\frac{x-\mu}{\sigma}{\color{magenta}\right]}\right)\right\}. \tag{7} \] } \frame{ \frametitle{The Generalised Extreme Value distribution} The parameters $\mu$ ($-\infty<\mu<\infty$), $\sigma$ ($>0$) and $\xi$ ($-\infty<\xi<\infty$) are \textcolor{purple}{\bf location}, \textcolor{purple}{\bf scale} and \textcolor{purple}{\bf shape} parameters, respectively. \pause \vsp Shape parameter $\xi$: \pause \vspl \begin{itemize} \item $\xi=0$: \pause Gumbel (type I extreme value) distribution \pause \vspl \item $\xi>0$: \pause Fr\'{e}chet (type II extreme value) distribution \pause \vspl \item $\xi<0$: \pause Weibull (type III extreme value) distribution \pause \vspl \item Through inference for $\xi$, the data themselves determine the most appropriate type of tail behaviour -- no need for any \textit{a priori} judgements \pause \vspl \item The standard error for $\xi$ accounts for our uncertainty in choosing between the 3 EV distributions \end{itemize} } \frame{ \frametitle{The Generalised Extreme Value distribution} \begin{center} \psfrag{D}[][]{\tiny \sf Density} \psfrag{0.00}[][]{\tiny \sf 0.00} \psfrag{0.02}[][]{\tiny \sf 0.02} \psfrag{0.04}[][]{\tiny \sf 0.04} \psfrag{0.06}[][]{\tiny \sf 0.06} \psfrag{0.08}[][]{\tiny \sf 0.08} \psfrag{90}[][]{\tiny \sf 90} \psfrag{100}[][]{\tiny \sf 100} \psfrag{110}[][]{\tiny \sf 110} \psfrag{120}[][]{\tiny \sf 120} \psfrag{130}[][]{\tiny \sf 130} \psfrag{140}[][]{\tiny \sf 140} \psfrag{150}[][]{\tiny \sf 150} \includegraphics[width=8cm,height=4cm]{lecplot1.ps} \end{center} } \frame{ \frametitle{The Generalised Extreme Value distribution} But what about the \textcolor{purple}{\bf constants} $a_{n}$ and $b_{n}$? \pause \vsp 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{ \sf as } n\to \infty. \] \pause After some algebra, it turns out that \[ M_{n} \xrightarrow{\hspace{0.25cm} D \hspace{0.25cm}} \mathcal{G}(\mu^{*}, \sigma^{*}, \xi), \qquad \textrm{ \sf approximately, as }n \to \infty, \] \pause with $a_{n}$ and $b_{n}$ being absorbed into $\mu^{*}$ and $\sigma^{*}$. } \frame{ \frametitle{The Generalised Extreme Value distribution} 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}$. \pause \vsp 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. } \frame{ \frametitle{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}=\textsf{log }n$, show that the limit distribution of $(M_{n}-b_{n})/a_{n}$ is of extreme value type, and identify the distribution. } \sol \frame{ \frametitle{Example 2.1: Solution (1/2)} We want the distribution of $(M_{n}-b_{n})/a_{n}$, i.e. \pause \begin{eqnarray*} \textsf{Pr}\left\{\frac{M_{n}-b_{n}}{a_{n}} \leq z \right\} &=& \pause \textsf{Pr}\left\{\frac{\textsf{Max}(X_{1},\ldots, X_{n})-b_{n}}{a_{n}} \leq z \right\} \\ \pause & & \\ &=& \textsf{Pr}\left\{\frac{X_{1}-b_{n}}{a_{n}}\leq z, \ldots, \frac{X_{n}-b_{n}}{a_{n}}\leq z \right\} \pause \\ & & \\ &=&\textsf{Pr}\left\{X_{1}\leq a_{n}z+b_{n}\right\} \times \cdots \\ \pause & & \\ &=& \left[1-e^{-(a_{n}z+b_{n})}\right]^{n}, \end{eqnarray*} as $X_{1}, \ldots, X_{n}$ are I.I.D.} \frame{ \frametitle{Example 2.1: Solution (2/2)} Now using $a_{n}=1$ and $b_{n}=\textsf{log }n$, we get \pause \begin{eqnarray*} \left[1-e^{-(z+\textsf{log }n)}\right]^{n} &=& \pause \left[1-e^{-z}e^{\textsf{log }n^{-1}}\right]^{n} \pause \\ & & \\ &=& \left[1-n^{-1}e^{-z}\right]^{n}. \pause \end{eqnarray*} From Stage 1 methods courses, you should know that \pause \[ \textsf{exp}(y) = \textsf{lim}_{n \to \infty}\left(1+\frac{y}{n}\right)^{n}. \pause \] Thus, we have \pause \[ \left[1+\frac{-e^{-z}}{n}\right]^{n} \longrightarrow \textsf{exp}\left(-e^{-z}\right), \qquad \textsf{ as } n\to \infty. \] \pause This is the \textcolor{purple}{\bf Gumbel} (type I extreme value) distribution. } \normal \frame{ \frametitle{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. } \sol \frame{ \frametitle{Example 2.2: Solution (1/1)} Following the solution to Example 2.1, we find: \pause \begin{eqnarray*} \textsf{Pr}\left\{\frac{M_{n}-b_{n}}{a_{n}}\leq z \right\} &=& \textsf{Pr}\left\{X_{1} \leq a_{n}z+b_{n}\right\} \times \cdots \pause \\ & & \\ &=& \left[\textsf{exp}\left\{-\frac{1}{a_{n}z+b_{n}}\right\}\right]^{n}. \pause \end{eqnarray*} Letting $a_{n}=n$ and $b_{n}=0$ gives \pause \begin{eqnarray*} \left[\textsf{exp}\left\{-\frac{1}{nz}\right\}\right]^{n} &=& e^{-1/z}, \pause \end{eqnarray*} which is the \textcolor{purple}{\bf Fr\'{e}chet} (type II extreme value) distribution with $\alpha=1$, i.e. a unit Fr\'{e}chet distribution. } \normal \frame{ \frametitle{Typical application -- 1: Data pre--processing} \pause \begin{itemize} \item Choose your \textcolor{purple}{\bf block length} $n$ (usually the number of observations in a calendar year) \pause \vspl \item Discard all but the largest observation within each block \pause \vspl \item $n$ too small: limiting arguments will not hold \pause \vspl \item $n$ too large: not enough maxima to work with! \end{itemize} } \frame{ \frametitle{Typical application -- 2: Parameter estimation} \pause Fit the GEV to your set of block maxima $M_{n,i}$ -- numerical maximum likelihood estimation is the most common approach here. \pause \vsp Assuming independence, we form the likelihood in the usual way: \pause \[ L(\mu,\sigma,\xi;\boldsymbol{x}) = \prod_{i=1}^{m}g(x_{i};\mu,\sigma,\xi), \pause \] where $g$ is the GEV probability density function and can be found, after differentiation of the distribution function (6), to be \pause \[ \frac{1}{\sigma}\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]_{+}^{-(1/\xi+1)}\textrm{\sf exp}\left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}\right\}. \tag{8} \] (try this yourself!)} \frame{ \frametitle{Typical application -- 2: Parameter estimation} Use the GEV probability density function in (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. } \sol \frame{ \frametitle{Typical application -- 2: Parameter estimation} The likelihood is given by \pause \begin{eqnarray*} L(\mu,\sigma,\xi;\boldsymbol{x}) &=& \prod_{i=1}^{m}\frac{1}{\sigma}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+}^{-(1/\xi+1)}\\ & & \hspace{1cm} \times \textrm{\sf exp}\left\{-\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}\right\}\end{eqnarray*} } \frame{ \frametitle{Typical application -- 2: Parameter estimation} The log--likelihood is given by \pause \begin{eqnarray*} \ell(\mu,\sigma,\xi; \boldsymbol{x}) &=& \sum_{i=1}^{m}\textsf{log }\sigma^{-1} + \sum_{i=1}^{m}\textsf{log}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+}^{-(1/\xi+1)} \\ & & \hspace{1cm} - \sum_{i=1}^{m}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+}^{-1/\xi} \pause \\ & & \\ &=& -m\textsf{log }\sigma -(1/\xi +1)\sum_{i=1}^{m}\textsf{log }\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+} \\ & & \hspace{1cm} -\sum_{i=1}^{m}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}. \end{eqnarray*} } \frame{ \frametitle{Typical application -- 2: Parameter estimation} The (log)--likelihood equations are $\frac{\partial \ell}{\partial \mu} = \frac{\partial \ell}{\partial \sigma} = \frac{\partial \ell}{\partial \xi}=0$. For example: \pause \begin{eqnarray*} \frac{\partial \ell}{\partial \mu} &=& \frac{\xi+1}{\sigma}\sum_{i=1}^{m}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]^{-1} \\ & & \hspace{1cm} + \frac{1}{\sigma}\sum_{i=1}^{m}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]^{-(1/\xi+1)} = 0 \\ \pause \end{eqnarray*} The other two (log)--likelihood equations would be found in exactly the same way. } \normal \frame{ \frametitle{Typical application -- 2: Parameter estimation} How would you use the (log) likelihood equations in order to obtain maximum likelihood estimates of $\mu$, $\sigma$ and $\xi$? \vsp Why can't we obtain closed form solutions for $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\xi}$? \vsp How can we get around this? } \sol \frame{ \frametitle{Typical application -- 2: Parameter estimation} The next step would be to replace $\mu$, $\sigma$ and $\xi$ with their corresponding estimators $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\xi}$ (respectively), and then solve for $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\xi}$. \pause \vsp In fact, for the GEV there are no closed--form solutions for $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\xi}$ -- the (log) likelihood equations cannot be solved analytically. \pause \vsp We can get around this by adopting a numerical method to obtain (approximate) solutions to the (log) likelihood equations -- \textsf{R} uses a \textcolor{purple}{\bf Newton--Raphson type algorithm}. } \normal \frame{ \frametitle{Typical application -- 3: Model adequacy} \pause As with all statistical models, there are various \textcolor{purple}{\bf goodness--of--fit properties} that should be considered to check the overall adequacy of the fitted GEV. \pause \vsp These include \textcolor{purple}{\bf probability plots}, \textcolor{purple}{\bf quantile--quantile} ($Q$-$Q$ plots) and simply plotting a histogram of the data against the fitted density. \pause \vsp Again, these will be reviewed shortly via a real--life data demonstration in \textbf{\sf R}. } \frame{ \frametitle{Typical application -- 4: Return level estimation} 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 \pause \vspl \begin{itemize} \item The height of a sea wall to protect against the \textcolor{purple}{\bf once in a hundred year sea--surge}; \vspl \item The ``\textcolor{purple}{\bf fifty year wind speed}'' to provide new structures enough protection against wind damage. \end{itemize} } \frame{ \frametitle{Typical application -- 4: Return level estimation} Such quantities, in extreme value terminology, are usually referred to as \textcolor{purple}{\bf return levels}. \pause \vsp If 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 equal to $1-1/r$ and solving for $x=\hat{z}_{r}$ (provided we have annual maxima). } \frame{ \frametitle{Typical application -- 4: Return level estimation} \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). \pause \vsp 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.\pause \vsp Then we can write down the following probability statement: \begin{eqnarray*} \textrm{\sf Pr}(\textrm{\sf annual maximum}>z_{100}) &=& \frac{1}{100}, \end{eqnarray*} i.e. \[ 1 - \textrm{\sf Pr}(\textrm{\sf annual maximum}\leq z_{100}) =\frac{1}{100}. \tag{9} \] } \frame{ \frametitle{Typical application -- 4: Return level estimation} Now the left--hand--side of (9), in terms of our fitted GEV, is \[ 1-G(\hat{z}_{100};\hat{\mu},\hat{\sigma},\hat{\xi}), \] \pause giving \begin{eqnarray*} 1-\textrm{\sf 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{\sf i.e.} \\ \pause & & \\ \textrm{\sf 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*} } \frame{ \frametitle{Typical application -- 4: Return level estimation} Solving for $\hat{z}_{100}$ gives an estimate of the \textcolor{purple}{\bf 100--year return level} as\pause \[ \hat{z}_{100} = \hat{\mu}+\frac{\hat{\sigma}}{\hat{\xi}}\left[\left(-\textrm{\sf log}(0.99)\right)^{-\hat{\xi}}-1\right]; \] more generally, estimates of the $r$--year return level $z_{r}$ are given by \[ \hat{z}_{r} =\hat{\mu}+\frac{\hat{\sigma}}{\hat{\xi}}\left[\left(-\textrm{\sf log}\left(1-r^{-1}\right)\right)^{-\hat{\xi}}-1\right]. \tag{10} \] } \frame{ \frametitle{Typical application -- 4: Return level estimation} 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 (7)). \pause \vsp Of course, a point estimate of the return level alone is not enough; any estimate requires a corresponding standard error. \pause \vsp To obtain standard errors for return levels, we can use the \textcolor{purple}{\bf delta method}. \pause \vsp However, as we shall see later, we do not recommend using these standard errors to form confidence intervals for $z_{r}$. } \frame{ \frametitle{Sea surge at Wassaw Island} Recall the example from Section 1.3: \pause \vsp \begin{itemize} \item The southeastern coast of the USA is an area often hit by hurricanes \pause \vspl \item The historic city of Savannah, Georgia, has suffered direct hits from 22 Hurricanes since 1871 \pause \vspl \item The city can expect to `brushed' or directly hit once every other year \pause \vspl \item Predictions show that the city is three years overdue it's next direct hit \end{itemize} %\begin{table}[!h] {\small %\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}} } \frame{ \frametitle{Sea surge at Wassaw Island} Recall that, in Chapter 1, we tried to estimated quantities beyond the range of our observed data. \pause \vsp For example, using simple empirical arguments, \[ \textrm{\sf Pr}(\textrm{\sf Sea surge exceeds 14 feet}) = 0, \] simply because, over the time-frame we have data for, we have not observed such an extreme event. \vsp Does this really mean this event is \textcolor{purple}{\bf impossible}? } \frame{ \frametitle{1. Data pre--processing} None necessary --- 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). \pause \vsp Figure 6 shows a time series plot and histogram of the 50 annual maxima: \pause \vspl \begin{itemize} \item No obvious trend in our dataset \pause \vspl \item Issue of dependence? \end{itemize} } \frame{ \frametitle{1. Data pre--processing} \textbf{In \textsf{R}:} \vsp {\color{red} \texttt{> wassaw=scan('wassaw.txt')}} \pause \vsp Or: \vsp {\color{red} \texttt{> wassaw=c(8.5, 8.9, ..., 9.4)} } } \frame{ \frametitle{1. Data pre--processing} \textbf{In \textsf{R}:} \vsp {\color{red} \texttt{year = seq(1975, 2024, 1)}} \hspace{1cm}{\small\textcolor{orange}{$\#$ Sequence of values from 1975 to 2024}} \pause \vspace{0.5cm} Then, to produce the plots shown in Figure 6: \pause \vsp {\color{red} \texttt{> par(mfrow=c(1,2))}} \hspace{1cm} {\small \textcolor{orange}{$\#$ Partitions the plotting space}} \pause \vspl {\color{red} \texttt{> plot(wassaw\textrm{$\sim$}year,type='b')} \vspl \texttt{> hist(wassaw)}} } \frame{ \frametitle{1. Data pre--processing} %\begin{figure}[!h] \begin{center} \psfrag{0}[][]{\sf \tiny 0} \psfrag{2}[][]{\sf \tiny 2} \psfrag{4}[][]{\sf \tiny 4} \psfrag{6}[][]{\sf \tiny 6} \psfrag{1960}[][]{\sf \tiny 1980} \psfrag{1970}[][]{\sf \tiny 1990} \psfrag{1980}[][]{\sf \tiny 2000} \psfrag{1990}[][]{\sf \tiny 2010} \psfrag{2000}[][]{\sf \tiny 2020} \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=11cm,height=5cm]{was1.ps} %\caption{Time series plot (left) and histogram (right) of annual maximum sea--surges observed at Wassaw Island.} %\end{figure} \end{center}} \frame{ \frametitle{2. Fitting the GEV: Parameter estimation} We now use \textsf{R} to maximise \begin{eqnarray*} \ell(\mu,\sigma,\xi; \boldsymbol{x}) &=& -50\textrm{\sf log }\sigma - (1+1/\xi)\sum_{i=1}^{50}\textrm{\sf log}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+} \\ &&\hspace{1cm} - \sum_{i=1}^{50}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+}^{-1/\xi} \end{eqnarray*} with respect to $\mu$, $\sigma$ and $\xi$. \pause \vsp \textbf{In \textsf{R}:} \vsp {\color{red} \texttt{> theta=c(mean(wassaw), sd(wassaw), 0.1)}} \hspace{1cm} {\small\textcolor{orange}{$\#$ Sets up parameter vector $\theta=(\mu,\sigma,\xi)$}} } \frame{ \frametitle{2. Fitting the GEV: Parameter estimation} {\small {\color{red} \textcolor{purple}{\texttt{\bf gev.loglik=function(theta)$\boldsymbol{\{}$}} \\ \hspace{0.4cm} \texttt{mu=theta[1]}\\ \hspace{0.4cm} \texttt{sigma=theta[2]}\\ \hspace{0.4cm} \texttt{xi=theta[3]}\\ \hspace{0.4cm} \texttt{m=min((1+(xi*(dataset-mu)/sigma))) }\\ \hspace{0.4cm} \texttt{if(m<0.00001)return(as.double(1000000))} \\ \hspace{0.4cm} \texttt{if(sigma<0.00001)return(as.double(1000000))}\\ \hspace{0.4cm} \texttt{if(xi==0)$\{$}\\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)} \\ \texttt{+ \hspace{0.6cm}-sum((dataset-mu)/sigma)}\\ \texttt{+ \hspace{0.6cm}-sum(exp(-((dataset-mu)/sigma)))$\}$}\\ \hspace{0.4cm}\texttt{else$\{$} \\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)}\\ \texttt{+ \hspace{0.6cm}-(1/xi+1)*sum(log(1+(xi*(dataset-mu)/sigma)))}\\ \texttt{+ \hspace{0.6cm}-sum((1+(xi*(dataset-mu)/sigma))**(-1/xi))$\}$} \\ \hspace{0.8cm} \texttt{return(-loglik)$\}$} }} } \frame{ \frametitle{2. Fitting the GEV: Parameter estimation} {\small {\color{red} \texttt{ > gev.loglik=function(theta)$\{$} \\ \hspace{0.4cm} \textcolor{purple}{\texttt{\bf mu=theta[1]}}\\ \hspace{0.4cm} \textcolor{purple}{\texttt{\bf sigma=theta[2]}}\\ \hspace{0.4cm} \textcolor{purple}{\texttt{\bf xi=theta[3]}}\\ \hspace{0.4cm} \texttt{m=min((1+(xi*(dataset-mu)/sigma))) }\\ \hspace{0.4cm} \texttt{if(m<0.00001)return(as.double(1000000))} \\ \hspace{0.4cm} \texttt{if(sigma<0.00001)return(as.double(1000000))}\\ \hspace{0.4cm} \texttt{if(xi==0)$\{$}\\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)} \\ \texttt{+ \hspace{0.6cm}-sum((dataset-mu)/sigma)}\\ \texttt{+ \hspace{0.6cm}-sum(exp(-((dataset-mu)/sigma)))$\}$}\\ \hspace{0.4cm}\texttt{else$\{$} \\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)}\\ \texttt{+ \hspace{0.6cm}-(1/xi+1)*sum(log(1+(xi*(dataset-mu)/sigma)))}\\ \texttt{+ \hspace{0.6cm}-sum((1+(xi*(dataset-mu)/sigma))**(-1/xi))$\}$} \\ \hspace{0.8cm} \texttt{return(-loglik)$\}$} }} } \frame{ \frametitle{2. Fitting the GEV: Parameter estimation} {\small {\color{red} \texttt{ > gev.loglik=function(theta)$\{$} \\ \hspace{0.4cm} \texttt{mu=theta[1]}\\ \hspace{0.4cm} \texttt{sigma=theta[2]}\\ \hspace{0.4cm} \texttt{xi=theta[3]}\\ \hspace{0.4cm} \textcolor{purple}{\texttt{\bf m=min((1+(xi*(dataset--mu)/sigma))) }}\\ \hspace{0.4cm} \textcolor{purple}{\texttt{\bf if(m$<$0.00001)return(as.double(1000000))}} \\ \hspace{0.4cm} \texttt{if(sigma<0.00001)return(as.double(1000000))}\\ \hspace{0.4cm} \texttt{if(xi==0)$\{$}\\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)} \\ \texttt{+ \hspace{0.6cm}-sum((dataset-mu)/sigma)}\\ \texttt{+ \hspace{0.6cm}-sum(exp(-((dataset-mu)/sigma)))$\}$}\\ \hspace{0.4cm}\texttt{else$\{$} \\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)}\\ \texttt{+ \hspace{0.6cm}-(1/xi+1)*sum(log(1+(xi*(dataset-mu)/sigma)))}\\ \texttt{+ \hspace{0.6cm}-sum((1+(xi*(dataset-mu)/sigma))**(-1/xi))$\}$} \\ \hspace{0.8cm} \texttt{return(-loglik)$\}$} }} } \frame{ \frametitle{2. Fitting the GEV: Parameter estimation} {\small {\color{red} \texttt{ > gev.loglik=function(theta)$\{$} \\ \hspace{0.4cm} \texttt{mu=theta[1]}\\ \hspace{0.4cm} \texttt{sigma=theta[2]}\\ \hspace{0.4cm} \texttt{xi=theta[3]}\\ \hspace{0.4cm} \texttt{m=min((1+(xi*(dataset-mu)/sigma)))}\\ \hspace{0.4cm} \texttt{if(m<0.00001)return(as.double(1000000))} \\ \hspace{0.4cm} \textcolor{purple}{\texttt{\bf if(sigma$<$0.00001)return(as.double(1000000))}}\\ \hspace{0.4cm} \texttt{if(xi==0)$\{$}\\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)} \\ \texttt{+ \hspace{0.6cm}-sum((dataset-mu)/sigma)}\\ \texttt{+ \hspace{0.6cm}-sum(exp(-((dataset-mu)/sigma)))$\}$}\\ \hspace{0.4cm}\texttt{else$\{$} \\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)}\\ \texttt{+ \hspace{0.6cm}-(1/xi+1)*sum(log(1+(xi*(dataset-mu)/sigma)))}\\ \texttt{+ \hspace{0.6cm}-sum((1+(xi*(dataset-mu)/sigma))**(-1/xi))$\}$} \\ \hspace{0.8cm} \texttt{return(-loglik)$\}$}}}} \frame{ \frametitle{2. Fitting the GEV: Parameter estimation} {\small {\color{red} \texttt{ > gev.loglik=function(theta)$\{$} \\ \hspace{0.4cm} \texttt{mu=theta[1]}\\ \hspace{0.4cm} \texttt{sigma=theta[2]}\\ \hspace{0.4cm} \texttt{xi=theta[3]}\\ \hspace{0.4cm} \texttt{m=min((1+(xi*(dataset-mu)/sigma))) }\\ \hspace{0.4cm} \texttt{if(m<0.00001)return(as.double(1000000))} \\ \hspace{0.4cm} \texttt{if(sigma<0.00001)return(as.double(1000000))}\\ \hspace{0.4cm} \textcolor{purple}{\texttt{\bf if(xi==0)$\boldsymbol{\{}$}}\\ \hspace{0.6cm} \textcolor{blue}{\textsf{loglik=--length(dataset)*log(sigma)}} \\ \texttt{+ \hspace{0.6cm}\textcolor{blue}{\sf --sum((dataset--mu)/sigma)}}\\ \texttt{+ \hspace{0.6cm}\textcolor{blue}{\sf --sum(exp(--((dataset--mu)/sigma)))}\textcolor{purple}{\bf $\boldsymbol{\}}$}}\\ \hspace{0.4cm}\texttt{else$\{$} \\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)}\\ \texttt{+ \hspace{0.6cm}-(1/xi+1)*sum(log(1+(xi*(dataset-mu)/sigma)))}\\ \texttt{+ \hspace{0.6cm}-sum((1+(xi*(dataset-mu)/sigma))**(-1/xi))$\}$} \\ \hspace{0.8cm} \texttt{return(-loglik)$\}$} }}} \frame{ \frametitle{2. Fitting the GEV: Parameter estimation} {\small {\color{red} \texttt{ > gev.loglik=function(theta)$\{$} \\ \hspace{0.4cm} \texttt{mu=theta[1]}\\ \hspace{0.4cm} \texttt{sigma=theta[2]}\\ \hspace{0.4cm} \texttt{xi=theta[3]}\\ \hspace{0.4cm} \texttt{m=min((1+(xi*(dataset-mu)/sigma))) }\\ \hspace{0.4cm} \texttt{if(m<0.00001)return(as.double(1000000))} \\ \hspace{0.4cm} \texttt{if(sigma<0.00001)return(as.double(1000000))}\\ \hspace{0.4cm} \texttt{if(xi==0)$\{$}\\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)} \\ \texttt{+ \hspace{0.6cm}\texttt{-sum((dataset-mu)/sigma)}}\\ \texttt{+ \hspace{0.6cm}\texttt{-sum(exp(-((dataset-mu)/sigma)))$\}$}}\\ \hspace{0.4cm}\textcolor{purple}{\texttt{\bf else$\boldsymbol{\{}$}} \\ \hspace{0.6cm} \textcolor{blue}{\textsf{loglik=--length(dataset)*log(sigma)}}\\ \texttt{+ \hspace{0.6cm}\textcolor{blue}{\sf --(1/xi+1)*sum(log(1+(xi*(dataset--mu)/sigma)))}}\\ \texttt{+ \hspace{0.6cm}\textcolor{blue}{\sf --sum((1+(xi*(dataset--mu)/sigma))**(--1/xi))}\textcolor{purple}{\bf $\boldsymbol{\}}$}} \\ \hspace{0.8cm} \texttt{return(-loglik)$\}$} }}} \frame{ \frametitle{2. Fitting the GEV: Parameter estimation} {\small {\color{red} \textcolor{purple}{\texttt{\bf gev.loglik=function(theta)$\boldsymbol{\{}$}} \\ \hspace{0.4cm} \texttt{mu=theta[1]}\\ \hspace{0.4cm} \texttt{sigma=theta[2]}\\ \hspace{0.4cm} \texttt{xi=theta[3]}\\ \hspace{0.4cm} \texttt{m=min((1+(xi*(dataset-mu)/sigma))) }\\ \hspace{0.4cm} \texttt{if(m<0.00001)return(as.double(1000000))} \\ \hspace{0.4cm} \texttt{if(sigma<0.00001)return(as.double(1000000))}\\ \hspace{0.4cm} \texttt{if(xi==0)$\{$}\\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)} \\ \texttt{+ \hspace{0.6cm}-sum((dataset-mu)/sigma)}\\ \texttt{+ \hspace{0.6cm}-sum(exp(-((dataset-mu)/sigma)))$\}$}\\ \hspace{0.4cm}\texttt{else$\{$} \\ \hspace{0.6cm} \texttt{loglik=-length(dataset)*log(sigma)}\\ \texttt{+ \hspace{0.6cm}-(1/xi+1)*sum(log(1+(xi*(dataset-mu)/sigma)))}\\ \texttt{+ \hspace{0.6cm}-sum((1+(xi*(dataset-mu)/sigma))**(-1/xi))$\}$} \\ \hspace{0.8cm} \textcolor{purple}{\texttt{\bf return(--loglik)$\boldsymbol{\}}$}} }}} \frame{ \frametitle{2. Fitting the GEV: Parameter estimation} \textbf{In \textsf{R}:} \vsp {\color{red} \texttt{> dataset=wassaw}} \hspace{1cm}{\small\textcolor{orange}{\small $\#$ Attaches the sea--surge extremes to ``dataset''}} \pause \vsp {\color{red} \texttt{> nlm(gev.loglik,theta)} \hspace{0.5cm} \hspace{1cm}{\small \textcolor{orange}{\small $\#$ Minimisation routine}} \pause \vspl \texttt{\$minimum}\\ \texttt{[1] 89.52412} \pause \vspl \texttt{\$estimate}\\ \texttt{[1] 8.7112735 1.3114836 -0.1084451} \pause \vspl \texttt{\$gradient}\\ \texttt{[1] 3.350727e-06 2.316675e-05 2.145839e-06} \vspl \texttt{\$code}\\ \texttt{[1] 1} }} \frame{ \frametitle{2. Fitting the GEV: Parameter estimation} \textbf{In \textsf{R}:} \vsp {\color{red} \texttt{> dataset=wassaw}} \hspace{1cm}{\small\textcolor{orange}{\small $\#$ Attaches the sea--surge extremes to ``dataset''}} \vsp {\color{red} \texttt{> nlm(gev.loglik,theta)} \hspace{0.5cm} \hspace{1cm}{\small \textcolor{orange}{\small $\#$ Minimisation routine}} \vspl \texttt{\$minimum}\\ \texttt{[1] 89.52412} \vspl \texttt{\$estimate}\\ \texttt{[1] \textbf{8.7112735} \textbf{1.3114836} \textbf{-0.1084451}} \vspl \texttt{\$gradient}\\ \texttt{[1] 3.350727e-06 2.316675e-05 2.145839e-06} \vspl \texttt{\$code}\\ \texttt{[1] 1} }} \frame{ \frametitle{2. Fitting the GEV: Standard errors} \textbf{In \textsf{R}:} \vsp {\color{red} \texttt{> A=nlm(gev.loglik,theta,hessian=TRUE)}} \pause\\ \hspace{1cm}\textcolor{orange}{\small $\#$ Stores output, including Hessian, in \texttt{A}} \pause \vsp Now: \vspl {\color{red} \texttt{> A} \vspl \hspace{1cm}$\vdots$ \vspl \texttt{\$hessian} \begin{tabular}{r r r r} &\texttt{[,1]} &\texttt{[,2]} &\texttt{[,3]}\\ \texttt{[1,]} &\texttt{27.375570} &\texttt{-4.837481} &\texttt{18.06140} \\ \texttt{[2,]} &\texttt{-4.837481} &\texttt{56.151740} &\texttt{29.89004}\\ \texttt{[3,]} &\texttt{18.061397} &\texttt{29.890039} &\texttt{118.18086}\\ \end{tabular}}} \frame{ \frametitle{2. Fitting the GEV: Standard errors} The \textcolor{purple}{\bf Hessian} is the matrix of second--order partial derivatives. \pause \vsp In the case of our log--likelihood function, this is \pause \[ \textsf{Hessian} = \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) = -\pause I_{\textrm{\sf obs}}, \] evaluated at the MLEs. \pause \vsp However, since \texttt{gev.loglik} returns the \textit{negative} log--likelihood, in this example \[ \textsf{Hessian} = I_{\textrm{\sf obs}}!! \] } \frame{ \frametitle{2. Fitting the GEV: Standard errors} Recall from MAS2901 that inversion of the Information matrix gives the \textcolor{purple}{\bf Variance--Covariance matrix}.\pause \vsp Thus, we have a $3\times 3$ matrix to invert! \pause \vsp You should be able to do this (e.g. assignments), but let's use the \texttt{solve} command in \textsf{\bf R}. } \frame{ \frametitle{2. Fitting the GEV: Standard errors} In \textbf{\sf R}, \texttt{solve} will solve the system of equations: \pause \[ \boldsymbol{a}X = \boldsymbol{b}. \] \pause In our example, $\boldsymbol{a} = I_{\textsf{obs}}$. Omitting $\boldsymbol{b}$ in the execution of \texttt{solve} assumes the identity matrix for $\boldsymbol{b}$, giving: \pause \[ I_{\textsf{obs}}\cdot X = \left(\begin{array}{c c c} 1 &0 &0 \\ 0 &1 &0 \\ 0 &0 &1 \end{array} \right) \pause \] Thus, solving for $X$ will return the inverse of $I_{\textsf{obs}}$ -- the \textcolor{purple}{\bf variance--covariance matrix}! } \frame{ \frametitle{2. Fitting the GEV: Standard errors} \textbf{In \textsf{R}:} \vsp {\color{red} \texttt{> solve(A\$hessian)} \begin{tabular}{r r r r} &\texttt{[,1]} &\texttt{[,2]} &\texttt{[,3]}\\ \texttt{[1,]} &\texttt{0.043869792} &\texttt{0.008491490} &\texttt{-0.008852201}\\ \texttt{[2,]} &\texttt{0.008491490} &\texttt{0.022223135} &\texttt{-0.006918367}\\ \texttt{[3,]} &\texttt{-0.008852201} &\texttt{-0.006918367} &\texttt{0.011564254}\\ \end{tabular} \pause \vsp \texttt{> varcovar=solve(A\$hessian)} \pause \vsp \texttt{> sqrt(diag(varcovar))} \pause \vspl \texttt{[1] 0.2094512 0.1490743 0.1075372}}} \frame{ \frametitle{2. Fitting the GEV: Standard errors} Thus, we now have the following inference for our annual maximum sea surges, in terms of the GEV distribution: \pause \begin{displaymath} \hat{\mu}=8.711 \textrm{\sf (0.209)}\qquad \hat{\sigma}=1.311 \textrm{\sf (0.149)} \qquad \hat{\xi}=-0.108 \textrm{\sf (0.108)} \end{displaymath} \pause \noindent From this, we can construct \textcolor{purple}{\bf confidence intervals} in the usual way (``Wald'' intervals). \pause \vsp For example: $\textrm{\sf parameter estimate} \pm 1.96 \times \textrm{\sf standard error}$ for a 95\% CI, owing to the normality of maximum likelihood estimators --- giving \pause \begin{displaymath} (8.301, 9.121) \qquad (1.019, 1.603)\qquad \textrm{\sf and} \qquad (-0.320, 0.104) \end{displaymath} for $\mu$, $\sigma$ and $\xi$ (respectively). \pause \vsp Note that the confidence interval for $\xi$ \textcolor{purple}{\bf includes zero} --- a Gumbel--type tail for our data could be appropriate. } \frame{ \frametitle{3. Model adequacy: Probability plots} \textcolor{purple}{\bf Probability plots} compare \textit{empirical cumulative probabilities} with corresponding values from the fitted model. \pause \vsp For example, the ordered sea--surges at Wassaw, $x_{(i)}$, $i=1, \ldots, 50$, are: \begin{center} {\small \begin{tabular}{|c c c c c c c c c c|}\hline 6.7 &6.8 &7.3 &7.3 &7.3 &7.4 &7.6 &7.7 &7.8 &8.1 \\ 8.2 &8.3 &8.4 &8.5 &8.6 &8.7 &8.7 &8.8 &8.8 &8.9 \\ 8.9 &8.9 &8.9 &9.0 &9.1 &9.1 &9.3 &9.3 &9.4 &9.5 \\ 9.6 &9.6 &9.7 &9.7 &9.8 &10.0 &10.0 &10.2 &10.3 &10.4 \\ 10.5 &10.5 &10.8 &11.1 &11.5 &11.8 &12.2 &12.6 &12.7 &13.0 \\ \hline \end{tabular} \pause \vsp \begin{tabular}{|c|c|c|} \hline $x_{(i)}$ &Fitted: $G(x_{(i)};\hat{\mu},\hat{\sigma},\hat{\xi})$ &Empirical: $i/(n+1)$ \\ \hline 6.7 &0.016 &0.020\\ \pause 6.8 &0.021 &0.039 \\ \pause 7.3 &0.063 &0.059 \\ $\vdots$ &$\vdots$ &$\vdots$ \\ \hline \end{tabular}} \end{center} } \frame{ \frametitle{3. Model adequacy: Probability plots} \textcolor{purple}{\bf Probability plots} compare \textit{empirical cumulative probabilities} with corresponding values from the fitted model. \vsp For example, the ordered sea--surges at Wassaw, $x_{(i)}$, $i=1, \ldots, 50$, are: \begin{center} {\small \begin{tabular}{|c c c c c c c c c c|}\hline 6.7 &6.8 &7.3 &7.3 &7.3 &7.4 &7.6 &7.7 &7.8 &8.1 \\ 8.2 &8.3 &8.4 &8.5 &8.6 &8.7 &8.7 &8.8 &8.8 &8.9 \\ 8.9 &8.9 &8.9 &9.0 &9.1 &9.1 &9.3 &9.3 &9.4 &9.5 \\ 9.6 &9.6 &9.7 &9.7 &9.8 &10.0 &10.0 &10.2 &10.3 &10.4 \\ 10.5 &10.5 &10.8 &11.1 &11.5 &11.8 &12.2 &12.6 &12.7 &13.0 \\ \hline \end{tabular} \vsp \begin{tabular}{|c|c|c|} \hline $x_{(i)}$ &Fitted: $G(x_{(i)};\hat{\mu},\hat{\sigma},\hat{\xi})$ &Empirical: $i/(n+1)$ \\ \hline 6.7 &\textbf{0.016} &\textbf{0.020}\\ 6.8 &\textbf{0.021} &\textbf{0.039} \\ 7.3 &\textbf{0.063} &\textbf{0.059} \\ $\vdots$ &$\vdots$ &$\vdots$ \\ \hline \end{tabular}} \end{center} } \frame{ \frametitle{3. Model adequacy: Quantile plots} \textcolor{purple}{\bf Quantile plots} compare empirical quantiles with corresponding values from the fitted model: \pause \vsp \begin{center} {\small \begin{tabular}{|c|c|c|} \hline Empirical: $x_{(i)}$ &Prob: $i/(n+1)$ &Fitted: $G^{-1}(i/(n+1);\hat{\mu},\hat{\sigma},\hat{\xi})$ \\ \hline 6.7 &0.020 &6.78 \\ \pause 6.8 &0.039 &7.07 \\ \pause 7.3 &0.059 &7.27 \\ $\vdots$ &$\vdots$ &$\vdots$ \\ \hline \end{tabular} } \end{center} } \frame{ \frametitle{3. Model adequacy: Quantile plots} \textcolor{purple}{\bf Quantile plots} compare empirical quantiles with corresponding values from the fitted model: \vsp \begin{center} {\small \begin{tabular}{|c|c|c|} \hline Empirical: $x_{(i)}$ &Prob: $i/(n+1)$ &Fitted: $G^{-1}(i/(n+1);\hat{\mu},\hat{\sigma},\hat{\xi})$ \\ \hline \textbf{6.7} &0.020 &\textbf{6.78} \\ \textbf{6.8} &0.039 &\textbf{7.07} \\ \textbf{7.3} &0.059 &\textbf{7.27} \\ $\vdots$ &$\vdots$ &$\vdots$ \\ \hline \end{tabular} } \end{center} } \frame{ \frametitle{Probability plot in \textsf{R}} \textbf{In \textsf{R}:} \pause \vsp {\color{red} \texttt{ >ordered=sort(dataset)}} \\ \hspace{1cm} \textcolor{orange}{\small $\#$ Orders the data} \pause \vsp {\color{red} \texttt{> empirical=vector('numeric',length(ordered))} \pause \vspl \texttt{> for(i in 1:length(empirical))$\{$} \\ \hspace{0.4cm} \texttt{empirical[i]=i/(length(dataset)+1)$\}$}} } \frame{ \frametitle{Probability plot in \textsf{R}} \noindent The function \texttt{GEV.DF} defines the distribution function for the GEV: \pause \vsp {\color{red} \texttt{GEV.DF=function(data,mu,sigma,xi)$\{$} \\ \hspace{0.4cm} \texttt{if(xi==0)$\{$} \\ \hspace{0.45cm} \texttt{GEV=exp(-exp(-((data-mu)/sigma)))$\}$}\\ \hspace{0.4cm} \texttt{else$\{$} \\ \hspace{0.45cm} \texttt{GEV=exp(-(1+xi*((data-mu)/sigma))**(-1/xi))$\}$}\\ \hspace{0.4cm}\texttt{return(GEV)$\}$}} } \frame{ \frametitle{Probability plot in \textsf{R}} \noindent Then: \pause \vsp {\scriptsize {\color{red} \texttt{> model=vector('numeric',length(dataset))} \\ \texttt{> for(i in 1:length(model))$\{$}\\ \hspace{0.4cm} \texttt{model[i]=GEV.DF(ordered[i],A\$est[1],A\$est[2],A\$est[3])$\}$}} }\pause \vsp Plotting \texttt{model} against \texttt{empirical} produces the corresponding probability plot. } \frame{ \frametitle{Probability/Quantile plots in \textsf{R}} %\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=11cm,height=6cm]{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} } \frame{ \frametitle{4. Return level estimation} In question 2 of Section 1.3, 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 \vspl \begin{itemize} \item [(i)] once in ten years; \vspl \item [(ii)] once in a hundred years. \end{itemize} } \frame{ \frametitle{4. Return level estimation} Using the data alone we could not obtain a meaningful estimate of (ii) because we only have 50 years of data. \pause \vsp However, we can now use our fitted GEV to extrapolate beyond the range of our data to estimate such return levels. \pause \vsp Using Equation (10): \[ \hat{z}_{r} =\hat{\mu}+\frac{\hat{\sigma}}{\hat{\xi}}\left[\left(-\textrm{\sf log}\left(1-r^{-1}\right)\right)^{-\hat{\xi}}-1\right], \tag{10} \] we find that \[ \hat{z}_{10} = 8.711-\frac{1.311}{0.108}\left[\left(-\textrm{\sf log}\left(1-10^{-1}\right)\right)^{0.108}-1\right] = 11.33 \textrm{ \sf feet}; \pause \] similarly, \[ \hat{z}_{10} = 8.711-\frac{1.311}{0.108}\left[\left(-\textrm{\sf log}\left(1-100^{-1}\right)\right)^{0.108}-1\right] = 13.46 \textrm{ \sf feet}. \] } \frame{ \frametitle{4. Return level estimation} 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). \pause \vsp In fact, due to the \textcolor{purple}{\bf invariance property} of maximum likelihood estimators, our estimates of $z_{10}$ and $z_{100}$ are also the maximum likelihood estimators of these quantities. } \frame{ \frametitle{4. Return level estimation} \noindent As with our inference for the GEV parameters, it is preferable to quote estimates of return levels with their estimated standard errors. \pause \vsp Since $z_{r}$ is a function of the GEV parameters $\mu$, $\sigma$ and $\xi$, we can use the \textcolor{purple}{\bf delta method} to find the approximate variance of $\hat{z}_{r}$. \pause \vsp Specifically, \begin{eqnarray*} \textrm{\sf 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 \pause \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{\sf log }y_{r}\right], \end{eqnarray*} where $y_{r}=-\textrm{\sf log}(1-r^{-1})$, evaluated at $(\hat{\mu},\hat{\sigma},\hat{\xi})$. } \frame{ \frametitle{Return level standard errors in \textsf{R}} \noindent Recall that we previously stored $V$ in the matrix \texttt{varcovar} in \textsf{R}. \pause \vsp Also, recall that our estimates of the GEV parameters are stored in \texttt{A\$est}. \pause \vsp We can define $\nabla z_{10}$, for example, in \textsf{R} as: \pause \vsp {\scriptsize {\color{red} \texttt{>}\textcolor{purple}{\bf y10=--log(1--(1/10))}\\ \texttt{>del=matrix(ncol=1,nrow=3)}\\ \texttt{>del[1,1]=1}\\ \texttt{>del[2,1]=-((A\$est[3])**(-1))*(1-(y10**(-A\$est[3])))}\\ \texttt{>del[3,1]=((A\$est[2])*((A\$est[3])**(-2))*(1-((y10)**(-A\$est[3]))))}\\ \texttt{-((A\$est[2])*((A\$est[3])**(-1))*((y10)**(-(A\$est[3])))*log(y10))}\\ \texttt{>del.transpose=t(del)}}} } \frame{ \frametitle{Return level standard errors in \textsf{R}} \noindent Recall that we previously stored $V$ in the matrix \texttt{varcovar} in \textsf{R}. \vsp Also, recall that our estimates of the GEV parameters are stored in \texttt{A\$est}. \vsp We can define $\nabla z_{10}$, for example, in \textsf{R} as \vsp {\scriptsize {\color{red} \texttt{>y10=-log(1-(1/10))}\\ \texttt{>}\textsf{\textcolor{purple}{\bf del=matrix(ncol=1,nrow=3)}}\\ \texttt{>del[1,1]=1}\\ \texttt{>del[2,1]=-((A\$est[3])**(-1))*(1-(y10**(-A\$est[3])))}\\ \texttt{>del[3,1]=((A\$est[2])*((A\$est[3])**(-2))*(1-((y10)**(-A\$est[3]))))}\\ \texttt{-((A\$est[2])*((A\$est[3])**(-1))*((y10)**(-(A\$est[3])))*log(y10))}\\ \texttt{>del.transpose=t(del)}}} } \frame{ \frametitle{Return level standard errors in \textsf{R}} \noindent Recall that we previously stored $V$ in the matrix \texttt{varcovar} in \textsf{R}. \vsp Also, recall that our estimates of the GEV parameters are stored in \texttt{A\$est}. \vsp We can define $\nabla z_{10}$, for example, in \textsf{R} as \vsp {\scriptsize {\color{red} \texttt{>y10=-log(1-(1/10))}\\ \texttt{>}\textsf{\textcolor{purple}{\bf del=matrix(ncol=1,nrow=3)}}\\ \texttt{>}\textsf{\textcolor{purple}{\bf del[1,1]=1}}\\ \texttt{>del[2,1]=-((A\$est[3])**(-1))*(1-(y10**(-A\$est[3])))}\\ \texttt{>del[3,1]=((A\$est[2])*((A\$est[3])**(-2))*(1-((y10)**(-A\$est[3]))))}\\ \texttt{-((A\$est[2])*((A\$est[3])**(-1))*((y10)**(-(A\$est[3])))*log(y10))}\\ \texttt{>del.transpose=t(del)}}} } \frame{ \frametitle{Return level standard errors in \textsf{R}} \noindent Recall that we previously stored $V$ in the matrix \texttt{varcovar} in \textsf{R}. \vsp Also, recall that our estimates of the GEV parameters are stored in \texttt{A\$est}. \vsp We can define $\nabla z_{10}$, for example, in \textsf{R} as \vsp {\scriptsize {\color{red} \texttt{>y10=-log(1-(1/10))}\\ \texttt{>}\textsf{\textcolor{purple}{\bf del=matrix(ncol=1,nrow=3)}}\\ \texttt{>del[1,1]=1}\\ \texttt{>}\textsf{\textcolor{purple}{\bf del[2,1]=-((A\$est[3])**(--1))*(1--(y10**(--A\$est[3])))}}\\ \texttt{>del[3,1]=((A\$est[2])*((A\$est[3])**(-2))*(1-((y10)**(-A\$est[3]))))}\\ \texttt{-((A\$est[2])*((A\$est[3])**(-1))*((y10)**(-(A\$est[3])))*log(y10))}\\ \texttt{>del.transpose=t(del)}}} } \frame{ \frametitle{Return level standard errors in \textsf{R}} \noindent Recall that we previously stored $V$ in the matrix \texttt{varcovar} in \textsf{R}. \vsp Also, recall that our estimates of the GEV parameters are stored in \texttt{A\$est}. \vsp We can define $\nabla z_{10}$, for example, in \textsf{R} as \vsp {\scriptsize {\color{red} \texttt{>y10=-log(1-(1/10))}\\ \texttt{>}\textsf{\textcolor{purple}{\bf del=matrix(ncol=1,nrow=3)}}\\ \texttt{>del[1,1]=1}\\ \texttt{>del[2,1]=-((A\$est[3])**(-1))*(1-(y10**(-A\$est[3])))}\\ \texttt{>}\textsf{\textcolor{purple}{\bf del[3,1]=((A\$est[2])*((A\$est[3])**(--2))*(1--((y10)**(--A\$est[3]))))}}\\ \texttt{\textsf{\textcolor{purple}{\bf --((A\$est[2])*((A\$est[3])**(--1))*((y10)**(--(A\$est[3])))*log(y10))}}}\\ \texttt{>del.transpose=t(del)}}} } \frame{ \frametitle{Return level standard errors in \textsf{R}} \noindent Recall that we previously stored $V$ in the matrix \texttt{varcovar} in \textsf{R}. \vsp Also, recall that our estimates of the GEV parameters are stored in \texttt{A\$est}. \vsp We can define $\nabla z_{10}$, for example, in \textsf{R} as \vsp {\scriptsize {\color{red} \texttt{>y10=-log(1-(1/10))}\\ \texttt{>del=matrix(ncol=1,nrow=3)}\\ \texttt{>del[1,1]=1}\\ \texttt{>del[2,1]=-((A\$est[3])**(-1))*(1-(y10**(-A\$est[3])))}\\ \texttt{>del[3,1]=((A\$est[2])*((A\$est[3])**(-2))*(1-((y10)**(-A\$est[3]))))}\\ \texttt{-((A\$est[2])*((A\$est[3])**(-1))*((y10)**(-(A\$est[3])))*log(y10))}\\ \texttt{>}\textsf{\textcolor{purple}{\bf del.transpose=t(del)}}}} } \frame{ \frametitle{Return level standard errors in \textsf{R}} \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: \pause \vsp {\color{red} \texttt{> sqrt(del.transpose\%*\%varcovar\%*\%del)} \pause \vspl \begin{tabular}{r r} &\texttt{[,1]}\\ \texttt{[1,]} &\texttt{0.3614568} \pause \end{tabular}} \vsp \noindent Estimated standard errors for other return levels can be obtained in a similar way. } \frame{ \frametitle{Return level standard errors in \textsf{R}} For example, for the standard error for the 100--year return level, we would replace \vspl {\color{red} \texttt{> y10=-log(1-(1/10))}} \pause \vsp with \vspl {\color{red} \texttt{> y100=-log(1-(1/100))}} \pause \vsp and then \texttt{y10} would be replaced with \texttt{y100} throughout. } \frame{ \frametitle{4. Return level estimation} Table 3 shows a range of estimated return levels, with associated standard errors in parentheses. \pause \vsp 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 \textcolor{purple}{\bf can be misleading}. \pause %\begin{table}[!h] %\centering {\small \begin{center} \begin{tabular}{|c c c c|} \hline \multicolumn{4}{|c|}{Return level: MLE (s.e.)} \\ $z_{10}$ &$z_{100}$ &$z_{200}$ &$z_{1000}$ \\ \hline 11.33 (0.361)&13.46 (0.938) &13.99 (1.182)&15.09 (1.821) \\ \hline \end{tabular} \end{center}}} \frame{ \frametitle{Using the \textsf{ismev} package in \textsf{R}} \begin{center} \LARGE \textit{Class demonstration in \textsf{R}} \end{center} } \frame{ \frametitle{Example 2.3: Rainfall in the Lake District} The data shown in Table 4 are annual maximum rainfall accumulations, obtained from daily records, for a period of 21 years (2001--2021 inclusive) at Eskdale in the Lake District. \begin{center} %\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} \end{center} } \frame{ \frametitle{Example 2.3: Rainfall in the Lake District} \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). \] } \frame{ \frametitle{Example 2.3: Rainfall in the Lake District} \begin{itemize} \item [(a)] Estimate the standard errors for $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\xi}$, and use these to estimate the 95\% confidence intervals for the GEV parameters. Comment. \item [(b)] What is the estimated correlation between $\sigma$ and $\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} } \sol \frame{ \frametitle{Example 2.3: Solution to part(a) (1/8)} The standard errors are found from the variance--covariance matrix, which is $I_{O}^{-1}$. \pause \vsp \begin{itemize} \item[\textbf{1.}] \textbf{Find the determinant of $\boldsymbol{I_{O}}$:} \pause {\small \begin{eqnarray*} \textsf{det}(I_{O}) &=& 0.0062(0.0091\times48.4122 - 0.2114^{2}) \\ & &+ 0.0046(-0.0046\times48.4122 - 0.1962\times(-0.2114))\\ & &+ 0.1962(0.0046\times0.2114 - 0.\textbf{1}962\times0.0091) \\ &=& 0.00146. \end{eqnarray*} } \end{itemize} } \frame{ \frametitle{Example 2.3: Solution to part(a) (2/8)} \begin{itemize} \item[\textbf{2.}] \textbf{Find the transpose of $\boldsymbol{I_{O}}$:} \[ I_{O}^{T} = \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) \] \end{itemize} } \frame{ \frametitle{Example 2.3: Solution to part(a) (3/8)} \begin{itemize} \item [\textbf{3.}] \textbf{Find the determinants of each of the $\boldsymbol{2\times 2}$ minor matrices of $\boldsymbol{I}_{O}^{T}$:} \begin{eqnarray*} \left|\begin{array}{c c} 0.0091 &-0.2114 \\ -0.2114 &48.4122 \end{array} \right| &=& 0.395861\\ \pause & & \\ \left|\begin{array}{c c} -0.0046 &-0.2114 \\ 0.1962 &48.4122 \end{array} \right| &=& -0.181219 \\ & & \\ \left|\begin{array}{c c} -0.0046 &0.0091 \\ 0.1962 &-0.2114 \end{array} \right| &=& -0.000813 \end{eqnarray*} \end{itemize} } \frame{ \frametitle{Example 2.3: Solution to part(a) (4/8)} \begin{eqnarray*} \left|\begin{array}{c c} -0.0046 &0.1962 \\ -0.2114 &48.4122 \end{array} \right| &=& -0.181219 \pause\\ & & \\ \left|\begin{array}{c c} 0.0062 &0.1962 \\ 0.1962 &48.4122 \end{array} \right| &=& 0.261661 \pause \\ & &\\ \left|\begin{array}{c c} 0.0062 &-0.0046 \\ 0.1962 &-0.2144 \end{array} \right| &=& -0.000408 \end{eqnarray*} } \frame{ \frametitle{Example 2.3: Solution to part(a) (5/8)} \begin{eqnarray*} \left|\begin{array}{c c} -0.0046 &0.1962 \\ 0.0091 &-0.2114 \end{array} \right| &=& -0.000813\\ & &\\ \left|\begin{array}{c c} 0.0062 &0.1962 \\ -0.0046 &-0.2114 \end{array} \right| &=& -0.000408\\ & &\\ \left|\begin{array}{c c} 0.0062 &-0.0046 \\ -0.0046 &0.0091 \end{array} \right| &=& 0.000035 \end{eqnarray*} } \frame{ \frametitle{Example 2.3: Solution to part(a) (6/8)} \begin{itemize} \item [\textbf{4.}] This gives us the matrix of cofactors -- we multiply each term in this by the signs indicated below, to get the \textbf{adjoint matrix}: \[ \left(\begin{array}{c c c} 0.395861 & -0.181219 &-0.000813 \\ -0.181219 &0.261661 &-0.000408 \\ -0.000813 &-0.000408 &0.000035 \end{array} \right) \times \left(\begin{array}{c c c} + &- &+ \\ - &+ &- \\ + &- &+ \end{array} \right) \] \[ = \left(\begin{array}{c c c} 0.395861 &0.181219 &-0.000813 \\ 0.181219 &0.261661 &0.000408 \\ -0.000813 &0.000408 &0.000035 \end{array} \right) \] \end{itemize} } \frame{ \frametitle{Example 2.3: Solution to part(a) (7/8)} \begin{itemize} \item [\textbf{5.}] \textbf{Final step:} \begin{eqnarray*} I_{O}^{-1} &=& \frac{1}{\textsf{det}(I_{O})}\bigg(\textsf{adj}(I_{O})\bigg) \\ \pause & & \\ &=&0.00146^{-1} \times \left(\begin{array}{c c c} 0.395861 &0.181219 &-0.000813 \\ 0.181219 &0.261661 &0.000408 \\ -0.000813 &0.000408 &0.000035 \end{array} \right) \\ & & \\ &=& \left(\begin{array}{c c c} 271.1377 &124.1226 &-0.5568 \\ 124.1226 &179.2199 &0.2795 \\ -0.5568 &0.2795 &0.0240 \end{array} \right). \end{eqnarray*} \end{itemize} The standard errors are thus $\sqrt{271.1377} = 16.466$, $\sqrt{179.2199} = 13.387$ and $\sqrt{0.0240}=0.155$ for $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\xi}$, respectively. } \frame{ \frametitle{Example 2.3: Solution to part(a) (8/8)} Using the estimated standard errors, we can form confidence intervals for the GEV parameters in the usual way: \pause \vsp \begin{eqnarray*} \mu &:& 304.242\pm1.96\times16.466 \longrightarrow \pause (271.969, 336.515) \pause \\ & & \\ \sigma &:& 68.977 \pm 1.96 \times 13.387 \longrightarrow \pause (42.738, 95.216) \pause \\ & & \\ \xi &:& 0.249 \pm 1.96 \times 0.155 \longrightarrow \pause (-0.055, 0.553) \end{eqnarray*} \pause \textbf{Comment:} $\xi>0$ (i.e. we have Fr\'{e}chet tails), suggesting unbounded, heavy tails -- seems plausible for rainfall. However, the 95\% CI for $\xi$ passes through zero (only just, though!). } \frame{ \frametitle{Example 2.3: Solution to part(b) (1/1)} The correlation between $\sigma$ and $\xi$ can be estimated using \begin{eqnarray*} \textsf{corr}(\hat{\sigma},\hat{\xi}) &=& \frac{\textsf{cov}(\hat{\sigma},\hat{\xi})}{\sqrt{\textsf{var}(\hat{\sigma})\times \textsf{var}(\hat{\xi})}} \pause \\ & & \\ &=& \frac{0.2795}{\sqrt{179.2119 \times0.0240}} \pause \\ & & \\ &=& 0.1348. \end{eqnarray*} } \frame{ \frametitle{Example 2.3: Solution to part(c) (1/1)} Using Equation (10): \pause \[ \hat{z}_{r} =\hat{\mu} + \frac{\hat{\sigma}}{\hat{\xi}}\left[\left(-\textsf{log }\left(1-r^{-1}\right)\right)^{-\hat{\xi}}-1\right], \] \pause we get: \begin{eqnarray*} \hat{z}_{100} &=& 304.242 +\frac{68.977}{0.249}\left[\left(-\textsf{log }\left(1-100^{-1}\right)\right)^{-0.249}-1\right] \\ & & \\ &=&898.1133 \textsf{ mm;} \pause \qquad \textsf{similarly,} \pause \\ & & \\ \hat{z}_{1000} &=&1574.085 \textsf{ mm.} \end{eqnarray*} } \frame{ \frametitle{Example 2.3: Solution to part(d) (1/3)} From earlier, we know that \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{\sf log }y_{r}\right], \end{eqnarray*}\pause where $y_{r}=-\textrm{\sf log}(1-r^{-1})$; evaluated at the MLEs for $\mu$, $\sigma$ and $\xi$, we get \pause \[ \nabla z_{100}^{T} = \left[1,8.6097,1621.187\right] \] } \frame{ \frametitle{Example 2.3: Solution to part(d) (2/3)} Thus, \pause \begin{eqnarray*} \textsf{var}(\hat{z}_{100}) &=& \nabla z_{100}^{T}V\nabla z_{100} \\ & & \\ &=& \left[1,8.6097,1621.187\right] \\ & &\times \left(\begin{array}{c c c} 271.1377 &124.1266 &-0.5568 \\ 124.1266 &179.2199 & 0.2795 \\ -0.5568 &0.2795 &0.0240 \end{array} \right) \\ & &\times \left[\begin{array}{c} 1\\ 8.6097\\ 1621.187 \end{array}\right] \\ & & \\ &=& 84768.61 \end{eqnarray*} Therefore, $\textsf{s.e.}(\hat{z}_{100}) = \sqrt{84768.61}=291.1505$. } \frame{ \frametitle{Example 2.3: Solution to part(d) (3/3)} Similar calculations for $\hat{z}_{1000}$ give \[ \textsf{s.e.}(\hat{z}_{1000}) = 932.325. \] \pause From these standard errors, we can construct 95\% confidence intervals: \pause \begin{eqnarray*} \hat{z}_{100} &:&898.1133\pm1.96\times291.1505 \longrightarrow \pause (327.5, 1468.8) \textsf{ mm}; \pause \\ & & \\ \hat{z}_{1000} &:&1574.085 \pm 1.96\times932.325 \longrightarrow \pause (-253.3, 3401.4) \textsf{ mm.} \pause \end{eqnarray*} \textbf{Comment:} CI's are \underline{very} wide; also, the CI for $\hat{z}_{1000}$ is partly negative -- not sensible! } \normal \frame{ \frametitle{Example 2.4} The magnitudes of the vertical forces produced by seismic \textit{degassing bursts} are known to be associated with destructive volcanoes. \vsp Table 5 shows the maximum force produced by degassing bursts, every quarter in the years 2015--2021 (inclusive), for the Kilauea Volcano in Hawaii. \vsp \begin{center} {\scriptsize %\begin{table}[!h] %\centering \begin{tabular}{|c|c c c c|} \hline &Jan--Mar &Apr--Jun &Jul--Sep &Oct--Dec \\ \hline 2015 &99983 &100067 &99905 &100367 \\ 2016 &99980 &99970 &100086 &99988 \\ 2017 &99912 &100084 &100432 &99921 \\ 2018 &100123 &99913 &100240 &100448 \\ 2019 &100162 &100296 &100015 &100122 \\ 2020 &99930 &99997 &100710 &100136 \\ 2021 &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} } \end{center} } \frame{ \frametitle{Example 2.4} These data were stored the in the vector \texttt{haw} in \textsf{R}, and the following output obtained: \vsp \texttt{gev.fit(haw)} \vspl \texttt{\$ conv}\\ \texttt{[1] 0} \vsp \texttt{\$nllh}\\ \texttt{[1] 178.2502} \vsp \texttt{\$mle}\\ \texttt{[1] 9.998028e+04 8.710824e+01 5.921945e-01} \vsp \texttt{\$se}\\ \texttt{[1] 21.6185212 21.3598264 0.3027163} } \frame{ \frametitle{Example 2.4} 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? } \sol \frame{ \frametitle{Example 2.4: Solution (1/5)} From Equation (10) we know that \[ \hat{z}_{r} = \hat{\mu}+\frac{\hat{\sigma}}{\hat{\xi}}\left[\left(-\textsf{log}\left(1-r^{-1}\right)\right)^{-\hat{\xi}}-1\right]. \] The quantity $\hat{z}_{r}$ is an estimate of the level that is exceeded, on average, once every $r$ \textbf{observations} -- if we have only one observation per year, this is exactly the the same as the $r$--\underline{year} return level.\pause \vsp In this example, we have \textcolor{purple}{\bf four} observations per year! } \frame{ \frametitle{Example 2.4: Solution (2/5)} Thus, an estimate of the one--year return level is given by \pause {\small \begin{eqnarray*} \hat{z}_{4} &=& 99980.28 + \frac{87.10824}{0.5921945}\left[\left(-\textsf{log}(1-4^{-1})\right)^{-0.5921945}-1\right] \\ & & \\ &=& 100140.8 \textsf{ kg.} \end{eqnarray*} } \pause Similarly, an estimate of the 50--year return level is given by \pause {\small \begin{eqnarray*} \hat{z}_{200} &=& 99980.28 + \frac{87.10824}{0.5921945}\left[\left(-\textsf{log}(1-200^{-1})\right)^{-0.5921945}-1\right] \\ & & \\ &=& 103218.6 \textsf{ kg.} \end{eqnarray*} }} \frame{ \frametitle{Example 2.4: Solution (3/5)} Also, an estimate of the 100--year return level is given by \pause {\small \begin{eqnarray*} \hat{z}_{400} &=& 99980.28 + \frac{87.10824}{0.5921945}\left[\left(-\textsf{log}(1-400^{-1})\right)^{-0.5921945}-1\right] \\ & & \\ &=& 104940.6\textsf{ kg.} \end{eqnarray*} }} \frame{ \frametitle{Example 2.4: Solution (4/5)} We can obtain the empirical estimate of the one--year return level by looking at the \textcolor{purple}{\bf empirical distribution function} (see Section 2.2.3 on quantile--quantile plots). \pause \vsp For example, ordering the data gives: \pause {\small \begin{center} \begin{tabular}{|c| c|} \hline $x_{(i)}$ &$i/(n+1)$ \\ \hline $\vdots$ &$\vdots$ \\ 99946 &0.241 \\ 99962 &0.276 \\ $\vdots$ &$\vdots$ \\ \end{tabular} \end{center}} An empirical estimate of the one--year return level is the value which is exceeded once every 4 observations -- the value for which $i/(n+1)=0.25$. } \frame{ \frametitle{Example 2.4: Solution (5/5)} $i/(n+1)=0.25$ is 25.71\% of the way between 0.241 and 0.276; \pause 25.71\% of the way between 99946 and 99962 gives \pause \[ \hat{z}_{4,\textsf{emp}} = 99950.11 \textsf{ kg.} \] \pause The model gives a slight over--estimate. } \normal \frame{ \frametitle{Example 2.5} The table below shows the annual maximum number of sunshine hours (per day) for Santiago de Compostela, Spain, for the years 1993--2017. \begin{center} \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} \end{center} \begin{itemize} \item [(a)] Use the \textsf{ismev} package in 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} } \sol \frame{ \frametitle{Example 2.5: Solution} \begin{center} \Large \textit{See R demo} \end{center} } \normal \frame{ \frametitle{The GEV: Words of warning} \begin{itemize} \item Block size? \pause \vspl \begin{itemize} \item [--] Convention is to work with blocks equal in length to the calendar year, but what if you have hourly data collected over just a few years? \vspl \item [--] Block length too small -- limiting arguments will not hold (the GEV is a limiting result, which holds \textit{approximately} for \underline{large} $n$) \pause \vspl \item [--] Block length is too large -- not enough maxima to work with! \pause \vspl \item [--] Possible sensitivity of GEV parameter estimates to block length \pause \end{itemize} \end{itemize} } \frame{ \frametitle{The GEV: Words of warning} \begin{itemize} \item Extremely wasteful of data: \pause \vspl \begin{itemize} \item [--] Discard all but the block maxima \pause \vspl \item [--] Often results in throwing away tens of thousands of observations -- some of which might be `extreme' -- just not as extreme as the block maxima! \pause \vspl \item [--] Results in large standard errors and extremely wide confidence intervals for return levels \end{itemize} \end{itemize} } \frame{ \frametitle{The GEV: Words of warning} \begin{itemize} \item Standard asymptotic likelihood results not always applicable for the GEV: \pause \vspl \begin{itemize} \item [--] When $\xi<0$, the end--points of the GEV are functions of the parameter values: e.g. $\mu=\sigma/\xi$ is the upper end--point \pause \vspl \item [--] This violates the usual regularity conditions of MLEs \pause \vspl \item [--] Effect: When $-1<\xi\leq-0.5$, MLEs for $\xi$ can be obtained, but are \textcolor{purple}{\bf super--efficient} (i.e. have variance smaller than the Cram\'{e}r--Rao lower bound) \pause \vspl \item [--] Another consequence: Cannot obtain MLE for $\xi$ when $\xi\leq -1$ \end{itemize} \end{itemize} } \end{document} \normal \frame{ \frametitle{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. \pause \vsp Actually, it turns out that constructing confidence intervals in the standard way: \[ \textsf{estimate}\pm1.96 \times \textsf{s.e.} \] is not advisable for return levels. \pause \vsp This is because of the \textcolor{purple}{\bf 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}$. } \frame{ \frametitle{Profile likelihood} \noindent An alternative, and often more accurate, method for making inferences on a particular parameter can be found using the \textcolor{purple}{\bf profile log--likelihood}. \pause \vsp 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}$. } \frame{ \frametitle{Profile likelihood} The \textcolor{purple}{\bf profile} log--likelihood for $\boldsymbol{\theta}_{j}$ is defined as \pause \begin{eqnarray*} \ell_{p}(\boldsymbol{\theta}_{j}) &=& \textrm{\sf max}_{\boldsymbol{\theta}_{-j}} \ell(\boldsymbol{\theta}_{j},\boldsymbol{\theta}_{-j}). \end{eqnarray*}\pause 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}$. } \frame{ \frametitle{Profile likelihood} For return levels using the GEV model: \vsp \begin{enumerate} \item [1.] Re--parameterise the GEV model so that $z_{r}$ becomes one of the model parameters. \pause \vspl For example, re--arrange equation (10) to write $\mu$ in terms of $\sigma$, $\xi$ and $z_{r}$: \pause \begin{eqnarray*} \mu &=& z_{r}-\frac{\sigma}{\xi}\left[\left(-\textrm{\sf 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 (12) into \begin{eqnarray*} -m\textrm{\sf log}\sigma &-& (1+1/\xi)\sum_{i=1}^{m}\textrm{\sf 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}, \end{eqnarray*} \end{enumerate} } \frame{ \frametitle{Profile likelihood} \begin{enumerate} \item [2.] For some fixed lower value of $z_{r}=z_{r,\textrm{\sf low}}$, maximise the GEV log--likelihood \[ \ell(\sigma,\xi,z_{r}=z_{r,\textrm{\sf low}}) \] with respect to the two remaining parameters ($\sigma$ and $\xi$), to obtain $\ell_{p}(z_{r})$ at $z_{r,\textrm{\sf low}}$. \vspl Here, in terms of the more general notation above, $\boldsymbol{\theta}_{j} = z_{r}$ and $\boldsymbol{\theta}_{-j}=(\sigma,\xi)$. \end{enumerate} } \frame{ \frametitle{Profile likelihood} \begin{enumerate} \item [3.] Repeat step (2) for a range of values $z_{r}$ for $z_{r,\textrm{\sf low}}\leq z_{r}\leq z_{r,\textrm{\sf up}}$, and plot $\ell_{p}(z_{r})$ against $z_{r}$ to show the profile log--likelihood curve for $z_{r}$. \end{enumerate} \pause \vsp Figure 2.6 shows plots of the profile log likelihood for the 100 and 1000 year return levels ($\ell_{p}(z_{100})$ and $\ell_{p}(z_{1000})$ respectively) for the Lake District rainfall data in Example 2.3. } \frame{ \frametitle{Profile likelihood} \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 \tiny $z_{100}$} \psfrag{Z1000}[][]{\sf \tiny $z_{1000}$} \includegraphics[width=11.7cm,height=5.5cm]{profilelakes.ps} } \frame{ \frametitle{Profile likelihood} 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. \pause \vsp So how can we use these plots? } \frame{ \frametitle{Profile likelihood} \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 \pause \begin{eqnarray*} \boldsymbol{\theta}^{(1)}&=&z_{r} \qquad \textsf{and} \\ \boldsymbol{\theta}^{(2)}&=&(\sigma,\xi), \pause \end{eqnarray*} and the profile log--likelihood is now defined as \pause \begin{eqnarray*} \ell_{p}(\boldsymbol{\theta}^{(1)}) &=& \textrm{\sf max}_{\boldsymbol{\theta}^{(2)}}(\boldsymbol{\theta}^{(1)},\boldsymbol{\theta}^{(2)}). \end{eqnarray*} \pause The following result leads to a procedure for making inferences on the maximum likelihood estimator of $\boldsymbol{\theta}^{(1)}$. } \frame{ \frametitle{Profile likelihood: Result} 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. } \frame{ \frametitle{Profile likelihood: Result} In \textit{practice}: \pause \vspl \begin{itemize} \item Obtain a plausible range of values for $D_{p}(z_{r})$ from the $\chi_{1}^{2}$ distribution. \vsp For example, working at the 5\% level of significance, we would have values \[ D_{p}(z_{r}) \leq \chi_{1}^{2}(0.05) = 3.842. \pause \] \vspl \item Convert this into a range of plausible values for the profile log--likelihood: \pause \begin{eqnarray*} \ell(\hat{\boldsymbol{\theta}}_{0})-\ell_{p}(z_{r}) &\leq& 1.921 \\ \ell_{p}(z_{r}) &\geq& \ell(\hat{\boldsymbol{\theta}}_{0})-1.921 \end{eqnarray*} \vspl \item Use this to obtain a plausible range of values for $z_{r}$. \end{itemize} } \frame{ \frametitle{Profile likelihood: Result} \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 \tiny $z_{100}$} \psfrag{Z1000}[][]{\sf \tiny $z_{1000}$} \includegraphics[width=11.7cm,height=5.5cm]{profilelakes.ps} } \sol \frame{ \frametitle{Profile Likelihood} \begin{tabular}{|c|c|c|} \hline &$\hat{z}_{100}$ &$\hat{z}_{1000}$ \\ \hline Standard 95\% CI &(327.5, 1468.8)mm &(--253.3, 3401.4)mm \\ & & \\ \hline Profiled 95\% CI &(602.5, 2765)mm &(800, 10300)mm\\ & & \\ \hline \end{tabular} \pause \vsp \textbf{Comments} \begin{itemize} \item The profiled confidence intervals more accuratley capture the asymmetry in the log--likelihood for the return levels\pause \vspl \item This manifests in \underline{much} higher upper bounds \pause \vspl \item Notice that the lower bounds are also affected, and are much closer to the MLE for $z_{r}$ \pause \vspl \item More realistic -- and no negative rainfall values! \pause \vspl \item Problem: confidence intervals often \textit{extremely} wide, especially for long return periods (e.g. 1000 years) \end{itemize} } \normal \frame{ \frametitle{5 Models for minima} Let $\bar{M}_{n}=\textrm{\sf min}\left\{X_{1}, \ldots, X_{n}\right\}$. \vsp 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. \pause \vsp \noindent If there exist sequences of constants $a_{n}>0$ and $b_{n}$ such that, as $n \to \infty$, \begin{eqnarray*} \textrm{\sf 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: \pause \begin{eqnarray*} \bar{G}(x; \bar{\mu}, \sigma,\xi) &=& 1-\textrm{\sf exp}\left\{-\left[1-\xi\left(\frac{x-\bar{\mu}}{\sigma}\right)\right]_{+}^{-1/\xi}\right\}, \end{eqnarray*} where $a_{+}=\textrm{\sf max}(0,a)$, $-\infty<\bar{\mu}<\infty$, $\sigma>0$ and $-\infty<\xi<\infty$. } \frame{ \frametitle{5 Models for minima} This result can be useful where we are interested in modelling extremely small, rather than extremely large, observations (e.g. annual minimum air temperatures). \pause \vsp 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}$. } \frame{ \vspace{1cm} \huge \textbf{Chapter 3} \vspace{1cm} \Large \textbf{Threshold methods} } \frame{ \frametitle{.1 Background and theoretical motivation} Threshold methods use a more natural way of determining whether an observation is extreme -- {\it all} observations greater than some high value (\textcolor{purple}{\bf threshold}) are considered. \pause \vsp This allows more efficient use of data and avoids the problems that can arise as a result of blocking... \pause \vsp ...but brings its own problems (see later). \pause \vsp We must first go back and consider the asymptotic theory appropriate for this new situation. } \frame{ \frametitle{.1 Background and theoretical motivation} Suppose once more that $X_{1}, X_{2}, \ldots, X_{n}$ is a sequence of IID random variables having marginal distribution $F$, and -- once again -- let \[ M_{n} = \textrm{\sf max}\left\{X_{1}, \ldots, X_{n}\right\}. \] \pause We know from Chapter 2 that \[ \textrm{\sf Pr}\left\{M_{n}\leq x\right\} \approx G(x), \] where \pause \[ G(x) = \textrm{\sf exp}\left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}\right\} \] is the \textcolor{purple}{\bf Generalised Extreme Value} distribution with location, scale and shape $\mu$, $\sigma$ and $\xi$ respectively. } \frame{ \frametitle{.1 Background and theoretical motivation} \noindent {\large \textbf{Theorem} \textit{(Distribution of threshold excess)}} \pause \vsp For a large enough threshold $u$, the distribution function of $(X-u)$, conditional on $X>u$, is approximately \[ H(y) = 1-\left(1+\frac{\xi y}{\tilde{\sigma}}\right)_{+}^{-1/\xi}, \tag{3.1} \] defined on $y>0$, where \[ \tilde{\sigma} = \sigma +\xi(u-\mu). \tag{3.2} \] } \sol \frame{ \frametitle{Outline proof (1/5)} Denote an arbitrary term in the $X_{i}$ sequence by $X$. \pause Then it follows that a description of the behaviour of extreme events is given by \pause \begin{eqnarray*} \textsf{Pr}\left\{X>u+y|X>u\right\} &= &\frac{\textsf{Pr}(X>u+y,X>u)}{\textsf{Pr}(X>u)} \\ & & \pause\\ &=& \frac{\textsf{Pr}(X>u+y)}{\textsf{Pr}(X>u)} \pause \\ & & \\ &=& \frac{1-F(u+y)}{1-F(u)}. \qquad (\dagger) \end{eqnarray*} } \frame{ \frametitle{Outline proof (2/5)} From Section 2.1, we know that \[ \textsf{Pr}\left\{M_{n}\leq x\right\} = F^{n}(x) \approx G(x) = \textsf{exp}\left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}\right\}, \] for some parameters $\mu$, $\sigma$ and $\xi$. \pause Hence \pause \[ n\textsf{log }F(x) \approx -\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]_{+}^{-1/\xi} \qquad (\ddagger) \] } \frame{ \frametitle{Outline proof (3/5)} For large values of $x$, a Taylor series expansion implies that \pause \[ \textsf{log }F(x) \approx -\left\{1-F(x)\right\}. \] \pause Substitution into $(\ddagger)$ gives \pause \begin{eqnarray*} n \left(-\left\{1-F(u)\right\}\right) &\approx& -\left[1+\xi\left(\frac{u-\mu}{\sigma}\right)\right]_{+}^{-1/\xi} \pause \\ & & \\ 1-F(u) &\approx& \frac{1}{n}\left[1+\xi\left(\frac{u-\mu}{\sigma}\right)\right]_{+}^{-1/\xi} \end{eqnarray*} for large $u$. \pause Similarly, for $y>0$, \pause \[ 1-F(u+y) \approx \frac{1}{n}\left[1+\xi\left(\frac{u+y-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}. \] } \frame{ \frametitle{Outline proof (4/5)} Substitution into ($\dagger$) then gives \pause \begin{eqnarray*} \textsf{Pr}\left\{X>u+y|X>u\right\} &= & \frac{1-F(u+y)}{1-F(u)} \pause \\ & & \\ &\approx& \frac{\frac{1}{n}\left[1+\xi\left(\frac{u+y-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}}{\frac{1}{n}\left[1+\xi\left(\frac{u-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}} \pause \\ & & \\ &=& \left[\frac{1+\xi\left(\frac{u-\mu}{\sigma}\right) + \xi \frac{y}{\sigma}}{1+\xi\left(\frac{u-\mu}{\sigma}\right)}\right]_{+}^{-1/\xi} \pause \\ & & \\ &=&\left[\frac{1+\xi\left(\frac{u-\mu}{\sigma}\right)}{1+\xi\left(\frac{u-\mu}{\sigma}\right)} + \frac{\xi \frac{y}{\sigma}}{1+\xi\left(\frac{u-\mu}{\sigma}\right)}\right]_{+}^{-1/\xi} \\ \pause & & \\ &=&\left[1+\frac{\xi y}{\sigma+\xi(u-\mu)}\right]_{+}^{-1/\xi} \end{eqnarray*} } \frame{ \frametitle{Outline proof (5/5)} This gives \[ \textsf{Pr}\left\{X>u+y|X>u\right\} =\left[1+\frac{\xi y}{\tilde{\sigma}}\right]_{+}^{-1/\xi}, \] \pause where $\tilde{\sigma} = \sigma+\xi(u-\mu)$, as required (compare with Equation 3.1 on page 45). \pause \vsp The family of distributions defined by Equation (3.1) is known as the \textcolor{purple}{\bf Generalised Pareto family}; the distribution itself is often referred to as the \textcolor{purple}{\bf Generalised Pareto Distribution}, or GPD for short. } \normal \frame{ \frametitle{.1 Background and theoretical motivation} \begin{itemize} \item If block maxima have approximate distribution $G$, then threshold excesses have a corresponding distribution given by the Generalised Pareto family. \pause \vsp \item The parameters of the GPD are uniquely determined by those of the GEV: \pause \begin{itemize} \item [--] The parameter $\xi$ in Equation (3.1) is equal to that of the corresponding GEV \vspl \pause \item [--] The GPD scale parameter is a function of the GEV location and shape parameters \end{itemize} \end{itemize} } \frame{ \frametitle{.1 Background and theoretical motivation} \begin{itemize} \item Estimates of the GEV parameters are sensitive to the size of block chosen to identify extremes; estimates of the GPD parameters are `stable'.\pause \vsp \item The duality between the GEV and GPD means that the shape parameter $\xi$ is dominant in determining the qualitative behaviour of the GPD:\pause \vspl \begin{itemize} \item If $\xi<0$ the distribution of excesses has an upper bound \pause \vspl \item if $\xi>0$ the distribution has no upper limit \pause \vspl \item the case $\xi=0$ is also unbounded, and is taken as the limit $\xi \to 0$, giving \pause \[ H(y) = 1 -\textrm{\sf exp}\left(-\frac{y}{\tilde{\sigma}}\right), \qquad y>0; \] i.e. an exponential distribution with rate $1/\tilde{\sigma}$. \end{itemize} \end{itemize} } \frame{ \frametitle{.1 Background and theoretical motivation} \begin{itemize} \item Until this point, we have used the notation $\tilde{\sigma}$ to denote the scale parameter of the GPD, so as to distinguish it from the corresponding parameter of the GEV... \pause \vsp \item ... For notational convenience we now drop this distinction, using $\sigma$ to denote the scale parameter within either family. \end{itemize} } \frame{ \frametitle{Example 3.1} Suppose $X_{1}, X_{2}, \ldots, X_{n}$ is a sequence of independent \textit{exp}(1) random variables. \vsp Show that the limiting distribution of threshold excesses belongs to the generalised Pareto family. } \sol \frame{ \frametitle{Example 3.1: Solution} If $X_{i} \sim \textit{exp}(1)$, then $F(x) = 1-e^{-x}$ for $x>0$. \pause \vsp By direct calculation, \pause \begin{eqnarray*} \frac{1-F(u+y)}{1-F(u)} &=&\frac{1-\left(1-e^{-(u+y)}\right)}{1-\left(1-e^{-u}\right)} \pause \\ & & \\ &=& \frac{e^{-(u+y)}}{e^{-u}} = e^{-y}, \qquad y>0. \end{eqnarray*} \pause Thus, the limit distribution is an exponential distribution; i.e. a GPD with $\xi=0$. Further, we know that when $\xi=0$, we have \[ 1-H(y) = \textsf{exp}\left(-\frac{y}{\sigma}\right); \pause \] hence, we also have $\sigma=1$. That is, the limit distribution of threshold exceedances is $GPD(1,0)$. } \normal \frame{ \frametitle{Example 3.2} Suppose $X_{1}, X_{2}, \ldots, X_{n}$ is a sequence of independent $U(0,1)$ random variables. \vsp Show that the limiting distribution of threshold excesses belongs to the generalised Pareto family. } \sol \frame{ \frametitle{Example 3.2: Solution} If $X_{i} \sim U(0,1)$, then $F(x) = x$ for $0\leq x \leq 1$. \pause \vsp Hence \pause \[ \frac{1-F(u+y)}{1-F(u)} = \pause \frac{1-u+y}{1-u} = \pause 1-\frac{y}{1-u}. \] \pause For the GPD with $\xi\neq 0$, we have \pause \[ 1-H(y) = \left[1+\frac{\xi y}{\sigma}\right]^{-1/\xi}. \] Thus, we have $\xi=-1$ and $\sigma=1-u$, i.e. the limiting distribution for threshold exceedances is \textit{GPD}(1--$u$,--1). } \normal \frame{ \frametitle{.2 Application: rainfall extremes in NYC} The file \texttt{newyork.txt}, available to download from the course webpage, gives daily rainfall accumulations (in mm) for New York City for the years 1914--1961 (inclusive). \vsp Much of the northeastern United States is relatively low--lying and so prone to flooding; this is made much worse on the odd occasion that a hurricane travels this far north (for example, ``\textcolor{purple}{\bf Superstorm Sandy'}' in 2012). \vsp Thus, analysing extreme rainfall data has a real practical motivation here, in terms of \textcolor{purple}{\bf river} and \textcolor{purple}{\bf sea flood defence systems}. \vsp In this Section, we will illustrate a complete threshold--based analysis of the rainfall extremes observed at New York. } \frame{ \frametitle{.2 Application: rainfall extremes in NYC} The data have been scanned into \textbf{\sf R} from the course webpage and stored in the vector \texttt{rain}:\pause \vsp \textcolor{blue}{\texttt{rain=scan('newyork.txt')}} } \frame{ \frametitle{.2.1 Threshold choice} \noindent The \textcolor{purple}{\bf threshold stability property} of the GPD means that if the GPD is a valid model for excesses over some threshold $u_{0}$, then it is valid for excesses over all thresholds $u>u_{0}$. \pause \vsp Denoting by $\sigma_{u_{0}}$ the GPD scale parameter for excesses over threshold $u_{0}$, the expected value of our threshold excesses, conditional on being greater than the threshold, is \[ E[X-u|X>u] = \frac{\sigma_{u_{0}}+\xi u}{1-\xi}. \tag{3.3} \] \pause Thus, for all $u>u_{0}$, $E[X-u|X>u]$, is a linear function of $u$. } \frame{ \frametitle{.2.1 Threshold choice} Furthermore, $E[X-u|X>u]$ is simply the mean of the excesses of the threshold $u$, for which the sample mean of the threshold excesses of $u$ provides an estimate. \pause \vsp This leads to the \textcolor{purple}{\bf mean residual life plot}, a graphical procedure for identifying a suitably high threshold for modelling extremes via the GPD.\pause \vsp In this plot, for a range of candidate values for $u$ we identify the corresponding mean threshold excess; we then plot this mean threshold excess against $u$, and look for the value $u_{0}$ above which we can see linearity in the plot. } \frame{ \frametitle{.2.1 Threshold choice} We can easily do this from first principles in \textbf{\sf R}. \pause \vsp First of all, we set up a vector of possible thresholds, starting at zero and going up to the maximum value in our dataset: \vspl \textcolor{blue}{\texttt{> u=seq(0,max(rain),0.1)}} \pause \vsp The vector \texttt{x} is then set up to take the corresponding values for the mean excess over each value in \texttt{u}: \vspl \textcolor{blue}{\texttt{> x=vector('numeric', length(u))}} } \frame{ \frametitle{.2.1 Threshold choice} Then the following code computes the mean excess for each value in \texttt{u} and stores it in \texttt{x}: \vspl \textcolor{blue}{\texttt{> for(i in 1:length(x))}} \\ \hspace{1cm} \textcolor{red}{\texttt{\{}} \\ \hspace{1.4cm} \textcolor{red}{\texttt{threshold.exceedances=rain[rain>u[i]]}}\\ \hspace{1.4cm} \textcolor{red}{\texttt{x[i]=mean(threshold.exceedances-u[i])}} \\ \hspace{1cm} \textcolor{red}{\texttt{\}}} \pause \vsp The MRL plot is then produced using the following code, giving the plot in Figure 3.1: \vspl \textcolor{blue}{\texttt{plot(u,x,type='l', main='MRL plot',ylab='mean excess')}} } \frame{ \frametitle{.2.1 Threshold choice} %\begin{figure}[!h] %\centering \psfrag{MRL plot}[][]{\small \sf MRL plot} \psfrag{mean}[][]{\scriptsize \sf mean excess} \psfrag{u}[][]{\scriptsize \sf $u$} \psfrag{0}[][]{\tiny \sf 0} \psfrag{5}[][]{\tiny \sf 5} \psfrag{10}[][]{\tiny \sf 10} \psfrag{15}[][]{\tiny \sf 15} \psfrag{20}[][]{\tiny \sf 20} \psfrag{40}[][]{\tiny \sf 40} \psfrag{60}[][]{\tiny \sf 60} \psfrag{80}[][]{\tiny \sf 80} \includegraphics[width=10cm,height=7.5cm]{figure10.ps} %\caption{Mean residual life plot for the New York rainfall data.} %\end{figure} } \frame{ \frametitle{.2.2 Data pre--processing} In problems sheet 1, we considered how to use \textbf{\sf R} to extract the set of block maxima to be modelled by the generalised extreme value distribution. \pause \vsp Writing \textbf{\sf R} code to do this can be a time--consuming process, and the code needs to written specifically for the data being analysed. \pause \vsp In a threshold--based analysis the data pre--processing is far more straightforward. \pause \vsp Using $u_{0}=30$ as our threshold for identifying extremes (see Figure 3.1), we can easily obtain our set of threshold exceedances for modelling with the generalised Pareto distribution: \vspl \textcolor{blue}{\texttt{above.threshold=rain[rain>30]}}\\ \textcolor{blue}{\texttt{threshold.exceedances=above.threshold-30}} } \frame{ \frametitle{.2.3 Fitting the GPD} The most commonly--used approach to fit the GPD to the set of threshold excesses is that of maximum likelihood. \vsp The GPD log--likelihood function can be derived in the usual way; this is left as an exercise, but can be shown to be: \[ \ell(\sigma,\xi; \boldsymbol{y}) = -152\textrm{\sf log }\sigma - (1+1/\xi)\sum_{i=1}^{152}\textrm{\sf log}\left(1+\frac{\xi y_{i}}{\sigma}\right)_{+}, \tag{3.4} \] where $\boldsymbol{y}=(y_{1}, \ldots, y_{152})^{T}$ are the set of exceedances above threshold $u_{0}=30$. \pause \vsp For the case $\xi=0$, interpreted as $\xi \to 0$, we have the log--likelihood f or an exponential distribution with rate $1/\sigma$ (again, see problems sheet 2). } \frame{ \frametitle{.2.3 Fitting the GPD} \large \textit{R demo, including assessment of model adequacy} } \frame{ \frametitle{.2.5 Return level estimation} Figure 3.2 seems to indicate that the GPD is suitable for our set of threshold exceedances $y_{1}, \ldots, y_{152}$; that is, \pause \[ \textrm{\sf Pr}(X>u+y|X>u) \approx \left[1+\frac{\hat{\xi} y}{\hat{\sigma}}\right]_{+}^{-1/\hat{\xi}}, \tag{3.5} \] for $\xi \neq 0$. \pause Working with the LHS of (3.5), we see that \pause \[ \textrm{\sf Pr}(X>u+y|X>u) =\frac{\textrm{\sf Pr}(X>u+y)}{\textrm{\sf Pr}(X>u)}, \] \pause giving \[ \textrm{\sf Pr}(X>u+y) = \textrm{\sf Pr}(X>u)\textrm{\sf Pr}(X>u+y|X>u). \]} \frame{ \frametitle{.2.5 Return level estimation} After substitution of (3.5), we get \[ \textrm{\sf Pr}(X>u+y) \approx \hat{\lambda}_{u}\left[1+\frac{\hat{\xi} y}{\hat{\sigma}}\right]_{+}^{-1/\hat{\xi}}, \tag{3.6} \] \pause where $\lambda_{u}=\textrm{\sf Pr}(X>u)$ and is estimated as the empirical threshold exceedance rate $\hat{\lambda}_{u}$. \pause } \frame{ \frametitle{.2.5 Return level estimation} Now each $y_{i}$, $i=1,\ldots, 152$, are the raw rainfall observations (exceeding the threshold) minus the threshold ($x_{i}-u$), as the GPD models the magnitude of excess over $u$. \pause \vsp Substitution of $y_{i}=x_{i}-u$ into (3.6) gives \pause \[ \textrm{\sf Pr}(X>x) \approx \hat{\lambda}_{u}\left[1+\hat{\xi}\left(\frac{x-u}{\hat{\sigma}}\right)\right]_{+}^{-1/\hat{\xi}}. \tag{3.7} \] } \frame{ \frametitle{.2.5 Return level estimation} Thus, an estimate of the level $z_{t}$ that is exceeded on average once every $t$ \textcolor{purple}{\bf observations} is obtained as the solution of \begin{eqnarray*} \hat{\lambda}_{u}\left[1+\hat{\xi}\left(\frac{\hat{z}_{t}-u}{\hat{\sigma}}\right)\right]_{+}^{-1/\hat{\xi}} &=& \frac{1}{t}, \end{eqnarray*} \pause giving \[ \hat{z}_{t} = u +\frac{\hat{\sigma}}{\hat{\xi}}\left[(t\hat{\lambda}_{u})^{\hat{\xi}}-1\right]\tag{3.8} \] \pause for $\xi \neq 0$, and \[ \hat{z}_{t} = u +\hat{\sigma}\textrm{\sf log}(t \hat{\lambda}_{u}) \] when $\xi=0$ (see problems sheet 2). } \frame{ \frametitle{.2.5 Return level estimation} By construction, $z_{t}$ is the $t$--\textcolor{purple}{\bf observation} return level. \pause \vsp However, it is often more convenient to give return levels on an annual scale, so that the $r$--\textcolor{purple}{\bf year} return level is the level expected to be exceeded once every $r$ \textcolor{purple}{\bf years}. \pause } \frame{ \frametitle{.2.5 Return level estimation} If there are $n_{y}$ observations per year, this corresponds to the $t$--observation return level with $t=r\times n_{y}$. \pause \vsp Hence, an estimate of the $r$--year return level $z_{r}$ is defined by \[ \hat{z}_{r}=u+\frac{\hat{\sigma}}{\hat{\xi}}\left[(rn_{y}\hat{\lambda}_{u})^{\hat{\xi}}-1\right], \] \pause unless $\xi=0$, in which case \begin{eqnarray*} \hat{z}_{r}&=&u+\hat{\sigma}\textrm{\sf log}(rn_{y}\hat{\lambda}_{u}). \end{eqnarray*} } \frame{ \frametitle{.2.5 Return level estimation} Thus, for the New York City rainfall extremes, we have \[ \hat{z}_{50} = 30+\frac{7.44}{0.184}\left[(50 \times 365.25 \times 0.00867)^{0.184}-1\right] = 92.24 \textrm{\sf mm} \] as an estimate of the 50--year return level, where $n_{y}=365.25$ to account for leap years (we have daily observations). } \frame{ \frametitle{.2.5 Return level estimation} The delta method can, once again, be used to obtain estimated standard errors for such return levels. \pause \vsp We should, however, also include uncertainty in our estimate of $\lambda_{u}$ in the calculation (since $z_{r}$ is a function of $\lambda_{u}$). \pause \vsp Now the number of threshold exceedances follows a binomial distribution $\textrm{\sf Bin}(N, \lambda_{u})$, where $N$ is the total number of observations in the series. We know (from MAS2302) \pause that \[ Var(\hat{\lambda}_{u}) = \hat{\lambda}_{u}(1-\hat{\lambda}_{u})/N = 0.0007^{2} \] in our rainfall example. } \frame{ \frametitle{.2.5 Return level estimation} Recall that, by the delta method, \begin{eqnarray*} \textrm{\sf Var}(\hat{z}_{r}) &\approx& \nabla z_{r}^{T}V\nabla z_{r}. \end{eqnarray*} \pause \noindent Here, $V$ is now the variance--covariance matrix of the triple $(\hat{\lambda}_{u},\hat{\sigma},\hat{\xi})^{T}$; \pause in our rainfall example, this is \begin{eqnarray*} V &=& \left( \begin{array}{c c c} 0.0007^{2} & & \\ 0 &0.958^{2} & \\ 0 &-0.0655 &0.101^{2} \\ \end{array} \right), \end{eqnarray*} assuming that $Var(\hat{\lambda}_{u},\hat{\sigma})=Var(\hat{\lambda}_{u},\hat{\xi}) = 0$. } \sol \frame{ \frametitle{.2.5 Return level estimation (1/3)} As in Chapter 2, \pause \begin{eqnarray*} \nabla z_{r}^{T} &=& \left[ \frac{\partial z_{r}}{\partial \lambda_{u}},\frac{\partial z_{r}}{\partial \sigma},\frac{\partial z_{r}}{\partial \xi}\right]; \pause \end{eqnarray*} this can be shown to give \begin{eqnarray*} \nabla z_{r}^{T} &=& \left[\sigma (rn_{y})^{\xi} \lambda_{u}^{\xi-1}, \xi^{-1}\left\{(rn_{y} \lambda_{u})^{\xi}-1\right\} \right., \\ & &\left. -\sigma\xi^{-2}\left\{(rn_{y} \lambda_{u})-1\right\}+\sigma \xi^{-1}(rn_{y} \lambda_{u})^{\xi}\textrm{\sf log}(rn_{y} \lambda_{u})\right], \end{eqnarray*} which we evaluate at $(\hat{\lambda}_{u},\hat{\sigma},\hat{\xi})$. } \frame{ \frametitle{.2.5 Return level estimation (2/3)} This gives: \pause \[ \nabla z_{50}^{T}=(2179.096, 8.366,181.757), \] \pause and so \pause \begin{eqnarray*} Var(\hat{z}_{50}) &=& (2179.096, 8.366, 181.757) \\ & & \\ && \times \left(\begin{array}{c c c} 0.0007^{2} & & \\ 0 &0.958^{2} & \\ 0 &-0.0655 &0.101^{2} \\ \end{array} \right) \\ & &\\ & & \times \left(\begin{array}{c} 2179.096 \\ 8.366 \\ 181.757 \\ \end{array}\right) \\ \pause & & \\ &=& 397.1853 \end{eqnarray*} } \frame{ \frametitle{.2.5 Return level estimation (3/3)} And so \[ s.e. (\hat{z}_{50}) = \sqrt{397.1853} = 19.930. \] } \normal \frame{ \frametitle{.2.6 Profile likelihood} As discussed in Chapter 2, confidence intervals for return levels are better constructed via the method of profile likelihood, owing to the asymmetry in the likelihood surface often observed. \vsp A plot of the profile log--likelihood for the 50--year return level for daily rainfall accumulations at New York is shown in Figure 3.3 of the lecture notes. \vsp This gives \[ (74.1\textsf{ mm}, 143\textsf{ mm}). \] \textit{Interpretation: ``Once every fifty years, we might expect daily rainfall accumulations in New York City to reach up to about 143mm''}. } \frame{ \frametitle{.2.6 Profile likelihood} Compare the 95\% confidence interval for the 50--year return level obtained from profiling the log--likelihood to that you would obtain using the standard error. Comment. } \sol \frame{ \frametitle{.2.6 Profile likelihood} From the profiled log--likelihood: \[ (74.1\textsf{ mm}, 143\textsf{ mm}). \] Using the standard error, we get: \[ 92.24 \pm 1.96\times 19.93 \longrightarrow (53.2\textsf{ mm},131.3\textsf{ mm}). \] Comment: the confidence interval based on the profiled log--likelihood is much more conservative, giving a substantially higher upper bound for the return level. } \normal \frame{ \vspace{1cm} \huge \textbf{Chapter 4} \vspace{1cm} \Large \textbf{Dependent extremes} } \frame{ \frametitle{.1 Introduction} The threshold--based approach to modelling extremes has obvious practical advantages over the more traditional `block maxima' approach. \vsp 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. \vsp 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$. } \frame{ \frametitle{.1 Introduction} \textbf{Advantages} \pause \vspl \begin{itemize} \item Data pre--processing dead easy \pause \vspl \item Estimates with \textcolor{purple}{\bf smaller standard errors} (including return levels!) \pause \vspl \item \textcolor{purple}{\bf Narrower confidence intervals} (standard or profile--likelihood based) \end{itemize} \pause \vsp \textbf{Disadvantages} \pause \vspl \begin{itemize} \item \textcolor{purple}{\bf Temporal dependence} (e.g. temperature extremes, wind gusts,...)\pause \vspl \item \textcolor{purple}{\bf Non--stationarity} (e.g. trends and seasonal effects) \end{itemize} } \frame{ \frametitle{.2 Example: Hourly wind speeds at High Bradfield} Hourly maximum wind gusts (in \textcolor{red}{\bf knots}) were collected at High Bradfield, in the Peak District, over a period of 10 years from January 1st 2003 to December 31st 2012. \pause \vsp This gives total of 87,672 observations (including any missing values: 81,835 after the missing values have been removed). \pause \vsp 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}. } \frame{ \frametitle{.2 Example: Hourly wind speeds at High Bradfield} \begin{center} \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=12cm,height=4.5cm]{brad.ps} \end{center} } \frame{ \frametitle{.2 Example: Hourly wind speeds at High Bradfield} \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=12cm,height=7cm]{brad2.ps} \end{center} } \frame{ \frametitle{.2 Example: Hourly wind speeds at High Bradfield} It is clear that we have (fairly strong) dependence between consecutive wind speeds recordings -- even at extreme levels. \pause \vsp How can we proceed? \pause \vspl \begin{itemize} \item \textcolor{purple}{\bf Ignore} this dependence -- proceed as normal? \pause \vspl \textbf{Possible consequence:} Standard errors under--estimated as there are fewer independence observations than the likelihood assumes \pause \vspl \item \textcolor{purple}{\bf Model} this dependence -- see Chapter 6 \pause \vspl \item \textcolor{purple}{\bf Filter out} this dependence \end{itemize} } \frame{ \frametitle{Maxima of dependent (but stationary) series} \noindent The book by \textcolor{purple}{\bf Leadbetter \textit{et al.} (1983)} considers, in great detail, properties of extremes of dependent processes. \pause \vsp A key result often used is `\textcolor{purple}{\bf Leadbetter's $\boldsymbol{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. \pause \vsp This condition is stated more formally in the definition below. } \frame{ \frametitle{Maxima of dependent (but stationary) series} {\large\textbf{Definition}\textit{ \textcolor{brown}{(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$,\pause {\small \begin{eqnarray} \Big\arrowvert\textrm{\sf 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{\sf Pr}\left\{\tilde{X}_{i_{1}}\leq u_{n}, \ldots, \tilde{X}_{i_{p}}\leq u_{n}\right\}\textrm{\sf Pr}\left\{\tilde{X}_{j_{1}}\leq u_{n}, \ldots,\tilde{X}_{j_{q}}\leq u_{n}\right\}\Big\arrowvert&\leq& \alpha(n,l), \nonumber \end{eqnarray}} \vsp where $\alpha(n,l) \to 0$ for some sequence $l_{n}$ such that $l_{n}/n \to 0$ as $n \to \infty$. } \frame{ \frametitle{Maxima of dependent (but stationary) series} \noindent For sequences of independent variables, the difference in probabilities in the above expression is exactly zero for \textit{any} sequence $u_{n}$. \vsp 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$. \vsp 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. } \frame{ \frametitle{Maxima of dependent (but stationary) series} {\large \textbf{Theorem} \textit{\textcolor{brown}{(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{\sf max}\{\tilde{X}_{1}, \ldots, \tilde{X}_{n}\}$. \vsp Now let $X_{1}, X_{2}, \ldots$ be an \textcolor{purple}{\bf independent} series with $X$ having the same distribution as $\tilde{X}$, and let $M_{n}=\textrm{\sf max}\{X_{1}, \ldots, X_{n}\}$. Then if $M_{n}$ has a non--degenerate limit law given by $\textrm{\sf Pr}\left\{(M_{n}-b_{n})/a_{n} \leq x\right\} \to G(x)$, it follows that \[ \textrm{\sf Pr}\left\{(\tilde{M}_{n}-b_{n})/a_{n} \leq x\right\} \to G^{\theta}(x) \tag{4.2} \] for some $0 \leq \theta \leq 1$. } \frame{ \frametitle{Maxima of dependent (but stationary) series} The parameter $\theta$ is known as the \textcolor{purple}{\bf extremal index}, and quantifies the extent of extremal dependence \pause \vspl \begin{itemize} \item $\theta=1$: completely independent process \pause \vspl \item $\theta \to 0$: increasing levels of (extremal) dependence \end{itemize} } \frame{ \frametitle{Maxima of dependent (but stationary) series} Since $G$ in the above theorem is necessarily an extreme value distribution, and due to the \textcolor{purple}{\bf 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. \pause \vsp The powering of the limit distribution by $\theta$ only affects the location and scale parameters of this distribution. } \frame{ \frametitle{Maxima of dependent (but stationary) series} What does this all mean in \textcolor{purple}{\bf practical terms}? \pause \vspl \begin{itemize} \item If maxima of a stationary series converge in distribution (and we know they do -- the GEV), and \pause \vspl \item if Leadbetter's $D(u_{n})$ condition holds (i.e. long--range dependence is negligible), then \pause \vspl \item the limit distribution is related to that of the independent series... \pause \vspl \item ...in fact, it's $G^{\theta}(x)$, where $G$ is the GEV for the independent series (with parameters $\mu$, $\sigma$ and $\xi$) \end{itemize} } \frame{ \frametitle{Maxima of dependent (but stationary) series} Further, it can be shown that \begin{eqnarray*} G^{\theta}(x) &=& \textrm{\sf exp}\left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]^{-1/\xi}\right\}^{\theta} \\ &=& \textrm{\sf 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}$. } \frame{ \frametitle{.1 Modelling block maxima} We can proceed as we did in Chapter 2 -- even in the presence of dependence... \pause \vsp ... provided long range dependence is weak (i.e. \textcolor{purple}{\bf Leadbetter's $\boldsymbol{D(u_{n})}$ condition} can be assumed) \pause \vsp The only difference is the parameters of the GEV: $(\mu, \sigma, \xi)\longrightarrow (\mu^{*},\sigma^{*}, \xi)$. \pause \vsp Since we estimates these parameters anyway, who cares?! } \frame{ \frametitle{.2 Modelling threshold exceedances} 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. \pause \vsp 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. } \frame{ \frametitle{.2 Modelling threshold exceedances} In practice, several techniques have been developed to circumvent this problem, including: \vspl \begin{itemize} \item [\textbf{1.}] filtering out an (approximately) independent set of threshold exceedances \vspl \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 \vspl \item [\textbf{3.}] Explicitly modelling the temporal dependence in the process \end{itemize} } \frame{ \frametitle{.2 Declustering} Since the mid--1990s, various methods for \textcolor{purple}{\bf declustering} a series of extremes, to extract a set of independent extremes, have been discussed in the literature. \pause \vsp The most natural, commonly--used method of declustering is that of \textcolor{purple}{\bf runs declustering}. \pause This is how it works: \pause \vspl \begin{itemize} \item [1.] Choose an auxiliary `declustering parameter' (which we call $\kappa$) \vspl \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 \vspl \item [3.] Go through the entire series identifying clusters in this way \vspl \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} } \frame{ \frametitle{.2 Declustering} This approach is often referred to as the \textcolor{purple}{\bf peaks over threshold} approach (POT, \textcolor{brown}{\bf Davison and Smith, 1990}) and is widely accepted as the main pragmatic approach for dealing with clustered extremes. \pause \vsp Easy to implement, but if: \pause \vspl \begin{itemize} \item $\kappa$ is too small, the cluster peaks will not be far enough apart to safely assume independence \pause \vspl \item $\kappa$ is too large, there will be too few cluster exceedances on which to form our inference \end{itemize} } \frame{ \frametitle{.2 Declustering} It has also been shown (\textcolor{brown}{\bf 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. } \frame{ \frametitle{.2 Declustering} \textbf{Using $\boldsymbol{\kappa=2}$:} \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=10cm,height=6cm]{brad2base.ps} \end{center} } \frame{ \frametitle{.2 Declustering} \textbf{Using $\boldsymbol{\kappa=2}$:} \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=10cm,height=6cm]{brad2i.ps} \end{center} } \frame{ \frametitle{.2 Declustering} \textbf{Using $\boldsymbol{\kappa=2}$:} \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=10cm,height=6cm]{brad2ii.ps} \end{center} } \frame{ \frametitle{.2 Declustering} \textbf{Using $\boldsymbol{\kappa=2}$:} \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=10cm,height=6cm]{brad2iii.ps} \end{center} } \frame{ \frametitle{.2 Declustering} \textbf{Using $\boldsymbol{\kappa=4}$:} \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=10cm,height=6cm]{brad2base.ps} \end{center} } \frame{ \frametitle{.2 Declustering} \textbf{Using $\boldsymbol{\kappa=4}$:} \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=10cm,height=6cm]{brad2iv.ps} \end{center} } \frame{ \frametitle{.2 Declustering} \textbf{Using $\boldsymbol{\kappa=4}$:} \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=10cm,height=6cm]{brad2iii.ps} \end{center} } \frame{ \frametitle{.2 Declustering} \textbf{Using $\boldsymbol{\kappa=10}$:} \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=10cm,height=6cm]{brad2base.ps} \end{center} } \frame{ \frametitle{.2 Declustering} \textbf{Using $\boldsymbol{\kappa=10}$:} \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=10cm,height=6cm]{brad2v.ps} \end{center} } \frame{ \frametitle{.2 Declustering} Suppose we use $\kappa=10$. \vsp 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!). \vsp Unfortunately, there is no function in \textsf{ismev} to perform runs declustering. \vsp I have written the following \textbf{\sf R} code to perform this declustering, for $\kappa=10$. } \frame{ \frametitle{.2 Declustering} \textit{Class demonstration in R} } \frame{ \frametitle{Words of warning (\textit{Fawcett \& Walshaw, 2012})} \psfrag{20}[][]{\sf \tiny 20} \psfrag{25}[][]{\sf \tiny 25} \psfrag{30}[][]{\sf \tiny 30} \psfrag{35}[][]{\sf \tiny 35} \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{k=10}[][]{\sf \tiny $\boldsymbol{\kappa=10}$} \psfrag{k=30}[][]{\sf \tiny \textcolor{blue}{$\boldsymbol{\kappa=30}$}} \psfrag{Profile}[][]{\sf \scriptsize Profile log-likelihood} \psfrag{1000 year return level}[][]{\sf \scriptsize 1000 year return level} \psfrag{q_1000}[][]{\sf \tiny $z_{1000}$ (metres)} \centering \includegraphics[height=2.75in,width=4.25in]{kappa1030.ps} } \frame{ \frametitle{Words of warning (\textit{Fawcett \& Walshaw, 2012})} \psfrag{20}[][]{\sf \tiny 20} \psfrag{25}[][]{\sf \tiny 25} \psfrag{30}[][]{\sf \tiny 30} \psfrag{35}[][]{\sf \tiny 35} \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{k=10}[][]{\sf \tiny $\boldsymbol{\kappa=10}$} \psfrag{k=30}[][]{\sf \tiny \textcolor{blue}{$\boldsymbol{\kappa=30}$}} \psfrag{Profile}[][]{\sf \scriptsize Profile log-likelihood} \psfrag{1000 year return level}[][]{\sf \small 1000 year return level} \psfrag{q_1000}[][]{\sf \tiny $z_{1000}$ (metres)} \centering \includegraphics[height=2.75in,width=4.25in]{kappa1030ii.ps} } \frame{ \frametitle{Words of warning (\textit{Fawcett \& Walshaw, 2012})} \psfrag{20}[][]{\sf \tiny 20} \psfrag{25}[][]{\sf \tiny 25} \psfrag{30}[][]{\sf \tiny 30} \psfrag{35}[][]{\sf \tiny 35} \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{k=10}[][]{\sf \tiny $\boldsymbol{\kappa=10}$} \psfrag{k=30}[][]{\sf \tiny \textcolor{blue}{$\boldsymbol{\kappa=30}$}} \psfrag{Profile}[][]{\sf \scriptsize Profile log-likelihood} \psfrag{1000 year return level}[][]{\sf \small 1000 year return level} \psfrag{q_1000}[][]{\sf \tiny $z_{1000}$ (metres)} \centering \includegraphics[height=2.75in,width=4.25in]{kappa1030iii.ps} } \frame{ \frametitle{Words of warning (\textit{Fawcett \& Walshaw, 2012})} \psfrag{20}[][]{\sf \tiny 20} \psfrag{25}[][]{\sf \tiny 25} \psfrag{30}[][]{\sf \tiny 30} \psfrag{35}[][]{\sf \tiny 35} \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{k=10}[][]{\sf \tiny $\boldsymbol{\kappa=10}$} \psfrag{k=30}[][]{\sf \tiny \textcolor{blue}{$\boldsymbol{\kappa=30}$}} \psfrag{Profile}[][]{\sf \scriptsize Profile log-likelihood} \psfrag{1000 year return level}[][]{\sf \small 1000 year return level} \psfrag{q_1000}[][]{\sf \tiny $z_{1000}$ (metres)} \centering \includegraphics[height=2.75in,width=4.25in]{kappa1030iv.ps} } \frame{ \frametitle{Words of warning (\textit{Fawcett \& Walshaw, 2012})} \psfrag{20}[][]{\sf \tiny 20} \psfrag{25}[][]{\sf \tiny 25} \psfrag{30}[][]{\sf \tiny 30} \psfrag{35}[][]{\sf \tiny 35} \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{k=10}[][]{\sf \tiny $\boldsymbol{\kappa=10}$} \psfrag{k=30}[][]{\sf \tiny \textcolor{blue}{$\boldsymbol{\kappa=30}$}} \psfrag{Profile}[][]{\sf \scriptsize Profile log-likelihood} \psfrag{1000 year return level}[][]{\sf \small 1000 year return level} \psfrag{q_1000}[][]{\sf \tiny $z_{1000}$ (metres)} \centering \includegraphics[height=2.75in,width=4.25in]{kappa1030v.ps} } \frame{ \frametitle{Words of warning (\textit{Fawcett \& Walshaw, 2012})} \psfrag{20}[][]{\sf \tiny 20} \psfrag{25}[][]{\sf \tiny 25} \psfrag{30}[][]{\sf \tiny 30} \psfrag{35}[][]{\sf \tiny 35} \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{k=10}[][]{\sf \tiny $\boldsymbol{\kappa=10}$} \psfrag{k=30}[][]{\sf \tiny \textcolor{blue}{$\boldsymbol{\kappa=30}$}} \psfrag{Profile}[][]{\sf \scriptsize Profile log-likelihood} \psfrag{1000 year return level}[][]{\sf \small 1000 year return level} \psfrag{q_1000}[][]{\sf \tiny $z_{1000}$ (metres)} \centering \includegraphics[height=2.75in,width=4.25in]{kappa1030vi.ps} } \frame{ \frametitle{Words of warning (\textit{Fawcett \& Walshaw, 2012})} \psfrag{5}[][]{\sf \tiny 5} \psfrag{10}[][]{\sf \tiny 10} \psfrag{15}[][]{\sf \tiny 15} \psfrag{20}[][]{\sf \tiny 20} \psfrag{25}[][]{\sf \tiny 25} \psfrag{30}[][]{\sf \tiny 30} \psfrag{35}[][]{\sf \tiny 35} \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{kappa}[][]{\sf \tiny $\boldsymbol{\kappa}$} \psfrag{sdsd}[][]{\sf \tiny \textcolor{blue}{$\boldsymbol{\kappa=30}$}} \psfrag{Profile}[][]{\sf \scriptsize Profile log-likelihood} \psfrag{1000 year return level}[][]{\sf \small 1000 year return level} \psfrag{z_1000}[][]{\sf \tiny $z_{1000}$ (metres)} \centering \includegraphics[height=2.75in,width=4.25in]{kappaall.ps} } \sol \frame{ \frametitle{Words of warning} \begin{itemize} \item Similar to the block maxima approach, this is \textcolor{purple}{\bf wasteful of data}! We are throwing away observations we have identified as extreme! \pause \vspl \item How do we choose $\kappa$? Estimates of GPD parameters -- and return levels/confidence intervals for return levels -- can be sensitive to the choice of $\kappa$ (often hugely so) \pause \vspl \item What if we are interested in the dependence? \end{itemize} } \announce \frame{ \frametitle{Welcome back!} \begin{itemize} \item Assignments \pause \vsp \item Mid--semester test \pause \vspl \begin{itemize} \item [--]This Thursday, \textcolor{purple}{\bf prompt 4.05pm} start \pause \vspl \item [--]Specimen, with solutions, available online \pause \vspl \item [--]Can cover anything from \textcolor{purple}{\bf Chapters 1--4} \pause \vspl \item [--]Key formulae given on paper (and exam paper in May/June): CDF for GEV and GPD \pause \vspl \item [--]Tips: should know the \textcolor{purple}{\bf delta method}; should know the link between the \textcolor{purple}{\bf Hessian matrix}, the \textcolor{purple}{\bf information matrix} and the \textcolor{purple}{\bf variance--covariance matrix}; should be familiar with \textbf{\sf R} output for things like \texttt{gev.fit} and \texttt{gpd.fit}; worthwhile memorising formulae for \textcolor{purple}{\bf return levels} \end{itemize} \end{itemize} } \frame{ \frametitle{Welcome back!} \begin{itemize} \item Project \pause \vspl \begin{itemize} \item [--]Due in 4pm, \textcolor{purple}{\bf Friday 10th May} \pause \vspl \item [--]Based on your own unique dataset, and the paper I gave out before Easter \pause \end{itemize} \vsp \item Next few weeks \pause \vspl \begin{itemize} \item [--] Chapter 5: non--stationarity (e.g. trend and/or seasonality) \pause \vspl \item [--] Chapter 6: Multivariate considerations \end{itemize} \end{itemize} } \normal \frame{ \vspace{1cm} \huge \textbf{Chapter 5} \vspace{1cm} \Large \textbf{Non--stationary extremes} } \frame{ \frametitle{5.1 Introduction} In the context of environmental processes, it is common to observe \textcolor{purple}{\bf non--stationarity}: \pause \vspl \begin{itemize} \item different \textcolor{purple}{\bf seasons} having different climate patterns \pause \vspl \item \textcolor{purple}{\bf long term trends} owing to climate change \pause \end{itemize} \vsp The models from chapters 2 and 3 assume that the observations used are IID! } \frame{ \frametitle{5.1 Introduction} We examined the effects of \textcolor{purple}{\bf dependence} in Chapter 4: \pause \vspl \begin{itemize} \item Likely that consecutive extremes of wind or rain will be \textcolor{purple}{\bf correlated} \pause \vspl \item Block maxima: Doesn't matter -- provided we can assume long--range independence (Leadbetter's $D(u_{n})$ condition), then the GEV still appropriate \pause \vspl \item Threshold exceedances: Not quite so straightforward, although we can \textcolor{purple}{\bf decluster}, and fit the GPD to \textcolor{purple}{\bf cluster peak excesses}? Problems? \end{itemize} } \frame{ \frametitle{5.1 Introduction} To date, no general theory for \textcolor{purple}{\bf non--stationary extremes} has been established. \pause \vsp In practice, it is common to adopt pragmatic `workarounds' based on the type of non--stationarity observed. \pause \vsp For this reason, in this Chapter we will give some specific examples of how practitioners have dealt with non--stationarity in recent work and publications. } \frame{ \frametitle{5.2 Annual maximum sea levels in Venice} \noindent Recall \textcolor{purple}{\bf question 4} in \textcolor{purple}{\bf problems sheet 2}, which investigated the use of the $r$--largest order statistics model for extreme sea levels at Venice (1961--2011). \pause \vsp These data are available to download from the course webpage. } \frame{ \frametitle{5.2 Annual maximum sea levels in Venice} \vspace{2cm} \begin{center} {\Large \it Demonstration in \textbf{\sf R}} \end{center} } \frame{ \frametitle{5.2 Annual maximum sea levels in Venice} One way of capturing the trend observed in the Venice annual maxima is to allow the GEV location parameter $\mu$ to vary across time. \pause \vsp From Figure 5.1, a \textcolor{purple}{\bf simple linear trend} in time seems plausible for our annual maximum sea levels $X$, and so we could use the model \[ X_{t} \sim \textrm{\sf GEV}(\mu(t),\sigma,\xi), \] \pause where \[ \mu(t) = \beta_{0}+\beta_{1}t \tag{5.1} \] and $t$ is an indicator of year. } \frame{ \frametitle{5.2 Annual maximum sea levels in Venice} \[ \mu(t) = \beta_{0}+\beta_{1}t: \] \begin{itemize} \item Variations over time modelled as a linear trend in the \textcolor{purple}{\bf location parameter} \pause \vspl \item $\beta_{1}$ represents the \textcolor{purple}{\bf slope} -- in this case, the annual rate of change in sea--level maxima at Venice \pause \vspl \item The time--homogeneous model is a special case of this time dependent model, with $\beta_{1}=0$... \pause \vspl \item ...since this is \textcolor{purple}{\bf nested} within the model which allows for a time dependence, the deviance statistic can be used to formally compare models (see later) \end{itemize} } \frame{ \frametitle{5.2.1 Parameter estimation} Recall the log--likelihood function for the GEV: \begin{eqnarray*} \ell(\mu,\sigma,\xi;\boldsymbol{x}) &=& -m\textrm{\sf log}\sigma-(1/\xi+1)\sum_{i=1}^{m}\textrm{\sf log}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+} \\ & & \hspace{0.2cm} -\sum_{i=1}^{m}\left[1+\xi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}, \end{eqnarray*} where $m$ is the number of block maxima $x_{1}, x_{2}, \ldots, x_{m}$. } \frame{ \frametitle{5.2.1 Parameter estimation} We simply replace $\mu$ in the above expression with equation (5.1), giving \begin{eqnarray*} \ell(\beta_{0},\beta_{1},\sigma,\xi;\boldsymbol{x},t) &=& -m\textrm{\sf log}\sigma \\ & & \hspace{0.2cm} -(1/\xi+1)\sum_{i=1}^{m}\textrm{\sf log}\left[1+\xi\left(\frac{x_{i}-(\beta_{0}+\beta_{1}t_{i})}{\sigma}\right)\right]_{+}\\ & &\hspace{0.2cm} -\sum_{i=1}^{m}\left[1+\xi\left(\frac{x_{i}-(\beta_{0}+\beta_{1}t_{i})}{\sigma}\right)\right]_{+}^{-1/\xi}, \end{eqnarray*} with the usual replacement when $\xi=0$. } \frame{ \frametitle{5.2.1 Parameter estimation} \begin{itemize} \item \textit{Could} maximise $\ell(\beta_{0},\beta_{1},\sigma,\xi;\boldsymbol{x},t)$ `from first principles' using \texttt{nlm} \pause \vspl \item However, easier to do in the \textsf{ismev} package! \end{itemize} \vsp \begin{center} {\Large \it Demonstration in \textbf{\sf R}} \end{center} } \frame{ \frametitle{5.2.1 Parameter estimation} Thus, we have \[ \hat{\beta}_{0} = 96.986 (4.249) \quad \hat{\beta}_{1} = 0.564 (0.139) \] \[ \hat{\sigma} = 14.584 (1.578) \quad \hat{\xi} = -0.027 (0.083) \] (with standard errors in parentheses). \pause \vsp This gives: \[ \hat{\mu}(t) = 96.986+0.564t, \] giving an estimated increase in maximum sea levels at Venice of about 0.564cm per year. } \frame{ \frametitle{5.2.1 Parameter estimation} For example, the estimated value for $\mu$ in the year 2013 would be \[ \hat{\mu}(53) = 96.986+0.564\times 53 = 126.878; \pause \] we could, of course, use the delta method to obtain the corresponding standard error. \pause \vspl Using this simple linear model, we can estimate $\mu$ for $t=1,2, \ldots, 51$ to cover the years for which we have data (i.e. 1961, 1962, \ldots, 2011). } \frame{ \frametitle{5.2.1 Parameter estimation} \begin{center} \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{2010}[][]{\sf \tiny 2010} \psfrag{80}[][]{\sf \tiny 80} \psfrag{100}[][]{\sf \tiny 100} \psfrag{120}[][]{\sf \tiny 120} \psfrag{140}[][]{\sf \tiny 140} \psfrag{160}[][]{\sf \tiny 160} \psfrag{180}[][]{\sf \tiny 180} \psfrag{y}[][]{\sf \scriptsize Year} \psfrag{sl}[][]{\sf \scriptsize Sea level (cm)} \includegraphics[width=12cm,height=8cm]{venicectrend.ps} \end{center} } %\caption{Time series plot showing the annual maximum sea levels (in cm above the average) with the trend for $\mu$ superimposed (dashed red line)} %\end{figure}} \frame{ \frametitle{5.2.2 Model choice} Fitting a completely stationary model to the set of annual maxima, as you did in question 4(a) of problems sheet 2, gives: \vsp \texttt{gev.fit(venice.anmax)} \vspl \texttt{\$conv}\\ \texttt{[1] 0}\\ \vspl \texttt{\$nllh}\\ \texttt{[1] 222.7145} \vspl \texttt{\$mle}\\ \texttt{[1] 111.09925486 17.17548761 -0.07673265} \vspl \texttt{\$se} \\ \texttt{[1] 2.6280071 1.8033672 0.0735214} } \frame{ \frametitle{5.2.2 Model choice} Questions: \pause \begin{itemize} \item Is the non--stationary model \textcolor{purple}{\bf worthwhile?} \pause \vspl \item Is the trend we observe in Figure 5.1 \textcolor{purple}{\bf significant}? \pause \vspl \item In other words, does the non--stationary model provide an \textcolor{purple}{\bf improvement} in fit over the simpler model shown here? \end{itemize} } \frame{ \frametitle{5.2.2 Model choice} We can use a version of the result in Section 2.4 (page 41) to address this question. \pause \vsp Generally, maximum likelihood estimation of nested models leads to a simple test procedure of one model against the other. \pause \vsp With models $\mathcal{M}_{0} \subset \mathcal{M}_{1}$, we define the deviance statistic as \[ D = 2\left\{\ell_{1}(\mathcal{M}_{1})-\ell_{0}(\mathcal{M}_{0})\right\}, \] where $\ell_{1}(\mathcal{M}_{1})$ and $\ell_{0}(\mathcal{M}_{0})$ are the maximised log--likelihood under models $\mathcal{M}_{1}$ and $\mathcal{M}_{0}$ respectively. \pause \vsp Asymptotically, $D \sim \chi^{2}_{k}$, where $k=\textsf{dim}(\mathcal{M}_{1})-\textsf{dim}(\mathcal{M}_{0})$. } \frame{ \frametitle{5.2.2 Model choice} Formally, \[ H_{0}: \mathcal{M}_{1}=\mathcal{M}_{0}. \] \pause Thus, large values of $D$ suggest that the increase in model size has been worthwhile, and that model $\mathcal{M}_{1}$ explains substantially more of the variation in the data than $\mathcal{M}_{0}$. } \sol \frame{ \frametitle{5.2.2 Model choice} Suppose \pause \vspl \begin{itemize} \item $\mathcal{M}_{1}$: model allowing for \textcolor{purple}{\bf linear trend} in $\mu$ \pause \vspl \item $\mathcal{M}_{0}$: \textcolor{purple}{\bf stationary} model \end{itemize} } \frame{ \frametitle{5.2.2 Model choice} Thus, \[ D = 2\left\{-216.0626-(-222.7145)\right\} = 13.3038. \pause \] \pause The difference in dimensionality is 1, so the critical value for our test is \[ \chi^{2}_{1}(0.05) = 3.841. \pause \] \vsp Relative to this, our deviance statistic is large: the model which allows for a linear trend in $\mu$ explains substantially more of the variation in the data than does the simpler stationary model. } \normal \frame{ \frametitle{5.2.2 Model choice} We could, of course, use this method to check for a more complex association through time. \pause \vsp For example, to check for a \textcolor{purple}{\bf quadratic trend} we might use a model $\mathcal{M}_{2}$ with the following form for $\mu$: \pause \[ \mu(t) = \beta_{0}+\beta_{1}t + \beta_{2}t^{2}. \] \pause \begin{center} {\Large \it Demonstration in \textbf{\sf R}} \end{center} } \frame{ \frametitle{5.2.2 Model choice} Comparing model $\mathcal{M}_{2}$ with model $\mathcal{M}_{1}$ gives \[ D = 2\left\{-216.0555 - (-216.0626)\right\} = 0.0142, \] which is \textit{small} compared to $\chi^{2}_{1}(0.05)=3.841$. \pause \vsp Thus, allowing for a quadratic dependence in time does not improve on our model which allows for a linear trend through time, and so we would reject model $\mathcal{M}_{2}$. } \frame{ \frametitle{5.2.3 Model diagnostics} Before estimating return levels, we should check the goodness--of--fit of our model which allows for a linear trend in $\mu$. \pause \vsp The lack of homogeneity in the distributional assumptions for each observation, however, mean some modification of the standard procedures (e.g. \textcolor{purple}{\bf probability plots} and \textcolor{purple}{\bf quantile plots}) is required. } \frame{ \frametitle{5.2.3 Model diagnostics} For example, for the Venice annual maximum sea levels, we have \[ X_{t} \sim \textrm{\sf GEV}(\mu(t), \sigma, \xi), \quad t=1, 2, \ldots, 51, \] giving a different GEV in each year indicated by $t$. \pause \vsp What we need to do is standardise so that we can assume the $X_{t}$ \textit{are} IID across all years $t$. \pause \vsp Usually, the set of non--stationary annual maxima are transformed to a common Gumbel distribution with distribution function $F(y) = \textrm{\sf exp}\left\{-e^{-y}\right\}$. } \sol \frame{ \frametitle{5.2.3 Model diagnostics} We can obtain the required transformation by equating $F(y_{t})$ to $\textrm{\sf GEV}(x_{t};\hat{\mu}(t),\hat{\sigma},\hat{\xi})$ and solving for $y_{t}$: \begin{eqnarray*} \textrm{\sf exp}\left\{-e^{-y_{t}}\right\} &=&\textrm{\sf exp}\left\{-\left[1+\hat{\xi}\left(\frac{x_{t}-(\hat{\beta}_{0}+\hat{\beta}_{1}t)}{\hat{\sigma}}\right)\right]^{-1/\hat{\xi}}\right\} \\ \pause && \\ e^{-y_{t}} &=& \left[1+\hat{\xi}\left(\frac{x_{t}-(\hat{\beta}_{0}+\hat{\beta}_{1}t)}{\hat{\sigma}}\right)\right]^{-1/\hat{\xi}} \\ \pause & & \\ y_{t} &=&\displaystyle\frac{1}{\hat{\xi}}\textrm{\sf log}\left[1+\hat{\xi}\left(\frac{x_{t}-(\hat{\beta}_{0}+\hat{\beta}_{1}t)}{\hat{\sigma}}\right)\right] , \end{eqnarray*} } \frame{ \frametitle{5.2.3 Model diagnostics} This gives gives the following transformation to common Gumbel margins for our Venice annual maxima: \[ y_{t} = -\frac{1}{0.027}\textrm{\sf log}\left[1-0.027\left(\frac{x_{t}-(96.986+0.564t)}{14.584}\right)\right]. \] \vsp \begin{center} {\Large \it Demonstration in \textbf{\sf R}} \end{center} } \normal \frame{ \frametitle{5.2.3 Model diagnostics} Now that we have transformed the original data to a single common distribution, we can apply the standard graphical diagnostics. \pause \vsp For example, we can compare the \textcolor{purple}{\bf empirical probabilities} and \textcolor{purple}{\bf quantiles} of $y_{t}$ to their \textcolor{purple}{\bf theoretical} counterparts from the Gumbel distribution. \pause \vsp Fortunately, the usual command in \textsf{ismev} can be used to produce these plots for us. \pause \begin{center} {\Large \it Demonstration in \textbf{\sf R}} \end{center}} \frame{ \frametitle{5.2.4 Return level estimation} Recall Equation (2.10) (\textcolor{red}{\bf memorise!}) from Chapter 2 for estimating return levels from the GEV: \[ \hat{z}_{r} = \hat{\mu}+\frac{\hat{\sigma}}{\hat{\xi}}\left[\left(-\textrm{\sf \sf log}\left(1-r^{-1}\right)\right)^{-\hat{\xi}}-1\right]. \] \pause Since we have a time--varying location parameter $\hat{\mu}(t) = \hat{\beta}_{0}+\hat{\beta}_{1}t$, we will clearly have time--varying estimates of return levels $\hat{z}_{r}(t)$. \pause \vsp For example, an estimate of the sea level we might expect to see in Venice once every 100 years is given by \[ \hat{z}_{r}(t) = (96.986+0.564 t)-\frac{14.584}{0.027}\left[\left(-\textrm{\sf log}\left(1-100^{-1}\right)\right)^{0.027}-1\right]. \] } \frame{ \frametitle{5.2.4 Return level estimation} Figure 5.4 shows how we might expect this estimate to vary for $t=52, 53, \ldots$, i.e. for the years 2012, 2013, \ldots. \pause \vsp We could treat these as forecasts of the 100--year return levels as we move through time. \pause \vsp Obviously, such forecasts will assume the linear trend for $\mu$ continues beyond the range of data we have observed and will, of course be subject to error (which we can estimate by constructing point--wise 95\% confidence intervals using the profile log--likelihood, for example). } \frame{ \frametitle{5.2.4 Return level estimation} \begin{center} %begin{figure}[!h] %\centeri \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{2010}[][]{\sf \tiny 2010} \psfrag{2020}[][]{\sf \tiny 2020} \psfrag{2030}[][]{\sf \tiny 2030} \psfrag{2040}[][]{\sf \tiny 2040} \psfrag{2050}[][]{\sf \tiny 2050} \psfrag{80}[][]{\sf \tiny 80} \psfrag{100}[][]{\sf \tiny 100} \psfrag{120}[][]{\sf \tiny 120} \psfrag{140}[][]{\sf \tiny 140} \psfrag{160}[][]{\sf \tiny 160} \psfrag{180}[][]{\sf \tiny 180} \psfrag{190}[][]{\sf \tiny 190} \psfrag{195}[][]{\sf \tiny 195} \psfrag{200}[][]{\sf \tiny 200} \psfrag{205}[][]{\sf \tiny 205} \psfrag{210}[][]{\sf \tiny 210} \psfrag{y}[][]{\sf \scriptsize Year} \psfrag{sl}[][]{\sf \scriptsize Sea level (cm)} \includegraphics[width=12cm,height=7.5cm]{venice100.ps} %\caption{Time series plot showing the forecasted 100--year return levels at Venice, for the years 2012--2050} %\end{figure} \end{center} } \frame{ \frametitle{5.3 Sea levels and the SOI} A different situation which could use the same approach as that in the previous section is where the extremal behaviour of a series is related to another variable, rather than time. \pause \vsp For example, studies have revealed a link between annual maximum sea levels at Fremantle, Australia, and the mean value of the \textcolor{purple}{\bf Southern Oscillation Index} (SOI, an indicator of meteorological volatility due to effects such as El Ni\~{n}o). } \frame{ \frametitle{5.3 Sea levels and the SOI} \begin{center} \psfrag{1900}[][]{\sf \tiny 1900} \psfrag{1920}[][]{\sf \tiny 1920} \psfrag{1940}[][]{\sf \tiny 1940} \psfrag{1960}[][]{\sf \tiny 1960} \psfrag{1980}[][]{\sf \tiny 1980} \psfrag{-1}[][]{\sf \tiny -1} \psfrag{0}[][]{\sf \tiny 0} \psfrag{1}[][]{\sf \tiny 1} \psfrag{2}[][]{\sf \tiny 2} \psfrag{1.2}[][]{\sf \tiny 1.2} \psfrag{1.4}[][]{\sf \tiny 1.4} \psfrag{1.6}[][]{\sf \tiny 1.6} \psfrag{1.8}[][]{\sf \tiny 1.8} \psfrag{l}[][]{\sf \tiny Sea level (m)} \psfrag{s}[][]{\sf \tiny SOI} \psfrag{y}[][]{\sf \tiny Year} \includegraphics[width=11cm,height=6cm]{soi.ps} \end{center} } \frame{ \frametitle{5.3 Sea levels and the SOI} Thus, the following model for $X_{t}$, the annual maximum sea level at Fremantle in year $t$, might be suitable: \pause \[ X_{t} \sim \textrm{\sf GEV}(\mu(t), \sigma, \xi), \] where \[ \mu(t) = \beta_{0}+\beta_{1}\textrm{\sf SOI}(t), \tag{5.3} \] where $\textrm{\sf SOI}(t)$ denotes the mean value of SOI in year $t$. } \frame{ \frametitle{5.3 Sea levels and the SOI} However, the plot in the right--hand--side of Figure 5.5 also reveals a possible trend in sea levels through time, suggesting \[ \mu(t) = \beta_{0}+\beta_{1}t, \tag{5.4} \] where $t=1, 2, \ldots,$ as in Example 5.2. } \frame{ \frametitle{5.3 Sea levels and the SOI} We can combine Equations (5.3) and (5.4) to allow for a dependence on time \textit{and} SOI by letting \[ \mu(t) = \beta_{0}+\beta_{1}\textrm{\sf SOI}(t)+\beta_{2}t; \tag{5.5} \] however, a technique of \textcolor{purple}{\bf forward selection} should be used to check whether or not any of Equations (5.3), (5.4) or (5.5) give significant improvement over the stationary model. } \frame{ \frametitle{5.3 Sea levels and the SOI} \vspace{1.7cm} \begin{center} {\large \textit{Class demonstration in \textbf{\sf R}}} \end{center} } \frame{ \frametitle{5.4 Rainfall extremes in New York City} Recall Section 3.2 in which we modelled rainfall extremes in New York using a threshold--based approach. \vsp The Generalised Pareto distribution was applied to rainfall exceedances over a threshold of $u=30$mm, giving estimates of the scale and shape as \pause \[ \hat{\sigma}=7.44 (0.958) \quad \hat{\xi}=0.184 (0.101) \] (respectively). } \frame{ \frametitle{5.4 Rainfall extremes in New York City} Recall that the $\textrm{\sf GPD}(\sigma,\xi)$ arises from the $\textrm{\sf GEV}(\mu,\sigma,\xi)$, where the GPD scale parameter is a function of the GEV location and shape parameters. \pause \vsp Thus, attempting to model any trend in our threshold exceedances is usually done through linear modelling of the scale parameter $\sigma$ (the GPD doesn't have a location parameter \textit{per se}). } \frame{ \frametitle{5.4 Rainfall extremes in New York City} Since the scale parameter $\sigma$ must be positive, we might choose to model a trend through time as \[ \sigma(t) = \textrm{\sf exp}\{\beta_0+\beta_1t\}, \tag{5.6} \] where $t$ is, once again, an indicator of time. \begin{center} {\large \textit{Class demonstration in \textbf{\sf R}}} \end{center} } \frame{ \frametitle{5.5 Generalisation} With reference to the examples discussed so far, we could model non--stationarity through \textit{any} of the parameters in our extremal model. \pause \vsp For example, take a non--stationary GEV model to describe the distribution of $X_{t}$, for $t=1, 2, \ldots, m$: \[ X_{t} \sim \textrm{\sf GEV}(\mu(t), \sigma(t), \xi(t)), \] where each of the model parameters have an expression in terms of a parameter vector $\boldsymbol{\beta}$ and some covariates. \pause \vsp The likelihood is then \[ L(x_{t};\boldsymbol{\beta}) = \prod_{i=1}^{m}g(x_{t};\mu(t),\sigma(t),\xi(t)), \] where $g$ is the GEV density function. } \frame{ \frametitle{5.5 Generalisation} From this, we can form the log--likelihood, and then maximise in the usual way (for example, using \texttt{nlm} in \textbf{\sf R}). \pause \vsp In terms of threshold exceedances $Y_{t}$, $t=1, 2, \ldots, k$, we could replace the GEV with the GPD: \[ Y_{t} \sim \textrm{\sf GPD}(\sigma(t),\xi(t)), \] with $\sigma$ being defined as in Equation (\textcolor{red}{\bf 5.6}) to retain the positivity of the GPD scale. } \frame{ \frametitle{Chapters 4 \& 5: Departures from IID: The story so far} \begin{itemize} \item \textbf{Temporal dependence:} block maxima approach -- still OK; threshold approach -- decluster using POT \pause \vsp \item \textbf{Trend:} Use the deviance statistic to check for trends in GEV location parameter/GPD scale parameter \pause \vsp \item \textbf{Covariates:} Use the deviance statistic, as above \pause \vsp \item \textbf{Seasonal variability} \end{itemize} } \frame{ \frametitle{5.6 Wind speed extremes at High Bradfield} %\begin{figure}[!h] \begin{center} \psfrag{5000}[][]{\sf \tiny 5000} \psfrag{10000}[][]{\sf \tiny 10000} \psfrag{15000}[][]{\sf \tiny 15000} \psfrag{20000}[][]{\sf \tiny 20000} \psfrag{25000}[][]{\sf \tiny 25000} \psfrag{30000}[][]{\sf \tiny 30000} \psfrag{35000}[][]{\sf \tiny 35000} \psfrag{0}[][]{\sf \tiny 0} \psfrag{20}[][]{\sf \tiny 20} \psfrag{40}[][]{\sf \tiny 40} \psfrag{60}[][]{\sf \tiny 60} \psfrag{80}[][]{\sf \tiny 80} \psfrag{w}[][]{\sf \tiny Wind speed (knots)} \psfrag{h}[][]{\sf \tiny Hour} \includegraphics[width=10cm,height=6cm]{season.ps} %\caption{Time series plot showing the first four years of wind speed data at High Bradfield in the Peak District.} %\end{figure} \end{center}} \frame{ \frametitle{5.6 Wind speed extremes at High Bradfield} %\begin{figure}[!h] \begin{center} \psfrag{0}[][]{\sf \tiny 0} \psfrag{20}[][]{\sf \tiny 20} \psfrag{40}[][]{\sf \tiny 40} \psfrag{60}[][]{\sf \tiny 60} \psfrag{80}[][]{\sf \tiny 80} \psfrag{100}[][]{\sf \tiny 100} \psfrag{1}[][]{\sf \tiny Jan} \psfrag{2}[][]{\sf \tiny Feb} \psfrag{3}[][]{\sf \tiny Mar} \psfrag{4}[][]{\sf \tiny Apr} \psfrag{5}[][]{\sf \tiny May} \psfrag{6}[][]{\sf \tiny Jun} \psfrag{7}[][]{\sf \tiny Jul} \psfrag{8}[][]{\sf \tiny Aug} \psfrag{9}[][]{\sf \tiny Sep} \psfrag{10}[][]{\sf \tiny Oct} \psfrag{11}[][]{\sf \tiny Nov} \psfrag{12}[][]{\sf \tiny Dec} \psfrag{w}[][]{\sf \tiny Wind speed (knots)} \psfrag{m}[][]{\sf \tiny Month} \includegraphics[width=10cm,height=6cm]{season2.ps} %\caption{Boxplots of wind speed data at High Bradfield, by month.} %\end{figure} \end{center}} \frame{ \frametitle{5.6.1 Single season approach} Under the single season approach, an extremal model is fitted to the extremes of an environmental process from the season which gives rise to the `most extreme' extremes. \pause \vsp \begin{itemize} \item Bradfield wind speeds -- use January only \vspl \item Some appeal: \textcolor{purple}{\bf easy to implement}, focuses on the \textcolor{purple}{\bf most extreme} extremes \vspl \item But: \textcolor{purple}{\bf wasteful} of data -- especially if we decluster as well! \end{itemize} } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} \noindent It is usual in strongly seasonal climates for the occurrence of extreme winds to be confined to a certain part of the yearly cycle. \vsp Example: U.K. $\longrightarrow$ unusual for wind damage to occur outside the period \textcolor{purple}{\bf October--March}. } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} The seasonal variation observed at Bradfield might be expected: \pause Prolonged, \textcolor{purple}{\bf anticyclonic} periods more prevalent during June, July and August. \vsp There is only a point to modelling the extremes which occur in summer months if we believe that they can help us understand what happens in winter months, where genuinely large events \textit{can} occur. } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} For this to be realistic, we must assume that the \textcolor{purple}{\bf same mechanism} is responsible for the generation of large gusts throughout the year \pause $\longrightarrow$ it is just the \textcolor{purple}{\bf scale} of this mechanism which changes. \vsp Indeed, in temperate climates (such as that of the U.K.), the same alternating sequence of anticyclones and depressions leads to most of the storms which occur throughout the year. \vsp The seasonal variability comes from the \textcolor{purple}{\bf severity} of these systems. } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} Taking the \textcolor{purple}{\bf calendar month} as our seasonal unit: \pause \vsp \begin{itemize} \item For the U.K., experience suggests this might be OK \vspl \item Strikes a good balance between: \vspl \begin{itemize} \item [--] reflecting the continuous nature of seasonal changes in climate \vspl \item [--] retaining a substantial amount of data for analysis within each season \end{itemize} \end{itemize} } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} {\small \begin{center} \begin{tabular}{c c c c c} \hline Month ($m$) &$u_{m}$ &$n_{m}$ &$\hat{\sigma}_{m}$ &$\hat{\xi}_{m}$ \\ \hline 1 &55.341 &28 &21.373 (5.358) &$-0.420$ (0.183) \\ 2 &41.531 &24 &15.130 (4.635) &$-0.226$ (0.233) \\ 3 &48.100 &29 &23.277 (6.316) &$-0.894$ (0.259) \\ 4 &39.910 &29 &14.853 (4.448) &$-0.440$ (0.249) \\ 5 &31.943 &46 &9.456 (1.990) &$-0.158$ (0.147) \\ 6 &35.670 &35 &12.329 (2.592) &$-0.409$ (0.143) \\ 7 &32.290 &36 &12.517 (2.609) &$-0.605$ (0.161) \\ 8 &32.639 &34 &10.199 (2.361) &$-0.203$ (0.159) \\ 9 &33.232 &49 &18.772 (3.668) &$-0.255$ (0.138) \\ 10 &44.914 &34 &11.669 (3.533) &$-0.274$ (0.254) \\ 11 &48.394 &33 &14.991 (3.381) &$-0.225$ (0.149) \\ 12 &49.341 &35 &18.681 (4.229) &$-0.416$ (0.166) \\ \hline \end{tabular} \end{center} } } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} In terms of \textcolor{purple}{\bf return level inference}, it would not make practical sense to have monthly varying estimates of the $r$--year return level $z_{r}$. \pause \vsp To include information from \textit{all} months in our return level estimation procedure, we solve \begin{eqnarray*} \prod_{m=1}^{12} H(\hat{z}_{r};\hat{\lambda}_{u_{m}},\hat{\sigma}_{m},\hat{\xi}_{m}) &=& 1-\frac{1}{rn_{y}} \end{eqnarray*} for $\hat{z}_{r}$, where $H$ is the GPD distribution function and $n_{y}$ is the (average) number of observations per year. } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} Thus, we need to solve \begin{eqnarray*} \prod_{m=1}^{12}\left\{1-\hat{\lambda}_{u}\left[1+\hat{\xi}\left(\frac{\hat{z}_{r}-u_{m}}{\hat{\sigma}_{m}}\right)\right]_{+}^{-1/\hat{\xi}_{m}}\right\} &=& 1-\frac{1}{rn_{y}} \end{eqnarray*} \pause for $\hat{z}_{r}$, that is \begin{eqnarray*} \prod_{m=1}^{12}\left\{1-\hat{\lambda}_{u}\left[1+\hat{\xi}\left(\frac{\hat{z}_{r}-u_{m}}{\hat{\sigma}_{m}}\right)\right]_{+}^{-1/\hat{\xi}_{m}}\right\} -\left(1-\frac{1}{rn_{y}}\right) &=&0. \end{eqnarray*} } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} This equation cannot be solved analytically. \pause \vsp Rather, a numerical procedure must be used. \pause \vsp This is easy to do in \textbf{\sf R} using the \texttt{uniroot} procedure. } \begin{frame} \frametitle{5.6.2 Seasonal piecewise approach} {\scriptsize \texttt{f = }\textcolor{blue}{\rm function}\texttt{(z)}$\{$\\ \hspace{0.4cm} \texttt{r} \textbf{=} \texttt{10}\\ \hspace{0.4cm} \texttt{ny} \textbf{=} \texttt{365.25*24}\\ \textcolor{brown}{\hspace{0.4cm} \texttt{sigma} \textbf{=} \texttt{c(21.373,15.130,23.277,..., 18.681)}}\\ \textcolor{brown}{\hspace{0.4cm} \texttt{xi} \textbf{=} \texttt{c(-0.420,-0.226,-0.894,-0.440,..., -0.416)}}\\ \hspace{0.4cm} \texttt{u} \textbf{=} \texttt{c(55.341,41.531,48.1,...,49.341)}\\ \hspace{0.4cm} \texttt{nobs} \textbf{=} \texttt{c(31*24*10,28*24*7+29*24*3,...,31*24*10)}\\ \textcolor{brown}{\hspace{0.4cm} \texttt{lambda} \textbf{=} \texttt{c(28/nobs[1],24/nobs[2],...,35/nobs[12])}}\\ \hspace{0.4cm} \texttt{component} \textbf{=} \texttt{\textrm{\color{red} vector}}\texttt{(`numeric',12)} \\ \hspace{0.4cm} \texttt{inner} \textbf{=} \texttt{\textrm{\color{red} vector}}\texttt{(`numeric',12)} \\ \vspace{0.1cm} \hspace{0.9cm} \textcolor{blue}{\rm for}\texttt{(m in 1:12)}$\{$ \\ \hspace{1.5cm} \texttt{inner[m]} \textbf{=} \texttt{\textrm{\color{red}max}((1+xi[m]*((z-u[m])/sigma[m])),0)} \\ \hspace{1.5cm} \texttt{component[m]} \textbf{=} \texttt{1-lambda[m]*((inner[m])**(-1/xi[m]))}$\}$\\ \vspace{0.1cm} \hspace{0.4cm} \texttt{answer} \textbf{=} \textrm{\color{red} prod}\texttt{(component)-(1-1/(r*ny))} \\ \hspace{0.4cm} \textrm{\color{blue} return} \texttt{(answer)}$\}$ } \end{frame} \begin{frame}[fragile] \frametitle{5.6.2 Seasonal piecewise approach} \begin{verbatim} > uniroot(f,lower=0,upper=200) $root [1] 102.3291 $f.root [1] -6.242917e-11 $iter [1] 12 $estim.prec [1] 6.103516e-05 \end{verbatim} \end{frame} \frame{ \frametitle{5.6.2 Seasonal piecewise approach} %\begin{table} {\tiny \begin{center} \begin{tabular}{|c|c c c c|} \hline &\multicolumn{4}{c|}{Return period ($r$ years)} \\ & 10 &50 &200 &1000 \\ \hline $\hat{z}_{r}$ (st. err.) &102.33 (3.970) &104.59 (15.951) &106.21 (23.793) &108.89 (44.865) \\ \hline \end{tabular} \end{center} } } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} For the standard errors, we use the delta method. However, we now have: \pause \begin{eqnarray*} V&=& \left( \begin{array}{c c c c c c} \frac{\hat{\lambda}_{u_{1}}(1-\hat{\lambda}_{u_{1}})}{N_{1}} &0 &\ldots &0 &\ldots &0 \\ 0 &\ddots &\ddots &\vdots &\ddots &\vdots \\ \vdots &\ddots &\frac{\hat{\lambda}_{u_{12}}(1-\hat{\lambda}_{u_{12}})}{N_{12}} &0 &\ldots &0 \\ 0 &\ldots &0 &v_{1,1} &\ldots &v_{1,24} \\ \vdots &\ddots &\vdots &\vdots &\ddots &\vdots \\ 0 &\ldots &0 &v_{24,1} &\ldots &v_{24,24} \end{array} \right), \end{eqnarray*} where $v_{i,j}$ denotes the $(i,j)$--th term of the variance--covariance matrix of $\hat{\sigma}_{m}$ and $\hat{\xi}_{m}$, $m=1, \ldots, 12$. } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} Hence, by the delta method, \begin{eqnarray*} \textrm{\sf Var}(\hat{z}_{r}) \approx \nabla z_{r}^{T}V\nabla z_{r}, \end{eqnarray*} where \begin{eqnarray*} \nabla z_{r}^{T} &=& \left[\frac{\partial z_{r}}{\partial \lambda_{u_{1}}}, \ldots, \frac{\partial z_{r}}{\partial \lambda_{u_{12}}}, \frac{\partial z_{r}}{\partial \sigma_{1}}, \ldots, \frac{\partial z_{r}}{\partial \sigma_{12}}, \frac{\partial z_{r}}{\partial \xi_{1}}, \ldots, \frac{\partial z_{r}}{\partial \xi_{12}}\right], \end{eqnarray*} evaluated at $(\hat{\lambda}_{u_{1}},\hat{\sigma}_{1},\hat{\xi}_{1}, \ldots, \hat{\lambda}_{u_{12}},\hat{\sigma}_{12},\hat{\xi}_{12})$. } \frame{ \frametitle{5.6.2 Seasonal piecewise approach} As an aside, notice how the estimate of the 50--year return level wind speed here differs to that when we assumed stationarity in Section 4.3.2: \pause \begin{eqnarray*} \textsf{Assuming stationarity} &:& \hat{z}_{50} = 101.53 \textsf{ knots} \\ \pause & & \\ \textsf{Seasonal piecewise} &:& \hat{z}_{50} = 104.59 \textsf{ knots} \end{eqnarray*} \pause Such discrepancies become more pronounced for longer--period return levels. } \frame{ \frametitle{5.6.3 Smoothly varying seasonal parameters} Various authors (e.g. Fawcett and Walshaw (2006)) have investigated the use of continuously varying parameters for the GPD when seasonal variation is present. \pause \vsp For example, \textcolor{purple}{\bf Fourier forms} can be used to allow the GPD scale and shape to vary smoothly through time. \pause \vsp However, such analyses for the Bradfield wind speed data yielded little, if any, improvement over the seasonal piecewise approach.} \announce \frame{ \frametitle{Mid-semester test} \begin{center} \begin{tabular}{c c c c} \hline Mean &Median &Max &St. dev. \\ \hline 38.9 &40.0 &46.0 &7.1 \\ \hline \end{tabular} \end{center} \pause \noindent \textbf{Question 1} \begin{itemize} \item [(a)] TRUE / FALSE questions: \begin{itemize} \item [--] Generally well done \vspl \item [--] \textcolor{red}{"The 50 year return level for daily temperatures in Seville, Spain, is $44.7^{\textrm{o}}$C. Thus, in the next 100 years we will see temperatures of at least $44.7^{\textrm{o}}$C on two days"} \vspl \item [--] \textcolor{red}{"A Fr\'{e}chet distribution might be appropriate for data on sunshine hours as the limiting distribution has a finite upper bound"} \end{itemize} \end{itemize} } \frame{ \frametitle{Mid-semester test} \begin{itemize} \item [(b)] Suspicious results: {\scriptsize \begin{center} \begin{tabular}{|c|c c c|c c c|} \hline &\multicolumn{3}{c|}{Estimated parameters} &\multicolumn{3}{c|}{Return level estimates $\hat{z}_{r}$} \\ & Location &Scale &Shape &$r=10$ &$r=50$ & $r=200$ \\ \hline GEV &67.2 &7.3 &--1.2 &73.1 &72.5 &85.3 \\ & (13.3) &(2.2) &(0.2) &(10.0) &(11.3) &(13.5) \\ \hline GPD & --- &6.8 &0.1 &73.6 &78.3 &86.2 \\ & &(1.9) &(0.1) &(15.2) &(7.1) &(6.9) \\ \hline \end{tabular} \end{center}} \end{itemize} } \frame{ \frametitle{Mid-semester test} \hspace{-0.5cm}\noindent \textbf{Question 2} \pause \begin{itemize} \item [(a)(i)] Typical answer: \textrm{--10.7145} \pause \\ Correct answer: $\textrm{exp}\{\textrm{--10.7145}\} = \textrm{0.0000222}$ \pause \vspl \item [(a)(ii)] Must mention the satisfactory convergence code of 0 here \pause \vspl \item [(b)(ii)] Common mistakes: \pause \vspl \begin{itemize} \item [--]Finding the MLE for the exponential rate $\zeta$ but not the GPD scale $\sigma=1/\zeta$ \pause \vspl \item [--] Not showing you have found the maximum of the log-likelihood function \pause \vspl \item [--] Using the exponential distribution function instead of the density function to construct the likelihood \end{itemize} \pause \vspl \item [(b)(iii)] All weird and wonderful things... \pause \vspl \item [(c)\textcolor{white}{(iii)}] Should really mention \textcolor{purple}{dependence} and \textcolor{purple}{seasonal variability} \end{itemize} } \announce \frame{ \frametitle{Stuff} \begin{itemize} \item Assignment 3 - last one - due in by 4pm, Thurs 12th May. Section A \underline{or} Section B!\pause \vspl \item Mid-semester test: feedback last week, slides on webpage/BB \pause \vspl \item Story so far: \pause \begin{itemize} \item [1.] Practical motivation \pause \vspl \item [2.] Limiting models for IID block maxima \pause \vspl \item [3.] Limiting models for IID threshold excesses \pause \vspl \item [4.] Departures from IID: \pause \textcolor{purple}{\bf Dependence} (wind speeds in the Peak District), \pause \textcolor{purple}{\bf trend} (Venice sea levels, New York rainfall), \pause \textcolor{purple}{\bf covariate effects} (sea surge and SOI), \pause \textcolor{purple}{\bf seasonality} (back to the wind speeds...) \pause \vspl \item [5.] Bayesian inference for extremes \underline{or} multivariate extremes... \end{itemize} \end{itemize} } \normal \frame{ \vspace{1cm} \huge \textbf{Chapter 6} \vspace{1cm} \Large \textbf{Multivariate extremes} } \frame{ \frametitle{6.1 Introduction} In this section we consider the problems we face if we wish to model the extremal behaviour of two or more (dependent) processes \textcolor{purple}{\bf simultaneously}.\pause \vsp There are several reasons why we may wish to do this:\pause \vsp \begin{itemize} \item to model the extreme behaviour of a \textcolor{purple}{\bf particular variable} over \textcolor{purple}{\bf several nearby locations} (e.g. rainfall over a network of sites -- simultaneous flooding at several locations could cause devastation); \end{itemize} } \frame{ \frametitle{6.1 Introduction} \begin{itemize} \item to model the joint extremes of two or more \textcolor{purple}{\bf different variables} at a \textcolor{purple}{\bf 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); \pause \vsp \item to model the joint behaviour of extremes which occur as \textcolor{purple}{\bf consecutive observations} in a time--series (e.g. consecutive hourly maximum wind gusts during a storm). \end{itemize} } \frame{ \frametitle{6.1. Introduction} All of these problems suggest fitting an appropriate \textcolor{purple}{\bf limiting multivariate distribution} to the relevant data. \vsp However, as we shall see, the derivation of such a multivariate distribution is not as easy as we might hope. \pause \vsp The analogy with the Normal distribution as a model for means breaks down as we move into $n$ dimensions! \vsp It is not even clear what the ` \textcolor{purple}{\bf relevant data}' should be! \vsp Most of the increased complexity is apparent in the move from $1$ to $2$ dimensions, so we will focus largely on \textcolor{purple}{\bf bivariate problems}. } \frame{ \frametitle{6.2.1 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 in southwest Scotland shown in Figure 6.1. \pause \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} } \frame{ \frametitle{6.2.1 Example: Rainfall at 8 locations in Scotland} Such issues are of great interest, especially currently -- given the \textcolor{purple}{\bf severe flooding} experienced in the U.K. \vsp Suppose we have sequences of daily total rainfall at each location. \vsp There is likely to be \textcolor{purple}{\bf strong inter--site dependence} between extremes, in the sense that days with heavy rain are liable to occur simultaneously across locations. \vsp The raw multivariate observations are $8$--dimensional vectors of the daily rainfall over the eight sites. } \frame{ \frametitle{6.2.1 Example: Rainfall at 8 locations in Scotland} Now suppose we wish to take a block--maxima approach, with `blocks' being years. \vsp For any given year, the 8--dimensional vector of annual maxima is unlikely to be one of the raw multivariate observations. \vsp 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)$. } \frame{ \frametitle{6.2.1 Example: Rainfall at 8 locations in Scotland} Now consider the \textcolor{purple}{\bf componentwise block maxima} \[ M_{x,n} = \max_{i=1,\ldots,n}\{X_i\} \qquad \mbox{\sf and} \qquad M_{y,n} = \max_{i=1,\ldots,n}\{Y_i\}. \] \pause We define the \textcolor{purple}{\bf 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$.} \frame{ \frametitle{6.2.1 Example: Rainfall at 8 locations in Scotland} The first point to note is that \textcolor{purple}{\bf standard univariate extreme value results} apply in each margin. \pause \vsp When considering the dependence, this allows us to make a simplifying assumption. } \frame{ \frametitle{6.2.1 Example: Rainfall at 8 locations in Scotland} 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{\sf exp}(-1/z), \qquad z>0. \] This gives rise to a very simple normalisation of maxima: \[ \textrm{\sf Pr} (X_i < x) = \textcolor{white}{\textrm{\sf Pr} (M_{x,n}/n < x)} = \exp(-1/x), \qquad x>0, \] \textcolor{white}{ (and similarly for $Y_i$). So if we consider the re--scaled vector} \[ {\color{white}{\bf M_n^*} = \left(\max_{i=1,\ldots,n}\{X_i\}/n, \max_{i=1,\ldots,n}\{Y_i\}/n\right),} \] \textcolor{white}{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. } } \frame{ \frametitle{6.2.1 Example: Rainfall at 8 locations in Scotland} 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{\sf exp}(-1/z), \qquad z>0. \] This gives rise to a very simple normalisation of maxima: \[ \textrm{\sf Pr} (X_i < x) = \textrm{\sf Pr} (M_{x,n}/n < x) = \exp(-1/x), \qquad x>0, \pause \] (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. } \frame{ \frametitle{6.2.2 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. \pause \vsp Then if \[ \textrm{\sf 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 \tag{6.1} \end{equation} \pause 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) \tag{6.2} \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. \tag{6.3} \end{equation} } \frame{ \frametitle{6.2.2 Limiting distributions for bivariate extremes} Hence the class of bivariate extreme value distributions is in one-to-one correspondence with distribution functions $H$ satisfying the constraint (6.3). \pause \vsp If $H$ is differentiable with density $h$, then (6.2) 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*} \pause 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$. } \frame{ \frametitle{6.2.2 Limiting distributions for bivariate extremes} Since the GEV provides the complete class of \textcolor{purple}{\bf marginal} limit distributions, then the complete class of \textcolor{purple}{\bf bivariate} extreme value distributions is obtained as follows. \pause \vsp 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. } \frame{ \frametitle{6.2.2 Limiting distributions for bivariate extremes} Hence \[ G(x,y) = \exp\{-V(\tilde{x}, \tilde{y})\} \] is a bivariate extreme value distribution with the appropriate margins for valid $V(.)$.} \frame{ \frametitle{6.3.1 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. \pause \vsp Specifically, we consider \textcolor{purple}{\bf two} variables recorded concurrently at the \textcolor{purple}{\bf same site}.\pause \vsp A series of 3--hourly measurements on sea--surge were obtained from Newlyn, southwest England, giving, at each time point,\pause \vspl \begin{itemize} \item [(i)] measurements of the \textcolor{purple}{\bf wave height} (in metres) \pause \vspl \item [(ii)] measurements of the \textcolor{purple}{\bf surge height} (in metres) \end{itemize} } \frame{ \frametitle{6.3.1 Example: wave--surge data at Newlyn, Cornwall} Figure 6.2 shows these two variables plotted against each other. \vsp This plot suggests a tendency for extremes of one variable to coincide with extremes of the other. \vsp This dependence could be important -- the impact of an event that is \textcolor{purple}{\bf simultaneously extreme} in both variables is likely to be much greater than if extremes of either component occurred in isolation. } \frame{ \frametitle{6.3.1 Example: wave--surge data at Newlyn, Cornwall} %\begin{figure}[!h] \begin{center} \includegraphics[height=10cm,width=7cm,angle=-90]{wavesurge.ps} \end{center} %\caption{Concurrent wave and surge heights at Newlyn, Cornwall.} %\end{figure} } \frame{ \frametitle{6.3.2 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. \vsp For our bivariate observation $(X,Y)$, let's focus on $X$. \pause \vsp 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: \pause \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{ \sf 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{\sf Pr}\left(X>u_{x}\right)$. } \frame{ \frametitle{6.3.2 Componentwise maxima to threshold excesses} Now we can obtain a unit Fr\'{e}chet margin with the \textcolor{purple}{\bf 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}. \tag{6.4} \end{equation} \pause If we apply the analogous transformation in 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).} \frame{ \frametitle{6.4 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. \vsp Models can be fitted, e.g. by maximum--likelihood estimation, either: \pause \vspl \begin{itemize} \item in \textcolor{purple}{\bf two steps} (marginal components followed by dependence function), or \pause \vspl \item in a \textcolor{purple}{\bf single sweep} \end{itemize} \vspl 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. } \frame{ \frametitle{6.4.1 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. \pause \vsp The dependence structure must satisfy the conditions on $H\left(\omega\right)$. \vsp Possible choices are: \pause \vspl \begin{itemize} \item \textcolor{purple}{\bf Logistic Model} --- symmetric \vspl \item Negative Logistic Model \vspl \item \textcolor{purple}{\bf Bilogistic Model} --- asymmetric \vspl \item Dirichlet Model\\[-6mm] \end{itemize} \vspl Here we will focus on the logistic model and the bilogistic model as two commonly used but contrasting choices. } \sol \frame{ \frametitle{6.4.1 Modelling the dependence structure} \textbf{Symmetric dependence} \pause \vspl \begin{itemize} \item $X$ depends on $Y$ to exactly the same degree that $Y$ depends on $X$ \pause \vspl \item They both have the same \textcolor{purple}{\bf influence} over each other \end{itemize} \pause \vsp \textbf{Asymmetric dependence} \vspl \begin{itemize} \item $X$ has greater influene over $Y$ than $Y$ has over $X$ (or \textit{vice--versa}) \pause \vspl \item Example: Wind speeds at two nearby locations in the U.K. \end{itemize} } \normal \frame{ \frametitle{. 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)$, giving \[ G(x,y) = \textsf{exp}\left\{-\left(x^{-1/\alpha}+y^{-1/\alpha}\right)^{\alpha}\right\}. \pause \] \begin{itemize} \item $\alpha\rightarrow1$ corresponds to independent variables.\pause \item $\alpha\rightarrow0$ corresponds to perfectly dependent variables. \pause \item This model is symmetric --- the variables are \textit{exchangeable}. \end{itemize} } \frame{ \frametitle{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*} \pause \begin{itemize} \item Independence is obtained when $\alpha=\beta\rightarrow1$ and when one of $\alpha$ or $\beta$ is fixed and the other approaches 1. \pause \item When $\alpha=\beta$ the model reduces to the logistic model. \pause \item The value of $\alpha-\beta$ determines the extent of \textcolor{purple}{\bf asymmetry} in the dependence structure. \end{itemize}} \frame{ \frametitle{6.4.2 Likelihood calculations} After transformation to unit Fr\'{e}chet margins, we can obtain the probability density function of the chosen dependence model by \textcolor{purple}{\bf differentiation} of Equation (6.1) to give $g(x,y)$. \pause \vsp From this, the \textcolor{purple}{\bf likelihood} can be formed (and maximised) in the usual way. \pause } \frame{ \frametitle{6.4.2 Likelihood calculations} 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: \pause %\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.6}{\includegraphics{arrow}} \end{center} %\caption{Threshold classification of bivariate data.} %\end{figure} } \frame{ \frametitle{6.4.2 Likelihood calculations} 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: \pause \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{\sf if}\left(x,y\right)\in\text{\sf Region 1}\\ &\\ \left.\displaystyle\frac{\partial G}{\partial x}\right|_{\left(\tilde{x},\tilde{u}_{y}\right)}\quad&\text{\sf if}\left(x,y\right)\in\text{\sf Region 2}\\ &\\ \left.\displaystyle\frac{\partial G}{\partial y}\right|_{\left(\tilde{u}_{x},\tilde{y}\right)}\quad&\text{\sf if}\left(x,y\right)\in\text{\sf Region 3}\\ &\\ G\left(\tilde{u}_{x},\tilde{u}_{y}\right)\quad&\text{\sf if}\left(x,y\right)\in\text{\sf Region 4} \end{cases} \end{equation*} } \frame{ \frametitle{6.4.2 Likelihood calculations} Then, we have \[ L(\theta;x,y) = \prod_{i=1}^{n}g(\tilde{x}_{i},\tilde{y}_{i}). \] } \frame{ \frametitle{Example 6.1: Wave--surge analysis at Newlyn} Consider the wave--surge data from Newlyn, Cornwall, shown in Figures 6.2 and 6.3. \pause \vsp Flood defences in Newlyn have been designed to withstand a sea swell resulting from, at most, a \textcolor{purple}{\bf wave height} of $x=9 \textrm{\sf metres}$ combined with a \textcolor{purple}{\bf surge height} of $y=0.7 \textrm{\sf metres}$. \vsp A threshold--based approach to modelling is to be used for the wave--surge data shown in Figures 6.2 and 6.3. \vsp Mean residual life plots suggest marginal thresholds of $u_{x}=6.1 \textrm{\sf metres}$ and $u_{y}=0.32 \textrm{\sf metres}$ (as shown in Figure 6.2) for identifying wave height and surge height as extreme. } \frame{ \frametitle{Example 6.1: Wave--surge analysis at Newlyn} \begin{itemize} \item [(a)] Assuming extreme wave heights occur \textbf{independently} of extreme surge heights, find the probability that the flood defence system in Newlyn will be overwhelmed. \end{itemize} } \sol \frame{ \frametitle{Example 6.1: Solution to part (a)} Assuming independence, we have \vsp $\textrm{\sf Pr}(\textrm{\sf flood defence fails})$ \begin{eqnarray*} &=&\textrm{\sf Pr}(X>9) \times \textrm{\sf Pr}(Y>0.7) \\ \pause &=& \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}}\\ \pause & & \\ &=&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} \\ \pause \\ & & \\ &=&0.002995 \times 0.000583 \pause =0.00000175. \end{eqnarray*} } \normal \frame{ \frametitle{Example 6.1: Wave--surge analysis at Newlyn} \begin{itemize} \item [(b)] Now assume there is \textbf{extremal dependence} between wave height and surge height. \begin{itemize} \vspl \item [(i)] Assuming the logistic model for this dependence, obtain the likelihood contributions to $L(\alpha;x_{i},y_{i})$ if $x>u_{x}$ but $y 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} \end{document}