\documentclass[12pt]{article}
\textwidth=17.5cm
\addtolength{\textheight}{2.5in}
\topmargin=-2cm
\addtolength{\hoffset}{-2cm}
\usepackage{epsfig,amsmath,latexsym,amsfonts,enumerate,fancybox,psfrag,color}

\pagestyle{empty}

\newcommand{\announce}{\setbeamercolor{structure}{fg=green!60!black}}
\newcommand{\sol}{\setbeamercolor{structure}{fg=purple!80!black}}
\newcommand{\normal}{\setbeamercolor{structure}{fg=blue!60!green}}
\newcommand{\vsp}{\vspace{0.3cm}}
\newcommand{\vspl}{\vspace{0.15cm}}
\newcommand{\xvec}{\boldsymbol{x}}
\newcommand{\yvec}{\vec{y}}
\newcommand{\Xvec}{\boldsymbol{X}}
\newcommand{\Yvec}{\vec{Y}}
\newcommand{\Evec}{\vec{E}}
\newcommand{\tvec}{\bodsymbol{t}}
\newcommand{\Tvec}{\bodsymbol{T}}
\newcommand{\mutau}{\begin{pmatrix}\mu\\\tau\end{pmatrix}}
\newcommand{\thetavec}{\vec{\theta}}
\newcommand{\thetavechat}{\vec{\hat\theta}}
\newcommand{\thetahat}{\hat\theta}
\newcommand{\xvece}{\boldsymbol{x}=(x_1,x_2,\ldots,x_n)^T}
\newcommand{\Xvece}{\boldsymbol{X}=(X_1,X_2,\ldots,X_n)^T}
\newcommand{\thetavece}{\vec{\theta}=(\theta_1,\theta_2,\ldots,\theta_p)^T}
\newcommand{\exxvec}{$x_1,x_2,\ldots,x_n$\xspace}
\newcommand{\exxvecdot}{$x_1,x_2,\ldots,x_n$.\ }
\newcommand{\exXvec}{$X_1,X_2,\ldots,X_n$\ }

\begin{document}
% \enlargethispage{2cm}
%\begin{table}[h!]
%\begin{tabular}{|c|}
%\hline
%\hspace{0.25cm} Date Stamp \hspace{0.25cm} \\
%\\
%\\
%\hline
%\end{tabular}
%\end{table}
%
%\vspace{-2.5cm}

\begin{center}
\Large
\shadowbox{\sf \textbf{MAS8306: Problems sheet 1}} 

\vspace{0.3cm}

\Large
{\sf \textbf{Estimating Environmental Extremes}}
\end{center}
%\vspace{-0.5cm}
%
%{\flushleft{
%\fbox{Name: \hspace{8cm} Student Number: \hspace{5.6cm}} \\
%
%\vspace{0.1cm}
%
%\fbox{Stage \hspace{6cm} Lecturer: \hspace{5cm} Date: \hspace{3cm}}
%}} \\
%

\vspace{0.5cm}

\normalsize
%\noindent \textsc{Problems Sheet 1} is to be used in next week's problems class (4pm, Tues 14th Feb).  Further help will be available in the drop--in session the week after (4pm, Thurs 21st Feb).      

\vspace{0.5cm}

