
\section{The Peaks-Over-Threshold Method}\label{sec:pot}\index{peaks-over-threshold method}\index{extreme value distribution!peaks-over-threshold method} 

For many physical processes, the assumption of temporally independent observations is unrealistic. Stationarity is usually a more plausible assumption, and it says that even though observations or data may be dependent, their stochastic properties are temporally homogeneous. A GEV distribution remains an appropriate model for block maxima of stationary series, and a GP distribution can also be shown to remain appropriate for threshold excesses, see \cite{begosete04}. 

The degree of dependence between data points is in some sense quantified in the \textit{extremal index} $\theta$, which satisfies $0<\theta\leq 1$, cf. Section~\ref{sec:dss}. This parameter can be interpreted as a measure of the tendency of the process to cluster at extreme levels, and one can informally say that the inverse of the extremal index is the limiting mean cluster size. This means that if $\theta=0.5$, then extreme values would approximately arrive in groups of two. 

For the block maxima case, dependence is mostly absorbed in the parameters, which have to be estimated anyway. However, some change is necessary for threshold excess modeling. A common method used to overcome this issue is \textit{declustering}, where the generalized Pareto distribution is instead fitted to the maxima of clusters. These clusters are identified by some empirical rule. 

The easiest method of cluster identification is taking a threshold and saying that subsequent observations must be above this threshold to be part of the current cluster. A modification can be made by allowing one or more subsequent observations to be below the threshold before the current cluster is left. Selecting the amount of subsequent observations allowed before a cluster is left, is done by selecting a number for $r$, the allowed distance between observations above the threshold. For instance, this means that if observation number 50 and observation number 55 are both above the threshold $u$, they would be in the same cluster for $r\geq5$ since the distance between the observations is 5. For $r<5$, they would be considered to be from different clusters.

To be able to use a POT model in practice, the threshold $u$ must be selected, as discussed in Chapter~\ref{chap:pot}. In the mean residual life plots, cf. Eq.~\eqref{eq:pot.15}, for the three data sets that are analyzed in this chapter, the choice of $u_0$ has been clearly marked. Another method is available, based on fitting the model to a wide range of thresholds. In Eq.~\eqref{eq:pot.14} it was used that there is a linear relationship between scale parameters $\sigma_u$ above a valid threshold $u_0$. Furthermore, shape parameters $\gamma_u$ should be constant. By reparametrizing the scale parameter to $\sigma_u^* = \sigma_u - \gamma_u u$, one obtains a parameter which should also be constant above $u_0$. The plot of $(u,\gamma_u)$ and $(u,\sigma_u^*)$ is hereby called the \textit{stability plot}. Stability plots are also presented for the three available data sets.

Having chosen the reference level $u_0$, one would typically calculate the return period levels for several values of $u \geq u_0$. 
This is obtained by using the following formulas:
\begin{equation}
x_N = u + \sigma[(N n_y \zeta_u \theta)^{\gamma}-1]/\gamma,
\label{eq:returnxthetaN}
\end{equation}
for the case of $\gamma\neq 0$, and
\begin{equation}
x_N = u + \sigma \ln(N n_y \zeta_u \theta),
\label{eq:returnx0thetaN}
\end{equation}
for $\gamma=0$. Here $N$ = number of years, $n_y$ = the average number of observations per year, $\zeta_u = \Prob(X>u)$, $\theta$ = the extremal index.



The parameters of the Generalized Pareto distributions in this section are estimated by maximum likelihood. This was done with the maximum likelihood procedure in the \texttt{fpot} function in the \texttt{evd} library in R. Standard errors for the GP parameters come from this R function as well, and are extracted from a numerical approximation of the observed information. Confidence intervals for the parameters are then calculated by using the approximate normality of the maximum likelihood estimator. The uncertainty in $\zeta_u$ is ignored since it is usually small compared to the errors of the other parameters \citep{cole01}. 



\subsection{ Application to water level measurements}
%%

\subsubsection*{Oslo}\index{water level measurements!Oslo}

Figure~\ref{pot_meanexceed_oslo} shows the mean residual life plot for Oslo. To select an appropriate threshold, one would search for approximate linearity within the confidence bounds. The vertical line shows the choice made at 134 cm. The figure does not paint a completely clear picture, as is the case with most real data sets, but from around the indicated spot and up to about 200 cm there is a certain level of linearity. 

