
\section{Revised Joint Probabilities Method}\label{sec:rjpm}\index{extreme value distribution!revised joint probabilities method}\index{revised joint probabilities method}\index{RJP method} 


The revised joint probabilities (RJP) method is an attempt by \cite{tawn92} to improve the joint probabilities method employed by \cite{Pugh80}. Whereas they assume that hourly surge levels are independent, Tawn argues that this is clearly a false assumption. Instead, he uses that for a stationary sequence $Y_1,...,Y_n, \ldots$,
\begin{equation}
\Prob(\max \{Y_1,...,Y_n\} < y)\approx [\Prob(Y_1<y)]^{n\theta}
\label{eq:mixing}
\end{equation}
for large $y$, and $0 < \theta \leq 1$. $\theta$ is the extremal index, and $\theta^{-1}$ is defined as the limit of $\theta^{-1}(y)$ as $y$ tends towards the upper end point of the distribution of $Y$. $\theta^{-1}(y)$ is defined as the mean of the distribution of cluster sizes. The extremal index can be equal to 1 for both dependent and independent sequences, when the sequence behaves like an independent sequence at high levels. Unfortunately, as will be demonstrated in the section on the ACER method, for the data analysed in this chapter, $\theta(y)$ does not seem to be a robust parameter relative to its dependence on $y$. In fact, it will be seen that $\theta(y)$ may display significant dependence on $y$ while there is clear evidence that $\theta(y) \to 1$ when $y$ increases. However, in the present chapter, the extremal index is estimated as would be typically done in the RJP method.


The RJP method relies on a componentwise analysis of the sea level $Z$. It is divided into three components, $Z_t = M_t+X_t+Y_t,$
where $M_t$ is the mean sea level, $X_t$ the tidal level and $Y_t$ the surge level. For the purpose of this discussion, $M_t=0$ for all $t$. The tidal level $X_t$ is estimated by a numerical model. 


If $N$ is the number of hours in a year, the surge levels $Y_1,...Y_N$ are taken to be a realization of a stationary sequence. It is then assumed that \citep{cole01} $\Prob(\mbox{max}\{Y_1,...,Y_N\}\leq y) = F_s(y)^{N\theta_s} = G(y; \mu_s,\sigma_s,\gamma_s)$
for large values of $y$, say $y>u$. $F_s$ is the marginal distribution for surges, while $\theta_s$ ($0\leq\theta_s\leq 1$) is the extremal index for surges. $G$ refers to the GEV distribution model. It follows that, 
\begin{equation}
F_s(y)=\exp\left\{-(N\theta_s)^{-1}\left[1+\gamma_s\frac{y-\mu_s}{\sigma_s}\right]^{-1/\gamma_s}\right\} \ \mbox{ for $y>u$},
\label{eq:surgedist}
\end{equation}
on $\{y:1+\gamma_s(y-\mu_s)/\sigma_s>0\}$, with $\sigma_s>0$ and arbitrary $\gamma_s$ and $\mu_s$. 


If $T$ is the tidal cycle length in hours, it follows that \citep{tawn92},
\begin{equation}
\Prob(\mbox{max}\{Z_1,...,Z_T\}\leq z) = \Prob(\bigcap_{t=1}^T Y_t\leq z-X_t\} =  \left\{\displaystyle\prod_{t=1}^T{F_s(z-X_t)}\right\}^{\theta} ,
\label{eq:tidalcycle}
\end{equation}
where $\theta$ is an hourly sea level extremal index ($\theta_s\leq \theta \leq 1$). Combining Eqs.~\eqref{eq:surgedist} and \eqref{eq:tidalcycle}  gives the RJP distribution for annual maximum sea levels,
\begin{equation}
G(z)=\exp\left\{-\theta(T\theta_s)^{-1}\displaystyle\sum_{t=1}^T\left[1+\gamma_s\frac{z-\mu_s-X_t}{\sigma_s}\right]^{-1/\gamma_s}\right\},
\label{eq:rjpmdist}
\end{equation}
where $z>u+\max\{X_1,..., X_T\}$. 