\noindent Answers to all starred questions should be placed in the homework submission box in the foyer of the Maths \& Stats General Office (Herschel Building 3rd floor) by no later than \textcolor{red}{\bf 4pm, Thursday 22nd February}.  

 %\textit{Don't forget that CBA1 should be completed by then as well!}  %CBA1 goes live in practice mode Monday 9th February and in exam mode Monday 16th February.}      

\vspace{1cm}

\begin{itemize}

\item [\textbf{1.}]  Given the cumulative distribution function for the generalised extreme value (GEV) distribution, as provided by Equations (2.6) and (2.7) in the lecture notes (for $\xi\neq 0$ and $\xi=0$, respectively), obtain the corresponding probability density function for the GEV.  For the case $\xi=0$, show that the log--likelihood is given by
\[
\ell(\mu,\sigma; \boldsymbol{x}) = -m\textrm{log}\sigma - \sum_{i=1}^{m}\left\{\frac{x_{i}-\mu}{\sigma}\right\}-\sum_{i=1}^{m}\textrm{exp}\left\{-\left(\frac{x_{i}-\mu}{\sigma}\right)\right\},
\]
where $m$ is the number of block maxima.  

\item [\textbf{2*.}] A type I extreme value distribution, with distribution function given by
\[
G(x; \mu, \sigma) = \textrm{exp}\left\{-\textrm{exp}\left[-\left(\frac{x-\mu}{\sigma}\right)\right]\right\}
\]
for $-\infty < \mu < \infty$ and $\sigma>0$, is used to model annual maximum wind speeds (mph) at Corpus Christi, Texas, for the years 1943 and 1976 (inclusive).  The maximum likelihood estimates for $\mu$ and $\sigma$ are found to be $\hat{\mu} = 45.86$ and $\hat{\sigma} = 3.56$; the hessian matrix for the log--likelihood is also found as
\[
H = \left( \begin{array}{c c} 
          -2.542674 &1.591807 \\
          1.591807  &-4.486644
\end{array} \right).
\]
\begin{itemize}
\item [(a)] Estimate the standard errors for $\hat{\mu}$ and $\hat{\sigma}$, and use these to construct 95\% confidence intervals in the usual way.
\item [(b)] Find $\textsf{cor}(\hat{\mu},\hat{\sigma})$.  
\item [(c)] Obtain an expression for the $r$--year return level using this type I extreme value distribution.  Thus, estimate the 50--year and 500--year return level wind speeds at Corpus Christi, and obtain standard errors for these estimates via the delta method.  Why is it not advisable to construct confidence intervals for return level estimates in the usual way?
\end{itemize}

\item[\textbf{3*.}] Suppose $X_{1}, X_{2} \ldots, X_{n}$ is a sequence of independent Gumbel random variables.  By letting $a_{n}=1$ and $b_{n}=\textrm{log }n$, show that the distribution of $(M_{n}-b_{n})/a_{n}$ is of extreme value type, where
\[
M_{n} = \textrm{max}\left\{X_{1}, X_{2}, \ldots, X_{n}\right\}.
\]



\end{itemize}

\newpage
\LARGE
\noindent 
\shadowbox{\sf \textbf{Practical work}} 
\normalsize

\vspace{1cm}
\noindent Open \textbf{\sf R} or \textbf{\sf Rstudio}.  You will also need to access data for some of these questions from the course website, so have this open:
\begin{center}
\texttt{www.mas.ncl.ac.uk/$\sim$nlf8/teaching/mas8306}
\end{center}

%\noindent \textbf{No solution to question 4 is required in your submission}.  
\begin{itemize}
\item [\textbf{4.}] Recall the first step in a typical analysis of extremes: your data may need to pre--processed to extract the block maxima to which the generalised extreme value distribution can be fitted.  In the sea--surge example used throughout Chapter 2 of the lecture notes, no data pre--processing was necessary, as we already had the set of annual maximum sea surges for the years 1955--2004 (inclusive).  
\begin{itemize}
\item [(a)] Suppose you have daily data collected over a period of seven years -- from 2003 to 2009 (inclusive) -- and you want to extract the seven annual maxima from this series.  We can do this in \textbf{\sf R}.  First of all, set up a vector of years using the code:
\begin{verbatim}
> year=seq(2003,2009,1)
\end{verbatim}
The \textbf{\sf R} function \texttt{seq($a$,$b$,$c$)} creates a sequence of values from $a$ up to $b$ in steps of $c$.  Take a look at your vector \texttt{year} to make sure you get the desired result, by typing:
\begin{verbatim}
> year
\end{verbatim}
\item [(b)] If we use calendar years as our blocks, we should ensure that we take into account leap years.  The following code creates a vector the same length as \texttt{year} called \texttt{ind}; this will indicate whether or not the corresponding year is a leap year by taking the value 366 for leap years, and 365 elsewhere (a leap year has 366 days, a non--leap year 365).  Look through the code, and make sure you fully understand each line, before applying this in \textbf{\sf R}. 
\begin{verbatim}
> year.over.four=year/4
> rounded=round(year.over.four)
> ind=vector("numeric",length(year))
> for(i in 1:length(ind))
  {
    if(rounded[i]==year.over.four[i])
      {
        ind[i]=366
      }
    else{
      ind[i]=365
        }
   }
\end{verbatim}
\item [(c)] When obtaining the set of annual maxima, we want to find the maximum value between \texttt{lower} and \texttt{upper} in each block: for example, for our first year of data, \texttt{lower} would be 1 and \texttt{upper} 365 (since 2003 is not a leap year), for the second year of data, \texttt{lower} would be 366 and \texttt{upper} 731 (since 2004 \textit{is} a leap year).  

\newpage

We can compute \texttt{lower} and \texttt{upper} for each year using the following code -- read through this carefully, and then type it into \textbf{\sf R}.  
\begin{verbatim}
> upper=vector("numeric",length(ind))
> for(i in 1:length(upper))
  {
     upper[i]=sum(ind[1:i])
   }
> lower=vector("numeric",length(ind))
> lower[1]=1
> for(i in 2:length(lower))
  {
    lower[i]=upper[i-1]+1
  }
\end{verbatim}
\item [(d)] Now typing
\begin{verbatim}
> cbind(year,ind,lower,upper)
\end{verbatim}
gives
\begin{verbatim}
     year ind lower upper
[1,] 2003 365     1   365
[2,] 2004 366   366   731
[3,] 2005 365   732  1096
[4,] 2006 365  1097  1461
[5,] 2007 365  1462  1826
[6,] 2008 366  1827  2192
[7,] 2009 365  2193  2557
\end{verbatim}
\item [(e)] Suppose our daily measurements observed over the seven years are stored in the vector \texttt{mydata}; then the following code would extract the set of annual maxima:
\begin{verbatim}
> anmax=vector("numeric",length(year))
> for(i in 1:length(anmax))
   {
    anmax[i]=max(mydata[lower[i]:upper[i]])
   }
\end{verbatim}
Again, read through this code carefully and make sure you fully understand what will happen in the \texttt{for} loop.  \textit{Do not attempt to execute this part of the code, as you won't have anything stored in }\texttt{mydata}\textit{!}
\item [(f)] How would you amend the above procedure, as outlined in parts (a)--(e), if you had \textit{hourly} observations?
\end{itemize}

\newpage

\item [\textbf{5*.}] Daily rainfall accumulations have been recorded at 25 sites in England over a period of 74 years, from 1938 to 2011 inclusive.  You have been randomly allocated to one of these sites.  Go to the assignments section of the course website and click on the ``\textsf{rainfall data}'' link next to the material for \textsf{Problems sheet 1}; then right--click on your name and \texttt{Save Target As} ``\textsf{rainfall.txt''}, making sure you save your data file in a suitable folder on your \texttt{H:} drive.    
\\
\\
Now you can load the rainfall data into \textbf{\sf R} by typing
\begin{verbatim}
> rain=scan('rainfall.txt')
\end{verbatim}
This should read the rainfall data into \textbf{\sf R} and store them in the vector \texttt{rain}. If you get the following message:
\begin{verbatim}
  Cannot open file 'rainfall.txt': No such file or directory
\end{verbatim}
then you do not have \textbf{\sf R} open in the same directory as you saved your data file.  To change to the correct directory, click \texttt{File}$\rightarrow$\texttt{Change dir...} and then browse through your folders until you find the folder in which you saved your data file -- select this folder, and then click \texttt{OK}.  Repeating the \texttt{scan} code above should now successfully load your rainfall data into \textbf{\sf R}.  

\begin{itemize}
\item [(a)] By hand/using a calculator, work out how many observations there \textit{should} be in \texttt{rain}.  Now use the \textbf{\sf R} function \texttt{length}, or otherwise, to find out how many observations you \textit{actually} have.  What do you notice?
\item [(b)] Actually, data is missing for the year 1944 during the second world war.  How does this explain your answer to part (a)?
\item [(c)] Following the example code given in question 4(a), set up a \texttt{year} vector in \textbf{\sf R} for your rainfall data -- being careful to exclude the year 1944.  Now follow the example \textbf{\sf R} code shown in the rest of question 4 to extract the set of annual rainfall maxima for your site, storing your set of maxima in the vector \texttt{anmax}.  

Use the \textbf{\sf R} command \texttt{summary} to obtain the five number summary for your annual maximum rainfall accumulations (min, lower quartile, median, upper quartile, max), and include these in your solutions.     
\item [(d)] Produce a time series plot by plotting \texttt{anmax} against \texttt{year}, and include this plot in your solutions (using appropriate labels for the title/axes). Also produce a histogram of your annual maxima. Comment.
\end{itemize}

\vspace{0.5cm}

\item [\textbf{6*.}] This question will make use of the annual maximum rainfall accumulations you extracted from your rainfall dataset in question 5.  
\begin{itemize}
\item [(a)] Following the material in Section 2.2.2 of the lecture notes, use the \texttt{nlm} function in \textbf{\sf R} to fit the GEV to the set of annual maximum rainfall accumulations at your site.  In your solution, include the output \textbf{\sf R} produces when you apply the \texttt{nlm} routine, including the \texttt{minimum} of the function, the \texttt{estimates}, the convergence \texttt{code} and the number of \texttt{iterations} (see page 23 of the lecture notes).  What is the value of your maximised GEV log--likelihood?  
\item [(b)] Again, following the example in Section 2.2.2 of the lecture notes, use \textbf{\sf R} to invert the hessian matrix and thus obtain the variance--covariance matrix $V$.  In your solutions, write down $V$, and also write down your corresponding estimated standard errors for $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\xi}$.  Use these standard errors to construct 95\% confidence intervals for the GEV parameters in the usual way, and comment.  
\item [(c)] Estimate the 10--year, 100--year and 1000--year return levels for annual maximum rainfall accumulations at your site.  Also, use the delta method to estimate the standard errors for your return level estimates, following the example in Section 2.2.4 of the lecture notes.  It is not necessary to include any \textbf{\sf R} code here; just produce a table summarising your results similar to Table 2.2 in the lecture notes.  
\item [(d)] Use the function \texttt{gev.fit} within the \textsf{ismev} package in \textbf{\sf R} to confirm your answers to parts (a) and (b) -- include in your solutions the \textbf{\sf R} output given on execution of \texttt{gev.fit}.  What are the differences between the \texttt{gev.fit} fitting procedure and fitting the GEV ``from first principles'' using \texttt{nlm}?  
\end{itemize}

\vspace{0.5cm}

\item [\textbf{7.}] The \textbf{\sf R} package \textsf{evd} is another add--on package for extreme value analyses.  Load the package into \textbf{\sf R} by typing
\begin{verbatim}
> library(evd)
\end{verbatim}
We will use many of the functions provided by \texttt{evd} in the next problems sheet.  For now, one interesting suite of functions is that which calculates the GEV probability density (\texttt{dgev}), cumulative distribution function (\texttt{pgev}), quantile function (\texttt{qgev}) and allows random generation of extremes from a GEV (\texttt{rgev}).  Explore this suite of functions by typing
\begin{verbatim}
> ?pgev
\end{verbatim}
These functions work in exactly the same way as the analogous functions for the Normal distribution or the Poisson distribution -- for example, \texttt{dnorm}, \texttt{rnorm} or \texttt{ppois}.  
\begin{itemize}
\item [(a)] Generate a sequence of values between 0 and 50, in steps of 0.1, using the command \texttt{seq} (see question 4(a) for a reminder).  Store this sequence in \texttt{x}.   
\item [(b)] Now use the command \texttt{dgev} to find the GEV probability density for each value in \texttt{x}, using $\mu=20$, $\sigma=3.5$ and $\xi=0.3$ by typing
\begin{verbatim}
> y1=dgev(x,loc=20,scale=3.5,shape=0.3)
\end{verbatim}
\item [(c)] Produce a density plot for the GEV $\mathcal{G}(\mu=20,\sigma=3.5,\xi=0.3)$ distribution, by plotting \texttt{y1} against \texttt{x}.  Use \texttt{plot}, and join the points by using \texttt{type=l}.  
\item [(d)] Following parts (b) and (c), also produce a plot of the GEV densities when $\xi=0$ and $\xi<0$.  Superimpose these densities on the plot you produced in part (c) by using the \texttt{lines} command (instead of \texttt{plot}), and use different colours to distinguish between the three plots by using the argument \texttt{color='blue'}, for example. 
\item [(e)] What do you notice when you compare your three density plots? 
\item [(f)] Now experiment with three different values of $\mu$, keeping both $\sigma$ and $\xi$ constant.  
\item [(g)] Now experiment with three different values of $\sigma$, keeping both $\mu$ and $\xi$ constant.  
\end{itemize}
\end{itemize}


\end{document}