Figure~\ref{pot_stability_oslo} shows the effect of threshold selection on the model parameters. Ideally, the shape and modified scale parameters should be invariant to threshold change as long as the data are above the minimum threshold. One can see that this is basically the case from the beginning of the plot and up to about 170-175~cm. A threshold of 134~cm is within this range, and cannot be rejected based on Figure~\ref{pot_threshold_oslo}. 
\begin{figure}[!htb]
\centering
\subfigure[Mean of exceedances of increasing thresholds (cm).]{
\includegraphics[scale=0.35]{EWL/Figures/pot_meanexceed_oslocorrected.pdf}%
\label{pot_meanexceed_oslo}%
}
\subfigure[Stability of parameters for increasing threshold (cm).]{
\includegraphics[scale=0.35]{EWL/Figures/pot_stability_oslocorrected1.pdf}
\label{pot_stability_oslo}
}
\caption{Threshold selection plots for the POT model for Oslo.}
\label{pot_threshold_oslo}
\end{figure}

Having selected a threshold, one would go on to select the allowed distance between points in a cluster, $r$. Figure~\ref{pot_exiplot_oslo} shows the behaviour of the extremal index as the allowed distance between observations in a cluster is increased. The extremal index declines quite sharply until about $r=12$, then almost flattens out before declining faster again. A similar tapering in the steepness of the curve is seen at about $r=24$, and possibly at about $r=36$. 

With hourly observations, these changes in steepness coincide with the 12 hour cycle of the lunar tidal component. It is therefore safe to assume that extreme storms above $u=134$ may last from one high tide to another. This suggests selecting a value of $r$ larger than 12 to encapsulate the full length of such storms. But it does not seem plausible to have $r$ much larger, since if the tide was the dampening factor, then the subsequent rising tide should bring the storm with it up to extreme levels again. 

More insight into the choice of $r$ is granted by looking at Figure~\ref{pot_returnchange_oslo}, which shows the development of the 200 year return level as $r$ is increased from 1 to 40. The return level stays approximately constant from $r=1$ to $r=5$, but increases relatively sharply from there and up to $r=8$. It then decreases just as sharply, before staying relatively constant from $r=12$ up to about $r=22$. It then drops sharply, before staying constant up to $r=36$. For $r\in (1,5)$ and $r\in(12,22)$, the return level is approximately the same. Having argued that $r$ should not be much larger than 12, a value of $r=1$ is chosen since this gives approximately the same result while using more data. Note that $r=1$ means that a new cluster starts as soon as one value is below the threshold.

\begin{figure}[!ht]
\centering
\subfigure[Extremal index for $r$ between 1 and 40.]{
\includegraphics[scale=0.35]{EWL/Figures/pot_exiplot_oslocorrected.pdf}%
\label{pot_exiplot_oslo}%
}
\subfigure[200 year return level for $r$ between 1 and 40.]{
\includegraphics[scale=0.35]{EWL/Figures/pot_returnchange_oslocorrected.pdf}%
\label{pot_returnchange_oslo}%
}
\caption{Change in extremal index and 200 year return level as calculated by the POT model for Oslo.}
\label{pot_rchoice_oslo}
\end{figure}
With the threshold at 134 cm, the model parameters were estimated to be $\hat{\gamma}=0.019 (0.025)$ and $\hat{\sigma}=13.993 (0.494)$ by the \texttt{fpot} function in the \texttt{evd} package in $R$. Probability and quantile plots indicate that the model with estimated parameters fits well to the maxima of the identified clusters. From the 778199 points in the data series, 7364 points are above the chosen threshold, giving an estimate of the exceedance probability $\zeta$ of $\hat{\zeta}=0.0095$. From the 7364 points above the threshold, 1675 clusters were identified, giving an extremal index of $\hat{\theta}=0.2275$. All this gives the return level estimates shown in Table~\ref{pot_table}. 





The return level plot is shown in Figure~\ref{pot_rlevel_oslo}. All data points are contained within the confidence bounds of the fitted model, and most of them stay on or very close to the line. Points above 185 cm behave somewhat more erratically than those below, but not to a dramatic extent. This level corresponds approximately to the level at which the stability of parameters, as displayed in Figure~\ref{pot_stability_oslo}, begins to wear off. 


\begin{figure}[!htb]
\centering
\includegraphics[scale=0.35]{EWL/Figures/pot_rlevel_oslocorrected.pdf}
\caption{Return levels (cm) for increasing return periods (years) for Oslo.}
\label{pot_rlevel_oslo}
\end{figure}




\subsubsection*{Heimsjø}\index{water level measurements!Heimsjø}