Eq.~\eqref{eq:rjpmdist} shows the case where the surge distribution $F_s(\cdot)$, and therefore its parameters $\mu_s,\sigma_s$ and $\gamma_s$, are independent of the concomitant tidal level. Unfortunately, the assumption that tide and surge are independent processes is poor in shallow water areas, where turbulent frictional processes on the sea bed cause the tide and surge components to interact. This causes effects such as surge values at high tides being damped and surges on the rising tide being amplified \citep{Dixon1994}. These effects vary from site to site however, so it is attempted to model the interaction on the residuals from the observed surges. 


The tidal range from lowest observed tide (LAT) to highest observed tide (HAT) is split into $n_b$ equi-probable bands, i.e. each band has an equal amount of observed measurements. For each of the tidal observations, there is a concurrent surge observation. This means that if there are $n_b$ tidal bands and $n_{obs}$ observations in total, there are $n_{obs}/n_b$ observations of both the tide and surge in each band. 

If the tide and surge were independent, an equal amount of points should also be expected to exceed a given level $u$. But if there is interaction, then the least number of points should be in the top band where surges are damped. Similarly, the largest number should be in the middle bands where surges are magnified. The amount of discrepancy between bands will therefore be a quantifiable measure of the tide-surge interaction. 

The level $u$ is now chosen to be a high empirical quantile of the surge distribution, $z_q$ for $q=0.9975$. For the independent case, there should then be $(1-q)\cdot n_{obs}/n_b=v$ observations in each band. $n_b=5$ is chosen as in \cite{Dixon1994}, meaning $n_{obs}\cdot 0.0025/5=v$ observations per tidal band. In actuality, since independence is a flawed assumption, $N_i$ surges per band are observed, for $i=1,...,n_b$. To put the interaction into a quantifiable setting, a standard $\chi^2$ test statistic is used, $\chi^2=\sum_{i=1}^5 (N_i-v)^2 / v^2$. If $N_i\approx v$, then $\chi^2$ will be small. Tide and surge are deemed to interact with 95 \% confidence if the test statistic is above the associated 4 degrees of freedom table value of $\chi_{4,0.95}^2=9.488$. 


To account for the tide-surge interaction, the method used in \cite{Dixon1994} is adopted, where the surge series is location-scale normalized by
$S_t^*= (Y_t-a(X_t))/b(X_t)$, where $a(X_t)$ and $b(X_t)$ are some tide-dependent functions. $\{S_t^*\}$ is then supposed to be stationary, and established methods can be used to estimate the associated model parameters $\mu_{s^*}$, $\sigma_{s^*}$ and $\gamma_{s^*}$. The parameter estimates for the original surge series are then given by $\mu_s(X) = \mu_{s^*}b(X)+a(X)$, $\sigma_s(X) = \sigma_{s^*}b(X)$, and $\gamma_s(X) = \gamma_{s^*}$.

Equation \eqref{eq:rjpmdist} is then modified to 
\begin{equation}
G(z) = \exp\left\{-\theta(T\theta_s)^{-1}\displaystyle\sum_{t=1}^T\left[1+\gamma_{s^*}\frac{z-\mu_{s^*}b(X_t)-a(X_t)-X_t}{(\sigma_{s^*}b(X_t))}\right]^{-1/\gamma_{s^*}}\right\},
\label{eq:rjpmdistrev}
\end{equation}
for $z > u+\max\{X_1,..., X_T\}$. The estimation of $a(X)$ and $b(X)$ is discussed in \cite{Dixon1994}.



\subsection{ Estimating return levels with the RJP method}

In the 1992 and 1994 reports \citep{tawn92,Dixon1994}, the $r$-largest method is applied to estimate the parameters of the surge distribution. This method is a modification of the GEV annual maxima method, and uses the $r$ largest values per year to estimate model parameters. 

The threshold applied in the estimation of $\theta_s$ is now adopted, and the POT method is then used to estimate model parameters. This is a method which has the advantage of using much more data. The fact that for a POT model $\gamma_{POT}=\gamma_{GEV}=\gamma$ for the corresponding GEV model is used, with $\sigma_{POT}=\sigma_{GEV}+\gamma\cdot(u-\mu_{GEV})$ \citep{cole01}.  

