\chapter{Threshold methods} \section{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 ({\it threshold}) are considered. This allows more efficient use of data and avoids the problems that can arise as a result of blocking (see Section 2.6), but brings its own problems. We must first go back and consider the asymptotic theory appropriate for this new situation. \\ \\ 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{max}\left\{X_{1}, \ldots, X_{n}\right\}. \] We know from Chapter 2 that \[ \textrm{Pr}\left\{M_{n}\leq x\right\} \approx G(x), \] where \[ G(x) = \textrm{exp}\left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]_{+}^{-1/\xi}\right\} \] is the Generalised Extreme Value distribution with location, scale and shape $\mu$, $\sigma$ and $\xi$ respectively. \vspace{0.8cm} \noindent {\large \textbf{Theorem} \textit{(Distribution of threshold excess)}} For a large enough threshold $u$, the distribution function of $(X-u)$, conditional on $X>u$, is approximately \begin{eqnarray} H(y) = 1-\left(1+\frac{\xi y}{\tilde{\sigma}}\right)_{+}^{-1/\xi}, \end{eqnarray} defined on $y>0$, where \begin{eqnarray} \tilde{\sigma} = \sigma +\xi(u-\mu). \end{eqnarray} \newpage \vspace{0.6cm} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} {\large \textbf{Proof}} \huge \ding{46} \normalsize \vspace{22cm} }} \newpage \vspace{0.6cm} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} \huge \ding{46} {\large \textit{...proof continued...}} \normalsize \vspace{19cm} }} \vspace{0.8cm} The family of distributions defined by Equation (3.1) is known as the \textit{Generalised Pareto family}; the distribution itself is often referred to as the Generalised Pareto Distribution, or GPD for short. \subsection*{Comments} \begin{itemize} \item If block maxima have approximate distribution $G$, then threshold excesses have a corresponding distribution given by the Generalised Pareto family. \item The parameters of the GPD are uniquely determined by those of the GEV: \begin{itemize} \item [--] The parameter $\xi$ in Equation (3.1) is equal to that of the corresponding GEV; \item [--] The GPD scale parameter is a function of the GEV location and shape parameters. \end{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'. \item The duality between the GEV and GPD means that the shape parameter $\xi$ is dominant in determining the qualitative behaviour of the GPD: \begin{itemize} \item If $\xi<0$ the distribution of excesses has an upper bound; \item if $\xi>0$ the distribution has no upper limit; \item the case $\xi=0$ is also unbounded, and is taken as the limit $\xi \to 0$, giving \[ H(y) = 1 -\textrm{exp}\left(-\frac{y}{\tilde{\sigma}}\right), \qquad y>0; \] i.e. an exponential distribution with rate $1/\tilde{\sigma}$. \end{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. For notational convenience we now drop this distinction, using $\sigma$ to denote the scale parameter within either family. \end{itemize} \newpage %\vspace{0.6cm} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} {\large \textbf{Example 3.1}} \vspace{0.4cm} Suppose $X_{1}, X_{2}, \ldots, X_{n}$ is a sequence of independent \textit{exp}(1) random variables. Show that the limiting distribution of threshold excesses belongs to the generalised Pareto family. \huge \ding{46} \normalsize \vspace{8cm} }} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} {\large \textbf{Example 3.2}} \vspace{0.4cm} Suppose $X_{1}, X_{2}, \ldots, X_{n}$ is a sequence of independent $U(0,1)$ random variables. Show that the limiting distribution of threshold excesses belongs to the generalised Pareto family. \huge \ding{46} \normalsize \vspace{8cm} }} \section{Application: rainfall extremes in New York} 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). 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, ``Superstorm Sandy'' in 2012). Thus, analysing extreme rainfall data has a real practical motivation here, in terms of river and sea flood defence systems. In this Section, we will illustrate a complete threshold--based analysis of the rainfall extremes observed at New York. \\ \\ The data have been scanned into \textbf{\sf R} from the course webpage and stored in the vector \texttt{rain}: \begin{verbatim} > rain=scan('newyork.txt') \end{verbatim} \subsection{Threshold choice} \noindent The \textit{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}$. 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 \begin{eqnarray} E[X-u|X>u] &=& \frac{\sigma_{u_{0}}+\xi u}{1-\xi}. \end{eqnarray} Thus, for all $u>u_{0}$, $E[X-u|X>u]$, is a linear function of $u$. 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. This leads to the \textit{mean residual life plot}, a graphical procedure for identifying a suitably high threshold for modelling extremes via the GPD. 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. We can easily do this from first principles in \textbf{\sf R}. First of all, we set up a vector of possible thresholds, starting at zero and going up to the maximum value in our dataset: \begin{verbatim} > u=seq(0,max(rain),0.1) \end{verbatim} The vector \texttt{x} is then set up to take the corresponding values for the mean excess over each value in \texttt{u}: \begin{verbatim} > x=vector('numeric', length(u)) \end{verbatim} Then the following code computes the mean excess for each value in \texttt{u} and stores it in \texttt{x}: \newpage \begin{verbatim} > for(i in 1:length(x)) { threshold.exceedances=rain[rain>u[i]] x[i]=mean(threshold.exceedances-u[i]) } \end{verbatim} The MRL plot is then produced using the following code, giving the plot in Figure 3.1: \begin{verbatim} > plot(x~u,type='l', main='MRL plot',ylab='mean excess') \end{verbatim} \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=15cm,height=10cm]{figure10.ps} \caption{Mean residual life plot for the New York rainfall data.} \end{figure} Though interpretation of these plots can be subjective, linearity in Figure 3.1 might be suggested above $u_{0}=30\textrm{mm}$ (information in the far right--hand--side of these plots is unreliable; here, variability is high due to the limited amount of data above such high thresholds). Of course, there is a function in the \textsf{ismev} package that produces exactly the same plot, but including 95\% confidence bounds for the mean excess: \begin{verbatim} > mrl.plot(rain) \end{verbatim} \subsection{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. 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. In a threshold--based analysis the data pre--processing is far more straightforward. 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: \begin{verbatim} > above.threshold=rain[rain>30] > threshold.exceedances=above.threshold-30 \end{verbatim} We can look at our set of threshold exceedances by typing: \begin{verbatim} > threshold.exceedances [1] 1.8 2.5 1.8 14.5 0.5 13.2 5.6 8.1 2.0 1.8 3.0 9.1 0.5 1.8 2.3 [16] 3.0 0.5 2.5 18.5 5.3 10.6 0.5 4.3 2.8 0.5 15.7 1.8 3.5 3.5 1.8 [31] 4.8 5.3 7.8 46.7 2.3 4.0 3.8 6.6 0.5 15.7 56.6 5.6 17.8 17.5 4.3 [46] 18.5 0.7 13.4 29.4 5.1 23.3 3.5 0.5 0.2 10.9 12.7 53.3 24.9 29.2 1.8 [61] 7.3 2.5 4.0 37.3 1.2 0.2 6.1 6.8 8.4 1.0 3.3 17.0 2.0 3.0 8.1 [76] 0.5 42.4 4.3 7.1 3.0 10.9 9.9 17.0 6.3 0.5 0.5 25.9 1.8 21.3 55.3 [91] 11.9 0.5 3.0 5.6 25.9 14.2 8.1 4.3 1.8 2.0 1.8 5.6 15.2 0.5 9.4 [106] 0.2 14.5 1.8 3.8 21.6 5.3 29.4 3.5 5.3 0.5 6.8 17.8 12.9 7.6 25.4 [121] 5.3 12.4 3.0 3.0 10.1 4.8 8.1 9.4 4.0 5.6 4.3 3.5 1.0 6.6 6.3 [136] 8.4 8.1 17.0 1.0 0.5 1.2 5.6 18.8 11.9 1.7 1.2 21.3 3.5 7.6 9.4 [151] 9.4 15.7 \end{verbatim} Thus, we have identified 152 observations as being extreme. Compare this to an analysis of annual maxima in which we would have only 48 observations to work with -- one from each year. \subsection{Fitting the GPD} The most commonly--used approach to fit the GPD to the set of threshold excesses is that of maximum likelihood. The GPD log--likelihood function can be derived in the usual way; this is left as an exercise (actually, you are asked to do this in problems sheet 2!), but can be shown to be: \begin{eqnarray} \ell(\sigma,\xi; \boldsymbol{y}) &=& -152\textrm{log}\sigma - (1+1/\xi)\sum_{i=1}^{152}\textrm{log}\left(1+\frac{\xi y_{i}}{\sigma}\right)_{+}, \end{eqnarray} where $\boldsymbol{y}=(y_{1}, \ldots, y_{152})$ are the set of exceedances above threshold $u_{0}=30$. For the case $\xi=0$, interpreted as $\xi \to 0$, we have the log--likelihood for an exponential distribution with rate $1/\sigma$ (again, see problems sheet 2). In \textbf{\sf R}, we could now proceed as we did in Section 2.2.2 when fitting the GEV; that is, we could write a function which computes the (negative) log--likelihood for the GPD (using the 152 threshold excesses), and then use the \texttt{nlm} routine to minimise this with respect to $\theta=(\sigma,\xi)$. Actually, this is exactly what you will have to do in problems sheet 2! Here, we will obtain maximum likelihood estimates of the GPD parameters for our rainfall extremes by using the \texttt{gpd.fit} function in \textsf{ismev}: \begin{verbatim} > gpd.fit(rain,30) $threshold [1] 30 $nexc [1] 152 $conv [1] 0 $nllh [1] 485.0937 $mle [1] 7.4404389 0.1843305 $rate [1] 0.008669861 $se [1] 0.9583796 0.1011461 \end{verbatim} Thus, our fitted model is given by \[ \hat{\sigma}=7.440 (0.958) \qquad \hat{\xi}=0.184 (0.101) \] Note also that the output gives the number of exceedances (\texttt{nexc} = 152), the threshold exceedance rate (\texttt{rate} = 0.0087; just the number of exceedances divided by the total length of the series, that is, 152/17532), and the value of the minimised negative log--likelihood (\texttt{nllh}; thus, the maximised log--likelihood is --485.0937). \subsection{Model adequacy} We can use probability plots and quantile plots to check the suitability of the fitted GPD to the set of extracted threshold exceedances; as demonstrated in Chapter 2 (Section 2.2.3), such plots are easily done from first principles in \textbf{\sf R}; however, the function \texttt{gpd.diag} in \textsf{ismev} does for the GPD exactly what \texttt{gev.diag} did for the GEV in Chapter 2 (Section 2.2.5): \begin{verbatim} > A=gpd.fit(rain,30) > gpd.diag(A) \end{verbatim} The resulting plots are shown overleaf in Figure 3.2; all diagnostics seem to indicate a reasonable fit of the GPD to our rainfall extremes. \newpage \begin{figure} \centering \psfrag{0.0}[][]{\tiny \sf 0.0} \psfrag{0.2}[][]{\tiny \sf 0.2} \psfrag{0.4}[][]{\tiny \sf 0.4} \psfrag{0.6}[][]{\tiny \sf 0.6} \psfrag{0.8}[][]{\tiny \sf 0.8} \psfrag{1.0}[][]{\tiny \sf 1.0} \psfrag{7}[][]{\tiny \sf 7} \psfrag{8}[][]{\tiny \sf 8} \psfrag{9}[][]{\tiny \sf 9} \psfrag{10}[][]{\tiny \sf 10} \psfrag{11}[][]{\tiny \sf 11} \psfrag{12}[][]{\tiny \sf 12} \psfrag{13}[][]{\tiny \sf 13} \psfrag{14}[][]{\tiny \sf 14} \psfrag{0.00}[][]{\tiny \sf 0.00} \psfrag{0.02}[][]{\tiny \sf } \psfrag{0.04}[][]{\tiny \sf 0.04} \psfrag{0.06}[][]{\tiny \sf } \psfrag{0.08}[][]{\tiny \sf 0.08} \psfrag{0.10}[][]{\tiny \sf } \psfrag{0.05}[][]{\tiny \sf 0.05} \psfrag{0.15}[][]{\tiny \sf 0.15} \psfrag{0.20}[][]{\tiny \sf 0.20} \psfrag{0.25}[][]{\tiny \sf 0.25} \psfrag{30}[][]{\tiny\sf 30} \psfrag{40}[][]{\tiny\sf 40} \psfrag{50}[][]{\tiny\sf 50} \psfrag{60}[][]{\tiny\sf 60} \psfrag{70}[][]{\tiny\sf 70} \psfrag{80}[][]{\tiny\sf 80} \psfrag{100}[][]{\tiny\sf 100} \psfrag{200}[][]{\tiny\sf 200} \psfrag{300}[][]{\tiny\sf 300} \psfrag{400}[][]{\tiny\sf 400} \psfrag{x}[][]{\tiny $x$} \psfrag{1e-01}[][]{\tiny \sf 0.1} \psfrag{1e+00}[][]{\tiny \sf 1} \psfrag{1e+01}[][]{\tiny \sf 10} \psfrag{1e+02}[][]{\tiny \sf 100} \psfrag{1e+03}[][]{\tiny \sf 1000} \psfrag{Model}[][]{\tiny \sf Model} \psfrag{Empirical}[][]{\tiny \sf Empirical} \psfrag{z}[][]{\tiny \sf z} \psfrag{f}[][]{\tiny \sf f(x)} \psfrag{Return Level}[][]{\tiny \sf Return level} \psfrag{Return Period}[][]{\tiny \sf Return period} \psfrag{Probability Plot}[][]{\scriptsize \sf Probability plot} \psfrag{Quantile Plot}[][]{\scriptsize \sf Quantile plot} \psfrag{Return Level Plot}[][]{\scriptsize \sf Return level plot} \psfrag{Density Plot}[][]{\scriptsize \sf Density plot} \includegraphics[width=14.5cm,height=6.5cm]{newyorkgpd.ps} \caption{Diagnostic plots indicating the goodness--of--fit of the GPD to the New York rainfall extremes.} \end{figure} \subsection{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, \begin{eqnarray} \textrm{Pr}(X>u+y|X>u) \approx \left[1+\frac{\hat{\xi} y}{\hat{\sigma}}\right]_{+}^{-1/\hat{\xi}}, \end{eqnarray} for $\xi \neq 0$. Working with the left--hand--side of (3.5), we see that \[ \textrm{Pr}(X>u+y|X>u) =\frac{\textrm{Pr}(X>u+y)}{\textrm{Pr}(X>u)}, \] giving \[ \textrm{Pr}(X>u+y) = \textrm{Pr}(X>u)\textrm{Pr}(X>u+y|X>u). \] After substitution of (3.5), we get \begin{eqnarray} \textrm{Pr}(X>u+y) &\approx& \hat{\lambda}_{u}\left[1+\frac{\hat{\xi} y}{\hat{\sigma}}\right]_{+}^{-1/\hat{\xi}}, \end{eqnarray} where $\lambda_{u}=\textrm{Pr}(X>u)$ and is estimated as the empirical threshold exceedance rate $\hat{\lambda}_{u}$. 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$; substitution of $y_{i}=x_{i}-u$ into (3.6) gives \begin{eqnarray} \textrm{Pr}(X>x) &\approx &\hat{\lambda}_{u}\left[1+\hat{\xi}\left(\frac{x-u}{\hat{\sigma}}\right)\right]_{+}^{-1/\hat{\xi}}. \end{eqnarray} Thus, an estimate of the level $z_{t}$ that is exceeded on average once every $t$ 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*} giving \begin{eqnarray} \hat{z}_{t} &=& u +\frac{\hat{\sigma}}{\hat{\xi}}\left[(t\hat{\lambda}_{u})^{\hat{\xi}}-1\right] \end{eqnarray} for $\xi \neq 0$, and \[ \hat{z}_{t} = u +\hat{\sigma}\textrm{log}(t \hat{\lambda}_{u}) \] when $\xi=0$ (see problems sheet 2). By construction, $z_{t}$ is the $t$--\underline{observation} return level; however, it is often more convenient to give return levels on an annual scale, so that the $r$--\underline{year} return level is the level expected to be exceeded once every $r$ \underline{years}. If there are $n_{y}$ observations per year, this corresponds to the $t$--observation return level with $t=r\times n_{y}$. Hence, an estimate of the $r$--year return level $z_{r}$ is defined by \begin{eqnarray} \hat{z}_{r}&=&u+\frac{\hat{\sigma}}{\hat{\xi}}\left[(rn_{y}\hat{\lambda}_{u})^{\hat{\xi}}-1\right], \end{eqnarray} unless $\xi=0$, in which case \begin{eqnarray*} \hat{z}_{r}&=&u+\hat{\sigma}\textrm{log}(rn_{y}\hat{\lambda}_{u}). \end{eqnarray*} 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{mm} \] as an estimate of the 50--year return level, where $n_{y}=365.25$ to account for leap years (we have daily observations). \\ \\ \noindent The delta method can, once again, be used to obtain estimated standard errors for such return levels. We should, however, also include uncertainty in our estimate of $\lambda_{u}$ in the calculation (since $z_{r}$ is a function of $\lambda_{u}$). Since the number of threshold exceedances follows a binomial distribution $\textrm{Bin}(N, \lambda_{u})$, where $N$ is the total number of observations in the series, we know (from MAS2302) that \[ Var(\hat{\lambda}_{u}) = \hat{\lambda}_{u}(1-\hat{\lambda}_{u})/N = 0.0007^{2} \] in our rainfall example. Then, by the delta method, \begin{eqnarray*} \textrm{Var}(\hat{z}_{r}) &\approx& \nabla z_{r}^{T}V\nabla z_{r}. \end{eqnarray*} \noindent Here, $V$ is now the variance--covariance matrix of the triple $(\hat{\lambda}_{u},\hat{\sigma},\hat{\xi})^{T}$; 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$. The value --0.0656 is the estimated covariance between $\hat{\sigma}$ and $\hat{\xi}$, and is found in \texttt{A\$cov} (recall that we stored the fit of the GPD to our rainfall extremes in \texttt{A}): \begin{verbatim} > A$cov [,1] [,2] [1,] 0.91849139 -0.06546151 [2,] -0.06546151 0.01023053 \end{verbatim} We can now proceed as in Section 2.2.4 to obtain standard errors. \vspace{0.5cm} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} Obtain the standard error for the 50--year return level estimate of rainfall at New York. \huge \ding{46} \normalsize \vspace{18cm} }} \newpage \subsection{Profile likelihood} Recall that confidence intervals for return levels constructed using standard errors are usually not appropriate; rather, we should construct intervals using \textit{profile likelihood}. We have already considered the idea of using the profile likelihood to obtain more realistic confidence intervals for return levels in Section 2.4; we now implement the \textsf{ismev} command \texttt{gpd.prof} to construct the profile log--likelihood for our estimate of the 50 year return level at New York. \\ \\ Recall that we previously stored our results of fitting the GPD to the set of rainfall extremes observed at New York in \texttt{A}, that is: \begin{verbatim} > A=gpd.fit(rain,30) \end{verbatim} Now typing: \begin{verbatim} > gpd.prof(A,50,xlow=70,xup=150,npy=365.25) \end{verbatim} produces the plot shown in Figure 3.2. Notice that you need to provide the \texttt{gpd.prof} function with a range of values at which to fix the return level, before maximising the log--likelihood function with respect to the remaining parameters. Here, we have chosen a range from 70$\rightarrow$150, although this sometimes requires experimentation. Notice also that you need to tell the function the (average) number of observations per year (365.25 here, to account for leap years). \\ \\ Just to check, you can also superimpose your estimated return level on the plot, just to make sure all is well. Recall from Section 3.2.5 that this was \[ \hat{z}_{50}=92.24 \textrm{ mm}, \] and so the code \begin{verbatim} > abline(v=92.24) \end{verbatim} inserts the vertical line that can be seen on Figure 3.2. As you can see, this is at the mode of the profile log--likelihood -- as it should be! The default confidence interval is the 95\% confidence interval, and so the \texttt{gpd.prof} function automatically places a horizontal line at a distance of $\frac{1}{2}\chi^{2}_{1}(0.05)=1.921$ from the maximised value of the profile log--likelihood (see Section 2.4, and the supplementary handout on profile likelihood, for further details). This gives us a 95\% confidence interval, based on profiling the log--likelihood, of about (74.1mm, 143mm). Planners and civil engineers often design to the upper end--point of such confidence intervals, just to make sure they don't ``under--protect''. When feeding back such information to non--statisticians, we often say something like ``once every fifty years, we might expect daily rainfall accumulations in New York City to reach up to about 143mm''. \begin{figure}[!h] \centering \psfrag{-485}[][]{\sf \tiny --485} \psfrag{-486}[][]{\sf \tiny --486} \psfrag{-487}[][]{\sf \tiny --487} \psfrag{-488}[][]{\sf \tiny --488} \psfrag{-489}[][]{\sf \tiny --489} \psfrag{80}[][]{\sf \tiny 80} \psfrag{100}[][]{\sf \tiny 100} \psfrag{120}[][]{\sf \tiny 120} \psfrag{140}[][]{\sf \tiny 140} \psfrag{Profile Log-lik}[][]{\sf \tiny Profile log--likelihood} \psfrag{Return}[][]{\sf \tiny Return level} \psfrag{Z100}[][]{\sf \scriptsize $z_{100}$} \psfrag{Z1000}[][]{\sf \scriptsize $z_{1000}$} \includegraphics[width=16cm,height=12cm]{profileny.ps} \caption{Profile log--likelihood curves for the 50 year return level daily rainfall accumulation at New York.} \end{figure} \vspace{0.5cm} \framebox[17cm]{ \parbox[t]{16cm}{ \vspace{0.4cm} Compare the 95\% confidence interval for the 50--year return level obtained from profiling the log--likelihood to that you would obtain by using the standard error. Comment. \huge \ding{46} \normalsize \vspace{8cm} }}