Figure~\ref{pot_meanexceed_heimsjo} shows the mean residual life plot for the Heimsjø data. Approximate linearity is found around the vertical line, which is placed at $u=287$ cm. The stability plot in Figure~\ref{pot_stability_heimsjo} shows the effect on the model parameters of threshold change around the 287 cm mark. Above approximately 280 cm and up to around 325-330 cm, the parameters are approximately constant. 

\begin{figure}[!htb]
\centering
\subfigure[Mean of exceedances of increasing thresholds (cm).]{
\includegraphics[scale=0.35]{EWL/Figures/pot_meanexceed_heimsjo.pdf}%
%\caption{Mean of exceedances of increasing thresholds (cm)., Oslo.}%
\label{pot_meanexceed_heimsjo}%
}
\subfigure[Stability of parameters for increasing threshold (cm).]{
\centering
\includegraphics[scale=0.35]{EWL/Figures/pot_stability_heimsjo1.pdf}
%\caption{Stability of parameters for increasing threshold (cm).s for Oslo.}
\label{pot_stability_heimsjo}
}
\caption{Threshold selection plots for the POT model for Heimsjø.}
\label{pot_threshold_heimsjo}
\end{figure}

 An appropriate value of $r$ is sought for the identification of independent clusters and the estimation of the extremal index. Figure~\ref{pot_exiplot_heimsjo} shows how the extremal index estimate behaves. It is close to piecewise constant, with significant drops at around $r=10$ to 12 and at $r=24$. Figure~\ref{pot_returnchange_heimsjo} shows that there is quite little difference between the estimate for $r=1$ and larger $r$. The difference is somewhat larger than was found for Oslo, however, and $r$ is chosen to include storms crossing from one tide to another. This leads to the choice $r=13$.

\begin{figure}[!htb]
\centering
\subfigure[Extremal index for $r$ between 1 and 40.]{
\includegraphics[scale=0.35]{EWL/Figures/pot_exiplot_heimsjo.pdf}%
\label{pot_exiplot_heimsjo}%
}
\subfigure[200 year return level for $r$ between 1 and 40.]{
\includegraphics[scale=0.35]{EWL/Figures/pot_returnchange_heimsjo.pdf}%
\label{pot_returnchange_heimsjo}%
}
\caption{Change in extremal index and 200 year return level as calculated by the POT model for Heimsjø.}
\label{pot_rchoice_heimsjo}
\end{figure}

With the 287 cm threshold, model parameters are estimated as $\hat{\sigma}=13.709 (0.653)$ and $\hat{\gamma}=-0.116 (0.030)$. Probability and quantile plots show a fair agreement with the fitted model. However, for the highest sea levels, there is an issue with the model fit, which will be discussed below. From the 686867 points in the Heimsjø data set, 2363 points are above 287 cm, giving an estimate of the exceedance probability $\zeta$ of $\hat{\zeta}=0.0034$. From the points above the threshold, 733 clusters were identified, giving an extremal index of $\hat{\theta}=0.310$.

Return level estimates are presented in Table~\ref{pot_table}. They go from slightly overshooting to slightly undershooting the published return levels, but all are certainly within confidence bounds. 



As mentioned above, for the highest sea levels, there are some issues with the model fit. This is better seen by looking at the return level plot in Figure~\ref{pot_rlevel_heimsjo}, which shows irregularities after about 330 cm - corresponding to the level after which parameters in Figure~\ref{pot_stability_heimsjo} are no longer near-constant. The five highest points are above the model line, and it seems possible that return levels for long return periods can be underestimated. 


\begin{figure}[!htb]
\centering
\includegraphics[scale=0.35]{EWL/Figures/pot_rlevel_heimsjo.pdf}
\caption{Return levels (cm) for increasing return periods (years) for Heimsjø.}
\label{pot_rlevel_heimsjo}
\end{figure}



\subsubsection*{Honningsvåg}\index{water level measurements!Honningsvåg}

The mean residual plot in Figure~\ref{pot_meanexceed_honning} shows tendencies of linearity within confidence bounds slightly before the 280 cm mark. A threshold of 278 cm is selected, indicated by the vertical black line. The stability plot in Figure~\ref{pot_stability_honning} is far from constant up to around 275-280 cm, but from there it is reasonable to call the parameters near-constant up to perhaps 310 cm. The threshold $u=278$ is therefore barely within the acceptable area as far as the stability plot is concerned. 

