\documentclass[12pt,a4paper]{article} \usepackage{amsmath,fancybox,psfrag,graphicx,pifont} \usepackage{lscape} \usepackage[usenames,dvips]{color} \usepackage{xcolor} \usepackage{listings} %\usepackage{floatflt} %\usepackage{sectsty} %\allsectionsfont{\sffamily} \long\def\symbolfootnote[#1]#2{\begingroup% \def\thefootnote{\fnsymbol{footnote}}\footnote[#1]{#2}\endgroup} \topmargin=-2.0cm \textheight=25.2cm \oddsidemargin=-1.0cm \textwidth=17.5cm \parindent 0cm \parskip 0.4cm \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{gray}, % 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{brown}, % 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 } \newcommand{\cc}[1]{\texttt{#1}} \renewcommand{\vec}[1]{\boldsymbol{#1}} %\renewcommand{\vec}[1]{\mathbf{#1}} %\renewcommand{\vec}[1]{\underline{#1}} \newcommand{\xvec}{\vec{x}} \newcommand{\yvec}{\vec{y}} \newcommand{\Xvec}{\vec{X}} \newcommand{\Yvec}{\vec{Y}} \newcommand{\Evec}{\vec{E}} \newcommand{\tvec}{\vec{t}} \newcommand{\Tvec}{\vec{T}} \newcommand{\mutau}{\begin{pmatrix}\mu\\\tau\end{pmatrix}} \newcommand{\thetavec}{\vec{\theta}} \newcommand{\thetavechat}{\vec{\hat\theta}} \newcommand{\thetahat}{\hat\theta} \newcommand{\xvece}{\vec{x}=(x_1,x_2,\ldots,x_n)^T} \newcommand{\Xvece}{\vec{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$\ }d{\ca}{\ensuremath{\rightarrow}\ } %\title{\textbf{MAS2305/MAS3305 Statistical Inference -- Practical 1}} %\maketitle \begin{center} \section*{MAS3918/8918: Extremes \\ Problems sheet 2: Solutions to assessed questions} \end{center} \begin{itemize} \item [\textbf{1.}] \textcolor{magenta}{\textbf{(15 marks in total)}} \begin{itemize} \item [(a)] Assuming a Gumbel--type tail implies that $\xi=0$ when using the GPD for the set of threshold exceedances. From question 10 in Chapter 7, we know that the log-likelihood is therefore \[ \ell(\sigma; \boldsymbol{y}) = -k\textrm{log}\sigma - \sigma^{-1}\sum_{i=1}^{k}y_{i}; \] we now maximise $\ell(\sigma; \boldsymbol{y})$. Now \[ \frac{\partial \ell}{\partial \sigma} = -k/\sigma +\sum_{i=1}^{k}y_{i} \sigma^{-2}. \qquad \textrm{\textcolor{magenta}{\bf (2 marks)}} \] Equating to zero, and replacing $\sigma$ with $\hat{\sigma}$, gives \begin{eqnarray*} \frac{k}{\hat{\sigma}} &=&\frac{\sum_{i=1}^{k}y_{i}}{\hat{\sigma}^{2}} \\ & & \\ k\hat{\sigma}^{2} &=&\hat{\sigma}\sum_{i=1}^{k}y_{i} \\ \hat{\sigma} &=& k^{-1} \sum_{i=1}^{k}y_{i} = \bar{y}. \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}} \end{eqnarray*} Subtracting $u=40$ from all observations, and then discarding any resulting values which are less than zero, gives the set of threshold exceedances $\boldsymbol{y}$: \begin{center} \begin{tabular}{c c c c c c c} 8 &1 &6 &12 &17 &11 &3 \end{tabular} \end{center} Thus, for the pollution concentration extremes, we have \[ \hat{\sigma} = \frac{8+1+\ldots + 3}{7} = 8.286. \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}} \] \item [(b)] We have \[ \frac{\partial^{2}\ell}{\partial \sigma^{2}} = k\sigma^{-2}-2\sum_{i=1}^{k}y_{i}\sigma^{-3}; \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}} \] evaluating at $\sigma=\hat{\sigma} = 8.286$ gives \[ I_{O} = -\left(\frac{7}{8.286^{2}}-2\times \frac{58}{8.286^{3}}\right) = 0.1019; \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}} \] thus, $\textsf{var}(\hat{\sigma}) = I_{O}^{-1} = 9.8135$, giving $\textrm{e.s.e.}(\hat{\sigma}) = \sqrt{9.8135} = 3.1327$. \textcolor{magenta}{\bf (2 marks)} \newpage \item [(c)] For the $t$--observation return level, we need to solve \[ \text{Pr}(X>x) \approx \lambda_{u} \cdot S(x) = \frac{1}{t} \tag{*} \] for $x=z_{t}$, where $S$ is the survival function ($1-\text{CDF}$) for the Generalised Pareto Distribution. {\color{blue}In the case $\xi\neq 0$ we have \begin{eqnarray*} S(x) &=&1-H(x) \\ &=& 1-\left( 1-\left[1+\xi\left(\frac{x-u}{\sigma}\right)\right]_{+}^{-1/\xi}\right)\\ &=&\left[1+\xi\left(\frac{x-u}{\sigma}\right)\right]_{+}^{-1/\xi}, \end{eqnarray*} as in Eq. (19) (page 36) of the lecture notes; we then solve (*) for $x = \hat{z}_{t}$ after replacing $\sigma$, $\xi$ and $\lambda_{u}$ with their maximum likelihood estimates. This part is not needed for the solution to this question.} For the case $\xi=0$, we replace $S(x)$ with \begin{eqnarray*} S(x) &=&1-H(x) \\ &=& 1-\left(1-\exp\left\{-\frac{x-u}{\sigma}\right\}\right) \\ &=&\exp\left\{-\frac{x-u}{\sigma}\right\} \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}} \end{eqnarray*} as provided on page 31 of the lecture notes (since for the case $\xi=0$ we assume the limit as $\xi \to 0$, which gives the exponential distribution with rate $1/\sigma$). Again, we solve (*) for $x=\hat{z}_{t}$ after replacing $\sigma$ and $\lambda_{u}$ with their maximum likelihood estimates, giving: \begin{eqnarray*} \hat{\lambda}_{u} \cdot \exp\left\{-\frac{\hat{z}_{t}-u}{\hat{\sigma}}\right\} &=& \frac{1}{t} \\ -\frac{(\hat{z}_{t}-u)}{\hat{\sigma}} &=&\log\left(t\hat{\lambda}_{u}\right)^{-1} \\ \hat{z}_{t}-u &=&\hat{\sigma}\log(t\hat{\lambda}_{u}) \\ z_{t} &=&u+\hat{\sigma} \log(t\hat{\lambda}_{u}). \qquad \textrm{\textcolor{magenta}{\bf (2 marks)}} \end{eqnarray*} Now, putting onto an annual scale to get the $r$--\underline{year} return level requires a replacement of $t$ with $r \times n_{y}$ (as we did in the lecture notes for the case $\xi \neq 0$), where $n_{y}$ is the number of observation per year \textcolor{magenta}{\bf (1 mark)}. This gives \[ \hat{z}_{r} = u+\hat{\sigma}\log(rn_{y}\hat{\lambda}_{u}), \] as required. For our data, we have $\hat{\sigma} = 8.286$. Also, \[ \hat{\lambda} = 7/21 = 1/3, \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}} \] and $n_{y} = 365.25$ as we have daily data (thus an average of 365.25 observations per year to account for leap years). This gives \[ \hat{z}_{1} = 40+8.286\textrm{ log}(1\times365.25/3) = 79.789\mu\textrm{g}\textrm{ per }10^{4}\textrm{m}^{3}. \qquad \textrm{\textcolor{magenta}{\bf (2 marks)}} \] \end{itemize} \item [\textbf{2.}] \textcolor{magenta}{\textbf{(10 marks in total)}} Below are outline solutions for one site -- Twickenham. Note that students have been allocated to one of 13 such sites, so answers will vary by student; students must also choose an appropriate threshold after a subjective assessent of a mean residual life plot, meaning that even students allocated to the ame site might have slightly different results. Provided students' solutions look sensible, full marks can be obtained. For reference, indicative results for each site are provided in the table below. \subsubsection*{Twickenham} A mean residual life plot (not shown here) suggests a threshold of around 30mm might be appropriate. Thus: \textcolor{magenta}{\textbf{(1 mark)} for R output} \begin{lstlisting} > gpd.fit(Twickenham, 30) $threshold [1] 30 $nexc [1] 176 $conv [1] 0 $nllh [1] 527.05 $mle [1] 7.69831847 -0.04647912 $rate [1] 0.006601155 $se [1] 0.87720226 0.08552665 \end{lstlisting} \vspace{-0.5cm} Giving $\hat{\sigma} = 7.70 (0.877)$ \textcolor{magenta}{\bf (1 mark)} and $\hat{\xi} = -0.05 (0.086)$ \textcolor{magenta}{\bf (1 mark)}, with standard errors in parentheses. We also have $\hat{\lambda}_{u} = 0.0066 (0.0005)$ \textcolor{magenta}{\bf (1 mark)}, as we know from the notes that \[ \textsf{s.e.}(\hat{\lambda}_{u}) = \sqrt{\hat{\lambda}_{u}(1-\hat{\lambda}_{u})/N }= \sqrt{0.0066(1-0.0066)/26662} = 0.0005. \] For the 100-year return level, we use Eq. (21) in the lecture notes: \begin{eqnarray*} \hat{z}_{100} &=& u+\frac{\hat{\sigma}}{\hat{\xi}}\left[\left(rn_{y}\hat{\lambda}_{u}\right)^{\hat{\xi}}-1\right] \\ &=& 30+\frac{7.70}{-0.05}\left[(100\cdot 365.25\cdot 0.0066)^{-0.05}-1\right] \\ &=& 66.94 \text{mm}. \quad \text{\textcolor{magenta}{\bf (2 marks)}} \end{eqnarray*} For the profile log-likelihood plot, and after some experimentation/trial and improvement: \begin{lstlisting} > gpd.prof(A, 100, xlow=50, xup=100, npy = 365.25, conf = 0.95, nint = 100) > abline(v=86.4, col="red") > abline(v=59.9, col="red") \end{lstlisting} \vspace{-0.5cm} \begin{center} \includegraphics[width=9.5cm]{profplotass2.ps} \end{center} Giving a 95\% confidence interval of approximately $(59.9, 86.4)\text{mm}$. \textcolor{magenta}{\textbf{(2 marks)} for plot, \textbf{(2 marks)} for interval.} \subsubsection*{Indicative results for other sites} \begin{center} \begin{tabular}{|c|c c c|} \hline &$\hat{\sigma}$ &$\hat{\xi}$ &$\hat{z}_{100}$ \\ \hline Twickenham &7.698 &--0.046 &67.274 \\ Bishop Stortford &4.270 &--0.090 &35.079 \\ Boston &3.860 &0.009 &39.260 \\ Bradfield &6.335 &--0.069 &52.550 \\ Cranwell &2.968 &--0.006 &29.046 \\ Hebden Royd &4.887 &0.022 &49.834 \\ Immingham &4.105 &--0.016 &39.534 \\ Lancaster &15.882 &--0.039 &144.896 \\ Leeming &3.242 &0.039 &34.862 \\ Stratford &3.767 &0.009 &38.225 \\ Lyme Regis &3.744 &--0.024 &33.929 \\ Sheffield &7.358 &--0.001 &72.133 \\ Stalybridge &6.075 &0.039 &61.750 \\ \hline \end{tabular} \end{center} \end{itemize} \end{document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% We know that the variance--covariance matrix $V$ is given by \[ V = \left(I_{O}\right)^{-1}, \] where \[ I_{O} = -H. \] Thus, \begin{eqnarray*} V &=& \left(\begin{array}{c c} 2.542674 &-1.591807 \\ -1.591807 &4.486644 \end{array} \right)^{-1} \\ & & \\ &=&\frac{1}{2.542674\times 4.486644 - 1.591807^{2}}\left(\begin{array}{c c} 4.486644 &1.591807 \\ 1.591807 &2.542674 \end{array} \right) \\ & & \\ &=&\left(\begin{array}{c c} 0.5056 &0.1794 \\ 0.1794 &0.2865 \end{array} \right). \qquad \textrm{\textcolor{magenta}{\bf (2 marks)}} \end{eqnarray*} Thus, \begin{eqnarray*} \textrm{e.s.e.}(\hat{\mu}) &=& \sqrt{0.5056} = 0.711 \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}}\\ \textrm{e.s.e.}(\hat{\sigma}) &=& \sqrt{0.2865} = 0.5353. \qquad\textrm{\textcolor{magenta}{\bf (1 mark)}} \end{eqnarray*} So the 95\% confidence intervals for $\mu$ and $\sigma$ are given by \begin{eqnarray*} \mu &:& 45.86 \pm 1.96 \times 0.711 \longrightarrow (44.466,47.254) \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}}\\ \sigma &:& 3.56 \pm 1.96 \times 0.5353 \longrightarrow (2.511, 4.609) \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}} \end{eqnarray*} \item [(b)] \[ \textsf{cor}(\hat{\mu},\hat{\sigma}) = \frac{1.591807}{\sqrt{4.486644 \times 2.542674}} = 0.4713. \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}} \] \item [(c)] For the $r$--year return level $z_{r}$, we equate $G(z_{r};\mu,\sigma)$ to $1-1/r$ and solve for $z_{r}$: \begin{eqnarray*} \textrm{exp}\left\{-\textrm{exp}\left[-\left(\frac{z_{r}-\mu}{\sigma}\right)\right]\right\} &=& 1-\frac{1}{r} \\ && \\ \frac{z_{r}-\mu}{\sigma} &=&-\textrm{log}[-\textrm{log}(1-r^{-1})]\\ & & \\ z_{r} &=&\mu-\sigma\textrm{log}[-\textrm{log}(1-r^{-1})] \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}} \end{eqnarray*} Thus, \begin{eqnarray*} \hat{z}_{50} &=& 45.86-3.56\textrm{log}[-\textrm{log}(1-50^{-1})] = 59.75 \textrm{mph} \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}}\\ \hat{z}_{500}&=& 45.86-3.56\textrm{log}[-\textrm{log}(1-500^{-1})] = 67.98 \textrm{mph} \qquad \textrm{\textcolor{magenta}{\bf (1 mark)}} \end{eqnarray*} Now \textcolor{magenta}{\bf (2 marks) for working -- it's OK if students simply present R code here} \begin{eqnarray*} \frac{\partial z_{r}}{\partial \mu} &=& 1 \qquad \textrm{and} \\ \frac{\partial z_{r}}{\partial \sigma} &=& -\textrm{log}[-\textrm{log}(1-r^{-1})]. \qquad \end{eqnarray*} Thus, \begin{eqnarray*} \textsf{var}(\hat{z}_{r}) &=& \nabla z_{r}^{T} V \nabla z_{r}\\ &=&\left[1, -\textrm{log}[-\textrm{log}(1-r^{-1})]\right] \times V \times \left[\begin{array}{c} 1 \\ -\textrm{log}[-\textrm{log}(1-r^{-1})]\end{array} \right] , \end{eqnarray*} giving \[ \textsf{var}(\hat{z}_{50}) = \left[1,3.901939\right]\times \left(\begin{array}{c c} 0.5056 &0.1794 \\ 0.1794 &0.2865 \end{array} \right) \times \left[ \begin{array}{c} 1 \\ 3.901939 \end{array} \right] = 6.267615 \] and \[ \textsf{var}(\hat{z}_{500}) = \left[1,6.213607\right]\times \left(\begin{array}{c c} 0.5056 &0.1794 \\ 0.1794 &0.2865 \end{array} \right) \times \left[ \begin{array}{c} 1 \\ 6.213607 \end{array} \right] = 13.79650. \] So we have \begin{eqnarray*} \textrm{e.s.e.}(\hat{z}_{50}) &=&\sqrt{6.267615} = 2.5035 \quad \textrm{and} \qquad\textrm{\textcolor{magenta}{\bf (1 mark)}}\\ \textrm{e.s.e.}(\hat{z}_{500}) &=&\sqrt{13.79650} = 3.7144.\qquad\textrm{\textcolor{magenta}{\bf (1 mark)}} \end{eqnarray*} We would advise against using these standard errors to construct symmetric confidence intervals for the associated return levelsowing to the (severe) asymmetry often observed in the likelihood surface for return levels. \textrm{\textcolor{magenta}{\bf (1 mark)}} \end{itemize} \item [\textbf{2.}] \textcolor{magenta}{\textbf{(10 marks in total)}} All students, regardless of their dataset allocation, should present the following: \begin{enumerate} \item Data summaries -- simple histograms and/or time series plots of their data \textcolor{magenta}{\textbf{(1 mark)}} \item Fitting of the GEV: {\tiny \begin{center} \begin{tabular}{|l|l |l|l|c|l|} \hline &$\hat{\mu}$ &$\hat{\sigma}$ &$\hat{\xi}$ &$\hat{z}_{500}$ &Tail? \\ \hline Temperature &$32.96 (0.34)\longrightarrow (32.29,33.62)$ &$1.98 (0.25) \longrightarrow (1.49, 2.46)$ &$0.02 (0.12) \longrightarrow (-0.26, 0.21)$ &$44.37 (3.83)$ &Gumbel \\ Rainfall &$29.23 (1.10)\longrightarrow (27.08, 31.38)$ &$6.16 (0.91) \longrightarrow (4.38, 7.95)$ &$0.27 (0.16) \longrightarrow (-0.04, 0.57)$ &$127.51 (58.96)$ &Gumbel \\ Wind &$8.94 (0.18)\longrightarrow (8.59, 9.29)$ &$1.00 (0.14) \longrightarrow (0.73, 1.28)$ &$-0.34 (0.16)\longrightarrow (-0.65, -0.02)$ &$11.56 (0.58)$ &Weibull \\ Sunshine &$14.80 (0.06) \longrightarrow (14.7, 14.93)$ &$0.40 (0.05) \longrightarrow (0.30, 0.50)$ &$-0.65 (0.08) \longrightarrow (-0.82, -0.49)$ &$15.40 (0.02)$ & Weibull \\ \hline \end{tabular} \end{center} } Student should report: \begin{itemize} \item Estimates of $\mu$, $\sigma$ and $\xi$ with standard errors and/or 95\% confidence intervals\\ \textcolor{magenta}{\textbf{(3 marks)}, deduct 2 marks if standard errors/confidence intervals omitted.} \item Their estimate of the 500-year return level \textcolor{magenta}{\textbf{(1 mark)}} with standard error \textcolor{magenta}{\textbf{(2 marks)}} \item Comment on the tail behaviour \textcolor{magenta}{\textbf{(1 mark)}} \end{itemize} \item Diagnostic plots -- probability plots, quantile plots etc \textcolor{magenta}{\textbf{(2 marks)}, deduct 1 mark if no comments provided.} \end{enumerate} \end{itemize} \end{document}