In practice, the shape and scale are estimated with a maximum likelihood procedure by the function \texttt{fpot} in the $R$ package \texttt{evd}, while the location parameter is calculated by the \texttt{fgev} function. This is also the case for the confidence intervals of model parameters. The model is fitted to the maxima of clusters, where the clusters are identified by the same rule that governs the estimation of $\theta$. 


The extremal indices $\theta$ and $\theta_s$ require two choices each; a threshold above which clusters are counted and an empirical clustering rule which says how many non-exceedances are allowed before the current cluster is terminated. For the first choice, the quantile that was found in the POT analysis is used. There a thorough analysis was done on where the threshold should be placed, and it makes sense using this same threshold to estimate $\theta$. Applying the procedure of using the same quantiles for both extremal indices, as was done in \citep{Dixon1994}, corresponding thresholds are obtained for the surge series and the estimation of $\theta_s$ as well. 

For the second choice, it is necessary to decide what constitutes an independent storm. In \cite{Dixon1994}, they found that $r=30$ was a good choice, and that the ratio of extremal indices was not too dependent on this choice in any case. In the present analysis the results are more sensitive to this choice, since the POT method with clustering is used to estimate model parameters. Still this $r$ value of 30 is adopted.



A practical challenge with the RJP method is that Eq.~\eqref{eq:rjpmdistrev} is impossible to solve analytically for $z$. One may start out by defining $G(z_p)=1-p$ and get the relation,
\begin{equation}
-T\theta_s\theta^{-1}\ln(1-p)=\displaystyle\sum_{t=1}^T\left[1+\gamma_{s^*}\frac{z_p-\mu_{s^*}b(X_t)-a(X_t)-X_t}{(\sigma_{s^*}b(X_t))}\right]^{-1/\gamma_{s^*}},
\label{eq:rjpmreturn}
\end{equation}
but one cannot go further by algebraic manipulation. A numerical procedure is therefore developed where the right hand side of the equation is calculated for a range of relevant $z$ values and then matched to the left side of the equation. 

The tidal cycle length is 18.61 years, meaning that $T=18.61\cdot 8766 =163135.3\approx 163135$. Any span of consecutive 163135 observations should therefore approximately contain a full tidal cycle, and in the numerical procedure the last 163135 observations in the tidal series are taken. The numerical analysis code is written in R.



\subsection{ Application to water level measurements}



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

A scrutiny of the sea level data for Oslo reveals that a large percentage of the observed sea level rises stem from surges. Locations with such characteristics are called surge-dominant, as opposed to tide-dominant where most of an observed sea level stems from the current tidal level. Areas where surge dominates usually do not have as much tide-surge interaction as tide-dominant areas, but the $\chi^2$ test is still used to quantify the amount. 

By using 5 tidal bands and cutting the data at the 99.75\% quantile, 389 surges per band are expected. Instead 346, 380, 398, 431 and 386 data points per band are found, in order from the lowest tide level to the highest. This is not a very bad result, but there is still evidence of significant tide-surge interaction. The test statistic has a result of $\chi^2 = 9.778$, compared to the table value of $\chi_{4,0.95}^2 = 9.488$. 

Since the $\chi^2$ test shows a significant level of interaction, it is modelled by using 5 tidal bands. Figure~\ref{rjpm_tidesurgeint_oslo} shows the tide against surge data to the left, and tide against transformed surge data to the right, both showing only points above the 99.75\% quantile. There seems to be very little difference between the left and the right plots except for the scales on the $x$ axes. 