\begin{figure}[!htb]
\centering
\subfigure[Mean of exceedances of increasing thresholds (cm).]{
\includegraphics[scale=0.35]{EWL/Figures/pot_meanexceed_honning.pdf}%
%\caption{Mean of exceedances of increasing thresholds (cm)., Oslo.}%
\label{pot_meanexceed_honning}%
}
\subfigure[Stability of parameters for increasing threshold (cm).]{
\centering
\includegraphics[scale=0.35]{EWL/Figures/pot_stability_honning1.pdf}
%\caption{Stability of parameters for increasing threshold (cm).s for Oslo.}
\label{pot_stability_honning}
}
\caption{Threshold selection plots for the POT model for Honningsvåg.}
\label{pot_threshold_honning}
\end{figure}

To find an appropriate level of $r$, Figure~\ref{pot_exiplot_honning} is inspected. Honningsvåg is tide-dominant like Heimsjø, and the same pattern is shown here as was observed there; the extremal index plot is approximately piecewise constant. Figure~\ref{pot_returnchange_honning} shows a similar pattern as that for Heimsjø as well, but the differences in return levels are much larger. The values for $r\leq 10$ and $r>24$ are about the same, however. The difference is large enough to not be within the confidence intervals, and $r=13$ is chosen since it was argued earlier for a $r\geq 12$. 

\begin{figure}[!htb]
\centering
\subfigure[Extremal index for $r$ between 1 and 40.]{
\includegraphics[scale=0.35]{EWL/Figures/pot_exiplot_honning.pdf}%
%\caption{Extremal index for $r$ between 1 and 40., Oslo.}%
\label{pot_exiplot_honning}%
}
\subfigure[200 year return level for $r$ between 1 and 40.]{
\includegraphics[scale=0.35]{EWL/Figures/pot_returnchange_honning.pdf}%
%\caption{Extremal index for $r$ between 1 and 40., Oslo.}%
\label{pot_returnchange_honning}%
}
\caption{Change in extremal index and 200 year return level as calculated by the POT model, for Honningsvåg.}
\label{pot_rchoice_honning}
\end{figure}

A threshold of 278 cm gives estimated model parameters of $\hat{\sigma}=14.692 (0.530)$ and $\hat{\gamma}=-0.025 (0.026)$. Probability and quantile plots show a fair agreement with the fitted model. However, as for Heimsjø, for the highest sea levels, there is an issue with the model fit. There are 330883 sea level measurements in the Honningsvåg data, 9264 of which are above 278 cm. This makes for an estimate of the exceedance probability $\zeta$ of $\hat{\zeta}=0.0280$. From the points above the threshold, 1590 clusters were identified, giving an extremal index of $\hat{\theta}=0.172$.

Estimated return levels are presented in Table~\ref{pot_table}.



As indicated above, the quantile plot shows that a few of the most extreme data points are off the fitted model curve, with Figure~\ref{pot_rlevel_honning} giving a clearer impression. It looks likely that higher return levels may be overestimated, although all points are within confidence intervals.


\begin{figure}[!htb]
\centering
\includegraphics[scale=0.35]{EWL/Figures/pot_rlevel_honning.pdf}
\caption{Return levels (cm) for increasing return periods (years) for Honningsvåg.}
\label{pot_rlevel_honning}
\end{figure}



\begin{table}[!htb]
\begin{center}
	\caption{Return level estimates and 95\% confidence intervals for the POT model. }
	\label{pot_table}
	\vspace{1mm}
\begin{tabular}{c|c|c|c|c|c|c} 
\multicolumn{1}{c|}{}&\multicolumn{2}{c|}{Oslo}&\multicolumn{2}{c|}{Heimsjø} &\multicolumn{2}{c}{Honningsvåg}\\ 
 \hline 
\multicolumn{1}{c|}{$R = 1/p$}&\multicolumn{1}{c|}{$\hat{z}_p$}&\multicolumn{1}{c|}{CI}&\multicolumn{1}{c|}{$\hat{z}_p$}&\multicolumn{1}{c|}{CI}&\multicolumn{1}{c|}{$\hat{z}_p$}&\multicolumn{1}{c}{CI} \\ \hline
5 & 200.5 & (189.3, 213.9) & 329.5 & (321.7, 338.8) & 351.6 & (338.1, 368.3) \\ 
10& 211.1 & (197.0, 228.5) & 335.3 & (325.9, 347.0) & 360.5 & (344.2, 380.9) \\ 
20& 221.9 & (204.7, 243.8) & 340.7 & (329.5, 354.8) & 369.2 & (350.1, 393.8)) \\ 
100& 247.6 & (221.7, 282.2) & 351.7 & (336.5, 372.0) & 388.8 & (362.5, 424.6) \\ 
200& 258.9 & (228.8, 300.0) & 355.8 & (338.9, 378.9) & 397.0 & (367.4, 438.3) \\
\hline
\end{tabular}

\end{center}
\end{table}