Figure~\ref{rjpm_chiratio_oslo} shows how the ratio between the $\chi^2$ test statistic and the corresponding $\chi_{n_b-1,0.95}^2$ value develops after 2-30 bands have been used to transform the surge data. $\chi^2/\chi_{n_b-1,0.95}^2=1$ is indicated by the horizontal line, and ratios above this correspond to $\chi^2$ tests showing significant interaction. It is seen that for $n_b=2$ and $n_b=15$ the interaction is insignificant, but for $n_b=2$ the interaction is insignificant for the untransformed data as well, with a test statistic of 1.85 versus the corresponding table value of 3.84. Transforming the data with $n_b=15$ seems needlessly complex for a model where there was hardly significant interaction in the original test. Furthermore, there is little actual difference in return levels. The same is true if they are compared to the return levels achieved using no tide-surge correction. For parsimony the model without correction is selected, based on Eq.~\eqref{eq:rjpmdist} instead of Eq.~\eqref{eq:rjpmdistrev}. Estimated parameters are $\hat{\mu}_s=95.655 (2.313)$, $\hat{\sigma}_s=17.576 (1.053)$, $\hat{\gamma}_s=-0.049 (0.039)$, $\hat{\theta}=0.079$ and $\hat{\theta}_s=0.068$. 

Resulting return levels and confidence intervals for Oslo can be found in Table~\ref{rjpm_table}.  



\begin{figure}[!htb]
\centering
\subfigure[Tide levels and concurrent surges above the 99.75\% quantile.]{
\includegraphics[scale=0.35]{EWL/Figures/rjpm_tidesurgeint_oslocorrected.pdf}%
%\caption{Tide levels and concurrent surges above the 99.75\% quantile for Oslo.}%
\label{rjpm_tidesurgeint_oslo}
}
\subfigure[Ratio between $\chi^2$ test statistic and table value $\chi_{n_b-1,0.95}^2$.]{
\includegraphics[scale=0.35]{EWL/Figures/rjpm_chiratio_oslocorrected.pdf}%
%\caption{The interaction functions $\hat{a}_1(X)$, $\hat{a}_2(X)$ (red) for 5 and 10 tidal bands in Oslo.}%
\label{rjpm_chiratio_oslo}
}
\caption{Tide-surge interaction plots for Oslo.}
\label{rjpm_tidesurge_oslo}
\end{figure}




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

Heimsj� is tide-dominant, and it is observed from the data that there is quite a small distance between the maximum of its tidal series and the maximum of observed sea levels. This also holds for the mean of annual maxima from tide calculations and observations. Tide-dominance usually implies larger tide-surge interaction than surge-dominance, and the $\chi^2$ test for Heimsj� certainly lives up to this expectation. With 5 tidal bands and cutting at the 99.75\% quantile, about 343 surge observations are expected per band. Instead, 790, 479, 288, 112 and 75 points are found. As expected for a tide-dominant location, the least number of points are in the highest tidal band, where the high tide dampens surges. By far the most points are in the bottom band. 

The huge discrepancy between expected and observed number of surges in each band leads to a massively large value for the $\chi^2$ test statistic, with $\chi^2=1007.3$ compared to the table value of $\chi_{4,0.95}^2=9.488$. The results of transforming the data with 2 to 30 bands are shown in Figure~\ref{rjpm_chiratio_heimsjo}. The ratio never drops below 1 as one would wish, but instead $n_b=12$ is chosen, which gives a relatively low ratio compared to other choices. It gives $\chi^2=89.865$ versus the quantile value of $\chi_{11,0.95}=19.675$. Figure~\ref{rjpm_tidesurgeint_heimsjo} shows that the number of points in each band is much more similar in the transformed case to the right of the figure compared to the left, uncorrected side. 

Estimated parameters are $\hat{\mu}_{s^*}=3.494 (0.283)$, $\hat{\sigma}_{s^*}=1.671 (0.160)$, $\hat{\gamma}_{s^*}=0.074 (0.079)$, $\hat{\theta}=0.106$ and $\hat{\theta}_{s^*}=0.063$. 

Return level and confidence intervals for Heimsj� are presented in Table~\ref{rjpm_table}. The low return period levels are quite different from their published counterparts, but the upper confidence bound of the 20 year return level estimate contains the corresponding published value \citep{Hansen2000}. This could perhaps be explained by the tide-surge interaction still present, but estimation of return levels without any tidal bands also show this apparent underestimation of low return period levels. The standard error of the estimated $\gamma_s$ is larger than the estimate itself, and the estimate being positive means very large upper bounds for long return periods. 
%\vspace*{20mm}

%
\begin{figure}[!htb]
\centering
\subfigure[Tide levels and concurrent surges above the 99.75\% quantile.]{
\includegraphics[scale=0.35]{EWL/Figures/rjpm_tidesurgeint_heimsjo.pdf}%
%\caption{Tide levels and concurrent surges above the 99.75\% quantile for Oslo.}%
\label{rjpm_tidesurgeint_heimsjo}
}
\subfigure[Ratio between $\chi^2$ test statistic and table value $\chi_{n_b-1,0.95}^2$.]{
\includegraphics[scale=0.35]{EWL/Figures/rjpm_chiratio_heimsjo.pdf}%
%\caption{The interaction functions $\hat{a}_1(X)$, $\hat{a}_2(X)$ (red) for 5 and 10 tidal bands in Oslo.}%
\label{rjpm_chiratio_heimsjo}
}
\caption{Tide-surge interaction plots for Heimsj�.}
\label{rjpm_tidesurge_heimsjo}
\end{figure}




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

Honningsv�g is tide-dominant like Heimsj�, but the $\chi^2$ test shows much less interaction than for Heimsj�. It is still significant on a 95\% confidence level, with $\chi^2=23.652$. Figure~\ref{rjpm_chiratio_honning} shows how the $\chi^2$ ratio develops. After $n_b=2$, it does not go below 1 until $n_b=25$, but for $n_b=16$ it is quite close and this is chosen. Here, $\chi^2=26.106$ against $\chi_{15,0.95}^2=24.996$. Figure~\ref{rjpm_tidesurgeint_honning} shows the surge data before and after transformation with 16 tidal bands.
 
Estimated parameters are $\hat{\mu}_{s^*}=2.943 (0.274)$, $\hat{\sigma}_{s^*}=1.617 (0.124)$, $\hat{\gamma}_{s^*}=-0.041 (0.049)$, $\hat{\theta}=0.085$ and $\hat{\theta}_{s^*}=0.043$. 

Return levels and confidence intervals for Honningsv�g are shown in Table~\ref{rjpm_table}. 

 
%
\begin{figure}[!htb]
\centering
\subfigure[Tide levels and concurrent surges above the 99.75\% quantile.]{
\includegraphics[scale=0.35]{EWL/Figures/rjpm_tidesurgeint_honning.pdf}%
%\caption{Tide levels and concurrent surges above the 99.75\% quantile for Honningsvåg.}%
\label{rjpm_tidesurgeint_honning}
}
\subfigure[Ratio between $\chi^2$ test statistic and table value $\chi_{n_b-1,0.95}^2$.]{
\includegraphics[scale=0.35]{EWL/Figures/rjpm_chiratio_honning.pdf}%
%\caption{The interaction functions $\hat{a}_1(X)$, $\hat{a}_2(X)$ (red) for 5 and 10 tidal bands in Honningsvåg.}%
\label{rjpm_chiratio_honning}
}
\caption{Tide-surge interaction plots for Honningsv�g.}
\label{rjpm_tidesurge_honning}
\end{figure}


\begin{table}[!htb]
\begin{center}
	\caption{Return level estimates and 95\% confidence intervals for the RJP model. }
	\label{rjpm_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 & 194.5 & (184.1, 207.6) & 318.2 & (314.6, 323.5) & 332.4 & (331.3, 335.6) \\ 
10& 206.7 & (192.5, 225.2) & 325.4 & (321.1, 336.1) & 341.9 & (339.1, 347.9) \\ 
20& 217.9 & (200.0, 242.5) & 333.0 & (327.0, 354.0) & 350.9 & (346.1, 360.4)) \\ 
100& 242.0 & (215.0, 282.9) & 352.5 & (339.3, 422.4) & 370.6 & (360.2, 390.7) \\ 
200& 251.8 & (220.7, 300.7) & 362.1 & (344.1, 466.5) & 378.7 & (365.6, 404.4) \\
\hline
\end{tabular}
\end{center}
\end{table}