


\section{Introduction}

The title of this chapter is deliberately chosen to focus on the bivariate case instead of the general multivariate. The reason is mainly one of expediency, because the general multivariate case would easily embroil us in the necessity to roll out a heavy machinery of notation without contributing to a deeper understanding of the issues involved. For a discussion of the general multivariate case, the reader may consult the book by \cite{begosete04}.
The extension of extreme value statistics from the univariate to the multivariate case meets with several challenges. 
First of all, there is no direct simple generalization of the univariate extreme value types theorem to the multivariate case, and in particular, this also applies to the bivariate case. 

Developed on the basis of Gumbel's logistic and mixed models \citep{Gumbel1960,Gumbel1960b,Gumbel1961,GumbMusta1967}), 
the later results on possible asymptotic bivariate extreme value distributions became, in a sense, too general, which poses severe problems for practical applications. Significant efforts have been made to model and estimate a function which describes the dependence structure between extreme components, cf. e.g. \cite{CoTa1991,CoTa1994}.
However, there are no precise estimation tools that allow us to decide on the joint distribution 
of the bivariate extremes from a given set of bivariate data. 
Of course, the marginal data sets can be used to derive estimates of the marginal extreme value distributions,
as in \cite{Zachary1998273,dehaanjohn98}, but the joint distribution is still a long way off.

A popular method of trying to cope with the problem of bivariate extremes is to adopt a copula to represent the joint distribution structure. This copula is then usually combined with asymptotic extreme value distributions to represent the marginal distributions, typically of the GEV type \citep{cole01}.  
For this purpose, a range of different copulas have been proposed \citep{Tawn1988,waalgelder05}, see also \cite{cahabasa05} for a compendious treatment of bivariate copulas. 
Even in the case of the bivariate extreme value copula \citep{pikands81,BalakLai2009},
due to the properties of the dependence function, generally speaking, there are an infinite number of models.
Therefore, the main problem with this approach is that it is rather ad hoc. 
That is, there seems to be no appropriate theoretical justification for choosing one particular copula over the other. 

%Efforts have also been invested in trying to extend the peaks-over-threshold (POT) method to the multivariate case \cite{Kaufmann1995235}. This has %not yet resulted in a method of the same prediction capabilities as the univariate POT. 

It is therefore of considerable interest to note that the ACER method can easily be extended to several dimensions, in particular, to two \citep{naes11p01}. By this fact, a vehicle is obtained for providing a nonparametric statistical estimate of the bivariate extreme value distribution inherent in a bivariate time series. It will be seen that the bivariate ACER function is able to cover both spatial and temporal dependence characteristics of the given time series. Thus, it covers all simultaneous and non-simultaneous extreme events.
From a practical point of view, this makes it possible to investigate the true behavior of the bivariate extreme value distribution for a particular case, and at the same time check the validity of the proposed copula models for bivariate extremes. 

As a first effort in investigating the functional representation of the empirically estimated bivariate ACER surface, the bivariate extreme value copula approach will be adopted. Specifically, the Asymmetric logistic and Gumbel logistic models are used, combined with asymptotically consistent marginal extreme value distributions based on the univariate ACER functions.
Since the univariate ACER functions have proved to portray accurately the marginal tail behaviour, this will offer an opportunity to verify to some extent the viability of copula models to capture the dependence structure of bivariate extreme value distributions.

The performance of the bivariate ACER method will be illustrated by application to two measured bivariate time series. The first consists of simultaneous wind speed data measured at two stations off the coast of Norway. The second time series consists of simultaneous wind speed and wave height data measured at an offshore location in the North Sea. For some other examples on the application of the bivariate ACER method, see e.g. \cite{ganaka19aor,xiaosenxu22}



\section{Componentwise Extremes}\label{ce}\index{componentwise extremes}\index{bivariate extreme value distributions}

Assume that $(X_1, Y_1), (X_2, Y_2), \ldots$ is a sequence of iid bivariate random variables having distribution function $F(x, y)$. Let $M_{x,N} = \max_{1 \le i \le N}{X_i}$ and $M_{y,N} = \max_{1 \le i \le N}{Y_i}$. Then $M_N = (M_{x,N}, M_{y,N})$ is a vector of componentwise maxima, where the index $i$ for which the component time series $X_1, X_2, \ldots$ and $Y_1, Y_2, \ldots$ assume their extreme may be different between the two series. Hence, $M_N$ does not necessarily correspond to an observed vector of the original time series. 

The search for the limiting forms of the bivariate extreme value distributions follows more or less the pattern of the univariate case by studying $M_N$ as $N \to \infty$. By way of a first observation, the marginal distributions of the asymptotic bivariate extreme value distribution should by necessity be given by the asymptotic univariate extreme value distributions. Since the marginals are known, the representations may, in fact, be simplified a little by assuming that both $X_i$ and $Y_i$ have the standard Fr\'{e}chet distribution $F(z) = \exp(- 1/z)$, $z > 0$. By transformation of variables, any other marginal distribution can be obtained. With this specific choice of marginal, a special case of the GEV distribution with parameters $\mu = 0$, $\sigma = 1$ and $\gamma = 1$ is at hand. Now, it follows that $\Prob(M_{x,N} \le z) = \exp (-N/z)$, or, equivalently, $\Prob(M_{x,N}/N \le z) = \exp (-1/z)$, $z > 0$. The same result applies to $M_{y,N}$. Hence, to obtain standard univariate results for each margin, the re-scaled vector,  
\begin{equation}\label{eq:bv.1}
M_N^* = (M_{x,N}/N, M_{y,N}/N),
\end{equation}
should be considered. 

The main result of this section is the following:

\bigskip
Let $M_N^* = (M_{x,N}^*, M_{y,N}^*)$ be defined by Eq. \eqref{eq:bv.1} where the $(X_i, Y_i)$ are iid random vetors with standard Fr\'{e}chet marginal distributions. Then if, 
\begin{equation}\label{eq:bv.2}
\Prob(M_{x,N}^* \le x, M_{y,N}^* \le y ) \to G(x, y), \:\: \mbox{as} \:\: N \to \infty,
\end{equation}
where $G$ is a non-degenerate distribution function, $G$ has the form,
\begin{equation}\label{eq:bv.3}
G(x, y) = \exp \{- V(x,y) \}, \:\: x > 0, \: y > 0,
\end{equation}
where 
\begin{equation}\label{eq:bv.4}
V(x, y) = 2 \int_0^1 \max \left( \frac{w}{x},\frac{1-w}{y} \right) dH(w),
\end{equation}
where $H$ is a distribution function on $[0,1]$ satisfying the mean value constraint
\begin{equation}\label{eq:bv.5}
\int_0^1 w \, dH(w) = 1/2.
\end{equation}

The distributions obtained as limits in Eq. \eqref{eq:bv.2} are called the class of bivariate asymptotic extreme value distributions. According to this result, this class is in a 1-1 correspondence with the set of distributions $H$ on $[0,1]$ satisfying Eq.~\eqref{eq:bv.5}. So, this rather remarkable result tells us that for any such distribution function on $[0, 1]$, a valid bivariate asymptotic extreme value distribution is obtained.  

Two simple examples may serve as illustration. Denote by $W$ a random variable with distribution function $H$. First, let $W$ have the two possible outcomes 0 and 1, and $\Prob(W = 0) = \Prob(W = 1) = 0.5$. Hence the distribution function $H$ has jumps at 0 and 1, so that $dH(0) = 0.5$ and $dH(1) = 0.5$. The condition given by Eq. \eqref{eq:bv.5} is satisfied, and $V(x, y) = x^{-1} + y^{-1}$ by Eq. \eqref{eq:bv.4}, so that the corresponding bivariate extreme value distribution becomes,
\begin{equation}\label{eq:bv.6}
G(x, y) = \exp \{- (x^{-1} + y^{-1}) \} = \exp \{- x^{-1}\}  \exp \{ - y^{-1}) \}, \:\: x > 0, \: y > 0,
\end{equation}
which clearly illustrates the case of two independent variables with standard Fr\'{e}chet marginals. The opposite case of two fully dependent variables is obtained by considering a degenerate random variable with unit mass at $W = 0.5$. The corresponding distribution function then has one jump at $w=0.5$. i.e. $dH(0.5)=1$. Eq. \eqref{eq:bv.5} is satisfied, and the corresponding bivariate extreme value distribution is,
\begin{equation}\label{eq:bv.7}
 G(x, y) = \exp \{- \max(x^{-1}, y^{-1}) \}, \:\: x > 0, \: y > 0,
\end{equation}
where the marginals are still standard Fr\'{e}chet.

To obtain the general situation for any GEV marginal, it is now only necessary to transform the marginals from standard Fr\'{e}chet to the required members of the GEV family. Specifically, by defining,
\begin{equation}\label{eq:bv.8}
\tilde{x} = \left\{ 1 + \gamma_x \Big( \frac{x - \mu_x}{\sigma_x} \Big) \right\}^{ 1/ \gamma_x} \quad \mbox{and} \quad  \tilde{y} = \left\{ 1 + \gamma_y \Big( \frac{y - \mu_y}{\sigma_y} \Big) \right\}^{ 1/ \gamma_y},
\end{equation}
it follows that the complete set of bivariate asymptotic extreme value distributions is determined by distribution functions of the form,
\begin{equation}\label{eq:bv.9}
G(x, y) = \exp \{- V(\tilde{x}, \tilde{y}) \},  
\end{equation}
provided $[1 + \gamma_x (x - \mu_x)/\sigma_x] > 0$ and $[1 + \gamma_y (y - \mu_y)/\sigma_y] > 0$, and where the function $V$ satisfies Eq. \eqref{eq:bv.4} for some distribution function $H$, satisfying Eq.~\eqref{eq:bv.5}. The marginal distributions are GEV with parameters $(\mu_x, \sigma_x, \gamma_x)$ and $(\mu_y, \sigma_y, \gamma_y)$, respectively.







\section{Bivariate ACER Functions}\label{baf}\index{bivariate ACER functions}

Consider a bivariate stochastic process $Z(t) = \big(X(t), Y(t)\big)$ with dependent component processes, which has been observed over
a time interval, $(0,T)$ say. Assume that the sampled values $(X_1, Y_1), \ldots,
(X_N, Y_N)$ are allocated to the (usually equidistant) discrete times $t_1, \ldots, t_N$ in $(0,T)$.
Our goal in this section is to construct a methodology that allows us to accurately determine empirically the joint distribution function of the extreme value vector $M_N = \left(M_{x,N}, M_{y,N}\right)$, where $M_{x,N} = \max\left\{ X_j\,; j = 1, \ldots, N \right\}$, 
and with a similar definition of $M_{y,N}$.
Specifically, the goal is to estimate $P(\xi, \eta) = \Prob\left(M_{x,N} \le \xi, M_{y,N} \le \eta\right)$
accurately for large values of $\xi$ and $\eta$. 

%In the following we outline the implementation of a sequence of approximations based on
%conditioning, where the first is a one-step memory approximation, which may be considered a Markov-like approximation. This
%approximation concept is described in \cite{naes85aor,naes90jsv}. However, it is emphasized that it is not a Markov chain approximation.

For notational convenience, it is expedient to introduce the non-exceedance event $\mathcal{C}_{kj}(\xi, \eta) = \big\{ X_{j-1} \leq \xi, Y_{j-1} \leq \eta, \ldots, X_{j-k+1} \leq \xi, Y_{j-k+1} \le \eta \big\}$ for $2 \le k \le j \le N+1$. Then, from the definition of $P(\xi, \eta)$ it follows that,
%% 
\begin{align}\label{eq:M2}
&P(\xi, \eta) = \, \Prob \big( \mathcal{C}_{N+1,N+1}(\xi, \eta) \big) \nonumber \\
&= \, \Prob \big( X_N \leq \xi, Y_N \leq \eta \, |\, \mathcal{C}_{NN}(\xi, \eta) \big) \cdot \, \Prob\big(\mathcal{C}_{NN}(\xi, \eta) \big) \nonumber \\
&\vdots \nonumber \nonumber \\
&= \prod_{j=2}^N \, \Prob \big( X_j \leq \xi, Y_j \leq \eta \,|\, \mathcal{C}_{jj}(\xi, \eta) \big) \cdot \, \Prob \big( \mathcal{C}_{22}(\xi, \eta) \big).
\end{align}
%%

Following the pattern of the derivations of the univariate ACER functions in Chapter~\ref{chap:acer}, it will be shown in Section~\ref{ApA}, based on Eq.~\eqref{eq:M2} and the properties of conditional probability, that a sequence of approximations may be introduced which converges to the target distribution $P(\xi, \eta)$.
In practice, following the derivations in Section~\ref{ApA}, it is therefore assumed that the following representation applies for a suitably chosen $k$, and for large values of $\xi$ and $\eta$, cf. Eq.~\eqref{eq:M39},
%%
\begin{align}\label{eq:generalApprox}
P(\xi, \eta) \approx \exp \Big\{ - \sum_{j=k}^N \big( \alpha_{kj}(\xi; \eta) + \beta_{kj}(\eta; \xi) - \gamma_{kj}(\xi, \eta) \big)\Big\};  
\end{align}
%%
where $\alpha_{kj}(\xi; \eta) = \Prob \big( X_j > \xi \,|\,\mathcal{C}_{kj}(\xi, \eta) \big)$, $\beta_{kj}(\eta; \xi) = \Prob \big( Y_j > \eta \,|\,\mathcal{C}_{kj}(\xi, \eta) \big)$ and $\gamma_{kj}(\xi, \eta) = \Prob \big( X_j > \xi, Y_j > \eta\,|\,\mathcal{C}_{kj}(\xi, \eta) \big)$. Note that Eq.~\eqref{eq:generalApprox} applies equally well to stationary and nonstationary time series. This is due to the fact that the possible time dependence of the conditional exceedance probabilities $\alpha_{kj}(\xi; \eta)$, $\beta_{kj}(\eta; \xi)$ and $\gamma_{kj}(\xi, \eta)$ has been retained, which is reflected in the presence of the time parameter $j$.  

From Eq.~\eqref{eq:generalApprox} it emerges that for the estimation of the bivariate extreme value distribution, 
it is necessary and sufficient to estimate 
the sequence of functions $\big\{\alpha_{kj}(\xi; \eta) + \beta_{kj}(\eta; \xi) - \gamma_{kj}(\xi, \eta)\big\}_{j=k}^{N}$. 
To get a more compact representation, it is expedient to introduce the concept of a $k$'th order bivariate average conditional exceedance rate (ACER) function as follows, cf. Eq.~\eqref{eq:bA0},
%%
\begin{align}\label{eq:A0}
& \mathcal{E}_k (\xi, \eta)\,= \,\frac{1}{N-k+1}\,\sum\limits_{j=k}^N \big(\alpha_{kj}(\xi; \eta) + \beta_{kj}(\eta; \xi) - \gamma_{kj}(\xi, \eta)\big) ; \:\: k=1,\,2,\ldots
\end{align}
%%
Hence, when $N \gg k$, and for large values of $\xi$ and $\eta$,
\begin{align}\label{eq:A4}
P(\xi, \eta) \approx \exp \left\{\, - \,(N-k+1)\, \mathcal{E}_k (\xi, \eta) \right\} .
\end{align}
%%

The numerical estimation of the bivariate ACER function for the observed stationary or 
non-stationary time series consists in counting of the appropriate exceedance events.
The estimation procedure is derived in more details in Section~\ref{ApB}.




\section{Functional Representation of the Empirically Estimated Bivariate ACER Functions}\index{bivariate ACER functions!empirical estimation}

From the definition of $\mathcal{E}_k (\xi, \eta)$ follows that the product $\mathcal{E}_k (\xi, \eta)\cdot(N-k+1)$ represents 
the expected number of the bivariate observations $Z_j = (X_j,\, Y_j)$ such that their components exceed corresponding levels $\xi$ and $\eta$  
(both simultaneous and non-simultaneous) and follow after at least $k-1$ previous simultaneous non-exceedances. 
Therefore the bivariate ACER function $\mathcal{E}_k(\xi, \eta)$ is able to capture the temporal and spatial dependence structure
of the considered bivariate time series. Moreover, as is discussed in Section~\ref{ApA}, in practice the existence of an effective 
$k_e \ll N$ such that $P(\xi, \eta) = P_{k_e}(\xi, \eta) = \exp \left\{\, - \,(N-k_e+1)\,\mathcal{E}_{k_e}(\xi, \eta)\right\}$ can be assumed.

This provides a means to obtain high quantiles of the bivariate extreme value distribution. The joint $T$-year return period contour associated with the event that either $M_{x,N}$ or $M_{y,N}$ or both is exceeded, that is, $ \big\{\big(M_{x,N} > \xi^{\scriptscriptstyle T}\big) \cup \big(M_{y,N} > \eta^{\scriptscriptstyle T}\big) \cup \big(M_{x,N}>\xi^{\scriptscriptstyle T} \cap M_{y,N}>\eta^{\scriptscriptstyle T}\big)  \big\} $, is represented by,
%%
\begin{equation}\label{eq:Ty}
1 - F^{1yr}(\xi^{\scriptscriptstyle T},\, \eta^{\scriptscriptstyle T}) = \frac{1}{T}, 
\end{equation}
%%
where $F^{1yr}(\xi, \eta)$ is the joint distribution function of the annual maxima.
Assuming that the duration of the observation period of the bivariate process $Z(t)$ is $n_y$ years, then, with $k = k_e$,
%% 
\begin{equation}\label{eq:Foneyear}
F^{1yr}(\xi, \eta) = \exp \left\{\, - \,\frac{N-k+1}{n_y}\,\mathcal{E}_{k}(\xi, \eta)\right\}.
\end{equation}
%%
From Eqs.~(\ref{eq:Ty}) and~(\ref{eq:Foneyear}) follows that the joint $T$-year return levels 
$(\xi^{\scriptscriptstyle T},\,\eta^{\scriptscriptstyle T})$ are obtained as solution of the implicit equation:
%% 
\begin{equation}\label{eq:E_return}
\mathcal{E}_{k}(\xi^{\scriptscriptstyle T},\,\eta^{\scriptscriptstyle T}) = -\,\log\left(1-\frac{1}{T}\right)\frac{n_y}{N-k+1}.
\end{equation}
%%

It is evident, that the empirically estimated $k$-th order bivariate ACER rarely provides enough information for estimation of high quantiles of the joint extreme value distribution.
In addition, the exact behavior of the bivariate ACER as a continuous function of two variables cannot be decided using available statistical data.
Therefore, the sub-asymptotic functional form of the ACER surface ${\mathcal{E}}_k (\xi, \eta)$ can possibly be obtained approximately by the copula representation of a bivariate extreme value distribution.

From the result by \cite{Sklar59}, for any pair of random variables $(X,\, Y)$ with marginal 
distribution functions $F_{x}(\xi)$ and $G_{y}(\eta)$, the joint distribution function
$H_{xy}(\xi,\,\eta) = \Prob(X\le\xi,\,Y\le\eta)$ can be presented by a bivariate copula\index{bivariate copula} $C(u,v)$ as  follows:
$H_{xy}(\xi,\,\eta) = C\big(F_{x}(\xi),\ G_{y}(\eta)\big)$,  cf. e.g. \cite{Nelsen2006} or \cite{BalakLai2009}.
This result applies to any bivariate extreme value distribution as well.

Considering the above and using the result of \cite{pikands81}, any bivariate extreme value distribution $H_{xy}(\xi,\,\eta)$ with marginal univariate extreme value distributions $F_{x}(\xi)$ and $G_{y}(\eta)$ is given by the formula,
%%
\begin{align}\label{biEVD}
% C\left(u,\,v\right) = \exp\left\{\log\big(u v\big)\mathcal{D}\left(\frac{\log(u)}{\log(u v)} \right)  \right\}\,; \:\: 0<u<1,\,0<v<1\,,
H_{xy}(\xi,\,\eta) = \exp  \Bigg\{ \log\big(F_{x}(\xi)G_{y}(\eta)\big)
 \cdot\mathcal{D}\left(\frac{\log\big(F_{x}(\xi)\big)}{\log\big(F_{x}(\xi)G_{y}(\eta)\big)} \right) \Bigg\},
\end{align}
%%
where the Pickands dependence function\index{Pickands dependence function} $\mathcal{D}(\cdot)$ is a convex function and satisfies $\mathcal{D}(x):[0,1]\longmapsto [\,\max(x,\,1-x),\ 1]$, cf. \cite{jaworski2010copula}.

We assume that asymptotically consistent marginal extreme value distributions $ F_{x}(\xi)$ and $G_{y}(\eta)$ are represented by the corresponding univariate ACER functions, that is,
%%
\begin{align}\label{Marginals}
\begin{aligned}
F_{x}(\xi) & \approx \exp\big\{-(N-k+1) \varepsilon^{x}_{k}(\xi) \big\}\,, \:\: \xi \ge \xi_1 \,,\\
G_{y}(\eta)& \approx \exp\big\{-(N-k+1) \varepsilon^{y}_{k}(\eta) \big\}\,, \:\: \eta \ge \eta_1 \,,
\end{aligned}
\end{align}
%%
where the sub-asymptotic functional form of the univariate ACER function is represented as 
$\varepsilon^{x}_{k}(\xi) = q^{x}_k \,\exp \{- a^{x}_k (\xi - b^{x}_k )^{c^{x}_k} \}$ with a similar definition 
of $\varepsilon^{y}_{k}(\eta)$, cf. Chapter~\ref{chap:acer}.

Now, substituting Eq.~\eqref{Marginals} into Eq.~\eqref{biEVD} the following representation of the bivariate extreme value distribution applies:
%%
\begin{align}\label{biEVDviauniACER}
H_{xy}(\xi,\,\eta) =  \exp\Bigg\{\, - \,(N-k+1)\, \Big(\varepsilon^{x}_{k}(\xi) + \varepsilon^{y}_{k}(\eta)\Big) 
 \cdot\mathcal{D}\left(\frac{\varepsilon^{x}_{k}(\xi)}{\varepsilon^{x}_{k}(\xi) + \varepsilon^{y}_{k}(\eta)} \right) \Bigg\}.
\end{align}
%%  

On the other hand, as it has been discovered before in Eq.~\eqref{eq:A4}, the bivariate extreme value distribution can be expressed through the bivariate ACER function: $H_{xy}(\xi, \eta) = \exp \left\{\, - \,(N-k+1)\, \mathcal{E}_k (\xi, \eta) \right\}$.
Thereby, the functional form of the bivariate ACER surface can possibly be obtained by:
%%
\begin{align}\label{biACERviaD}
\begin{aligned}
\mathcal{E}_k (\xi, \eta) =  \Big(\varepsilon^{x}_{k}(\xi) + \varepsilon^{y}_{k}(\eta)\Big)
\mathcal{D}\left(\frac{\varepsilon^{x}_{k}(\xi)}{\varepsilon^{x}_{k}(\xi) + \varepsilon^{y}_{k}(\eta)} \right).
\end{aligned}
\end{align}
%% 
Consequently, our aim now is to find the dependence function $\mathcal{D}(\cdot)$ that would provide
optimal fit of the parametrical surface defined by Eq.~\eqref{biACERviaD} to the empirical bivariate ACER $\mathcal{\hat{E}}_k (\xi, \eta)$.
% mean square error function 

Subject to the form of the dependence function $\mathcal{D}(\cdot)$, different 
parametric differentiable and non-differentiable models can be considered. 
By setting $\mathcal{D}(x) = \theta x^2-\theta x + 1$, where $0\le\theta\le1$,
the {Type A} bivariate extreme value distribution\index{bivariate extreme value distribution!type A}, or Gumbel mixed model\index{Gumbel mixed model}\index{bivariate extreme value distribution!Gumbel mixed model}, is obtained, cf. e.g. \cite{Gumbel1960,Gumbel1960b}, \cite{GumbMusta1967}. 
Another differentiable model is acquired by setting $\mathcal{D}(x) = [x^m +(1-x)^m]^{1/m}$ for $m\ge1$. 
This is the {Type B} distribution\index{bivariate extreme value distribution!type B} or Gumbel logistic model\index{bivariate extreme value distribution!Gumbel logistic model} \citep{Gumbel1961,Hougaa1986}. 
The functional form of the bivariate ACER surface\index{bivariate extreme value distribution!bivariate ACER surface} in the Gumbel logistic case becomes, % ${\mathcal{A}}_k(\xi, \eta)$ 
%%
\begin{align}\label{gl}
\begin{aligned}
{\mathcal{G}}_k (\xi, \eta)\, = \,\Big[&\big(\varepsilon^{x}_{k}(\xi)\big)^m + \big(
\varepsilon^{y}_{k}(\eta)\big)^m \Big]^{\frac{1}{m}}.
\end{aligned}
\end{align}

 
The {Type C} distribution\index{bivariate extreme value distribution!type C}, also known as the biextremal model\index{biextremal model} \citep{Tiago84}, can be considered as an example of a non-differentiable model.
The dependence function in this case is $\mathcal{D}(x) = \max(x,\ 1-\theta x)$ for $0\le\theta\le1$.

In the literature, differentiable models have usually been of most interest, and they have been used to analyze bivariate environmental events. \cite{Yue1999} and \cite{Yue2000}, apply the Gumbel mixed model to rainfall data in order to provide storm frequency analysis. \cite{Yue2001b,Yue2001a} has also studied the Gumbel logistic model with application to flood peak -- flood volume pair of bivariate data.

In the work by \cite{YueWang04}, a comparison between the Gumbel mixed and the Gumbel logistic models has been made. The authors argued that both models are appropriate and give similar estimates of the joint distribution of two Gumbel-distributed random variables whose Pearson product moment correlation coefficient is: $0\leq \rho\leq 2/3$. When $\rho > 2/3$ the Gumbel mixed model cannot be applied, see also \cite{Tiago82}. For this reason, it was decided to consider the Gumbel logistic model as one that fits better the objectives of the present work. 

In the dependence function $\mathcal{D}(\cdot)$ for the Gumbel logistic model, \cite{Tawn1988} added extra parameters $\phi$ and $\theta$ to get further flexibility. This leads to the Asymmetric logistic model\index{bivariate extreme value distribution!asymmetric logistic model}, which sets
$\mathcal{D}(x) = [\phi^m x^m + \theta^m(1-x)^m]^{1/m} + (\theta-\phi)x +1 -\theta\ $ with $0\le\theta\le1,\ 0\le\phi\le1,\ m\ge1$. The functional form of the bivariate ACER surface in the Asymmetric logistic case is obtained as, % ${\mathcal{A}}_k(\xi, \eta)$ 
%%
\begin{align}\label{al}
\begin{aligned}
{\mathcal{A}}_k (\xi, \eta)\, = \,\Big[\big(\phi\varepsilon^{x}_{k}(\xi)\big)^m + \big(\theta
\varepsilon^{y}_{k}(\eta)\big)^m \Big]^{\frac{1}{m}} 
 + (1-\phi)\varepsilon^{x}_{k}(\xi) + (1-\theta)\varepsilon^{y}_{k}(\eta).
\end{aligned}
\end{align}
%%
The optimal parameters $m^{*}$ for $\mathcal{G}_k (\xi, \eta)$ and $\theta^{*},\ \phi^{*}$ and $m^{*}$ for $\mathcal{A}_k (\xi, \eta)$, can be found by minimizing a mean square error function. This is discussed only for the case of $\mathcal{A}_k (\xi, \eta)$, since the case of $\mathcal{G}_k (\xi, \eta)$ then follows easily. Specifically, the mean square error function for $\mathcal{A}_k (\xi, \eta)$ is defined as,
%%
\begin{align}\label{MSE}
F(m,\theta,\phi) = \sum\limits_{j=1}^{N_{\eta}}\sum\limits_{i=1}^{N_{\xi}} {w}'_{ij}\Big(\log \hat{\mathcal{E}}_k (\xi_{i}, \eta_{j}) - \log {\mathcal{A}}_k (\xi_{i}, \eta_{j})\Big)^2,
\end{align}
%%
where $N_{\xi}, N_{\eta}$ are numbers of levels $\xi$ and $\eta$, respectively, at which the ACER function have been empirically estimated, and 
$w'_{ij} = w_{ij} / \sum\sum w_{ij}$ with,
%%
\begin{align}
w_{ij} = \left(\log {CI}^{+}(\xi_{i}, \eta_{j}) - \log CI^{-}(\xi_{i}, \eta_{j})\right)^{-2}\,,
\end{align}
%%
denoting normalized weight factors that put more emphasis on the more reliable estimates.

The constrained optimization problem with the objective function $F$
defined in Eq.~\eqref{MSE}, is written as,
%%
\begin{align}\label{eq:17}
\left\{
\begin{aligned}
& F(m,\theta,\phi) \rightarrow \min\,,\\ %\min\limits_{\{q,b,a,c\}} 
& \{m,\theta,\phi\}\ \in\ \mathcal{S} \,,
\end{aligned} 
\right.
\end{align}
%%
with a constraints domain,
%%
\begin{align}\label{eq:18}
\mathcal{S}=\big\{ \{m,\theta,\phi\} \in \mathbb{R}^3\ \big\vert\ \, \theta,\,\phi\ \in[0,\ 1];\ m\in[1,\, +\infty) \big\}.
\end{align}
%% 



\section{Numerical Examples}


\subsection{Wind speed measured at two adjacent weather stations}\label{WS2}

In this section, the results obtained by \cite{naka15jweia} will be discussed. The simultaneous wind speed data measured along the Norwegian coast at Sula\index{wind speed measurements!Sula} and Nord{\o}yan Fyr\index{wind speed measurements!Nord{\o}yan Fyr} weather stations, were analyzed to obtain numerical estimates of bivariate extreme wind speeds. 
Figure~\ref{fig:ws1} shows the geographical locations of the measurement sites. The hourly maximum of the three seconds wind gust (10 meters above the ground) were recorded during 13 years (1999 - 2012).

%%
\begin{figure}[htb!]
	\centerline{\includegraphics[width=2.5in]{Figures/map_sula_nordoyan.pdf}}
	\caption{Map of a part of Norway with the marked weather stations: A -- Sula station, B -- Nord{\o}yan Fyr station}
	\label{fig:ws1}
\end{figure}
%%
Figure~\ref{fig:ws2} demonstrates the plot of the observed data. This plot reveals a rather strong dependence between the two time series.
%%
%%
\begin{figure}[htb!]
	\centerline{\includegraphics[width=2.75in]{Figures/biTS_all.pdf}} % scale=0.5
	\caption{Coupled observations of wind speed data observed at the Sula station ($\xi$ axis) and at the Nord{\o}yan Fyr station ($\eta$ axis).}
	\label{fig:ws2}
\end{figure}
%%


The Pearson product-moment correlation coefficient is found to be $\rho = 0.73$.
The Kendall's rank correlation coefficient is $\tau = 0.5$, while Spearman's $\rho$ equals $0.68$, which also indicates a nonlinear relationship between Sula and Nord{\o}yan wind speeds.

It was decided to divide the data series into 13 one-year records for the analysis. 
By this, also the standard deviation of the ACER function estimates can be calculated fairly accurately. 


The univariate ACER functions were estimated first, 
using the Matlab-based standalone downloadable application \citep{ACERmanual}.
In Figures~\ref{fig:ws3}~-~\ref{fig:ws4} the sequence of $\hat{{\varepsilon }}_1\,\ldots\,\hat{{\varepsilon }}_{96}$ are plotted versus different wind speed levels. Both figures reveal 
that there is significant temporal dependence between consecutive data. 
It is also seen that this dependence effect is largely accounted for by $k = 24$ since there is a marked degree of convergence in the tail of $\hat{{\varepsilon }}_k$ for $k \ge 24$ in both cases.
Here, for $k = 96$, which corresponds to conditioning on data
recorded up to 4 days earlier, $\hat{{\varepsilon }}_{96}$ is considered to represent the final converged results, since $\hat{{\varepsilon }}_{96}\approx\hat{{\varepsilon }}_{k}$ for $k>96$
in the tail. Therefore, there is no need to consider conditioning of
an even higher order than $96$. So, effectively, $k_e = 96$ for the recorded data. Also note that 4 days is a typical separation of wind speed data adopted in the declustering process to achieve independence between the data used in e.g. a peaks-over-threshold analysis. Figures~\ref{fig:ws3}~-~\ref{fig:ws4} 
also demonstrate that for extreme value estimation, $\hat{{\varepsilon}}_1$ 
can be used since the ACER functions all coalesce in the far tail. 
This makes it possible to choose $k = 1$, 
which makes much more data available for estimation, with a possible 
reduction of uncertainty in estimation as a result. 
%%
\begin{figure}[htb!]
	\centerline{\includegraphics[width=2.75in]{Figures/SULA_all_ACER.pdf}}
	\caption{ACER estimates for different degrees of conditioning. Wind speed data from the Sula station, cf. Figure~\ref{figP07}.}
	\label{fig:ws3}
\end{figure}
%%
%%
\begin{figure}[htb!]
	\centerline{\includegraphics[width=2.75in]{Figures/Nordoyan_all_ACER.pdf}}
	\caption{ACER estimates for different degrees of conditioning. Wind speed data from the Nord{\o}yan Fyr station.}
	\label{fig:ws4}
\end{figure}
%%

The sequence of estimated bivariate ACER surfaces $\hat{\mathcal{E}}_k (\xi, \eta)$ is shown in Figure~\ref{fig:ws5}. Matlab programs for ACER 2D analyses are available, cf. \cite{ACER2D14}.
$\hat{\mathcal{E}}_k (\xi, \eta)$ with $k = 1$ is the uppermost. As it is seen from the figure, the cross-section of the surfaces at a high value of the wind speed level $\eta$ gives the univariate ACER functions of the wind speed data from the Sula station, while the cross-section at a high level of $\xi$ represents the univariate ACER of the time series from the Nord{\o}yan Fyr.

Due to the observed convergence, the ACER surface for $k = 96$ is very close to the surface obtained by taking the logarithm of the exact bivariate extreme value distribution.
%%
\begin{figure}[htb!]
	\centering
	\includegraphics[width=3in]{Figures/surfaces_wind.pdf}
	\caption{Bivariate ACER surface estimates for different degrees of conditioning. $\hat{\mathcal{E}}_k (\xi, \eta)$ surfaces are plotted on a logarithmic scale.}
	\label{fig:ws5}
\end{figure}
%%

Parameters of the optimal Asymmetric logistic (AL) and Gumbel logistic (GL) surfaces are presented in Table~\ref{params}. 
%%
\begin{table}[h!]
	\caption{Optimal parameters of AL and GL fits.}
	\label{params}%
	\centering
	\begin{tabular}{c|c|c}
		\hline
		$k$      & AL    & GL \\
		\hline
		1 & $m_{A} = 2.44,\ \ \theta = 0.86,\ \ \phi = 0.92$ & $m_{G} = 2.01$ \\
		96 & $m_{A} = 4.53,\ \ \theta = 0.97,\ \ \phi = 0.97$ & $m_{G} = 3.87$ \\
		\hline
	\end{tabular}

\end{table}%
%%

Figures~\ref{fig:ws6}~-~\ref{fig:ws7} show the contour plots of the optimized Asymmetric logistic $\mathcal{A}_{k}(\xi, \eta)$ and the optimized Gumbel logistic ${\mathcal{G}}_{k}(\xi, \eta)$ fits to the data for $\hat{\mathcal{E}}_{k}(\xi, \eta)$ surface for $k = 1$ and $k = 96$, respectively. The contour lines of three surfaces are plotted for those levels of $\xi$ and $\eta$, 
where the bivariate ACER surface $\hat{\mathcal{E}}_{k}(\xi, \eta)$ have been empirically estimated. 
%%
\begin{figure}[htb!]
	\centering
	\includegraphics[width=2.75in]{Figures/V_levels.pdf}
	\caption{Contour plot of the empirically estimated $\hat{\mathcal{E}}_{1}$ surface, and the optimized Gumbel logistic $\mathcal{G}_{1}$ and
		optimized Asymmetric logistic $\mathcal{A}_{1}$ surfaces based on marginal univariate ACER. Boxes indicate levels on a logarithmic scale.}
	\label{fig:ws6}
\end{figure}
%%
%%
\begin{figure}[htb!]
	\centering
	\includegraphics[width=2.75in]{Figures/V_levels_96.pdf}
	\caption{Contour plot of the empirically estimated $\hat{\mathcal{E}}_{96}$ surface, and the optimized Gumbel logistic $\mathcal{G}_{96}$ and
		optimized Asymmetric logistic $\mathcal{A}_{96}$ surfaces based on marginal univariate ACER. Boxes indicate levels on a logarithmic scale.}
	\label{fig:ws7}
\end{figure}
%%

The figures reveal that the empirical bivariate ACER surface $\hat{\mathcal{E}}_{k}$ captures high correlation between the data, and so do the optimally fitted $\mathcal{G}_{k}$ and $\mathcal{A}_{k}$ surfaces. 
It is also seen that the behaviour of the estimated ACER surface in the case $k = 96$ in Figure~\ref{fig:ws7} affirms high uncertainty due to deficiency of data. However, the optimal surfaces $\mathcal{G}_{96}$ and $\mathcal{A}_{96}$ capture the statistical properties of the bivariate observations. 


It is noticeable that the level of agreement 
between the estimated bivariate ACER and both the optimized Asymmetric logistic and 
the Gumbel logistic surfaces is equally significant. 
%
Yet, it is also important to keep in mind that the empirical bivariate ACER $\hat{\mathcal{E}}_{k}$ 
is the only discrete surface. The \cite{MATLAB2009} built-in routine \textsf{contourc}
that has been used to obtain the figures, calculates the contour lines by
producing a regularly spaced grid determined by the dimensions of a surface.
Therefore, it evidently generates a certain spacing inaccuracy of the 
ACER $\hat{\mathcal{E}}_{k}$ surface level lines plot. 
In addition, the figures ascertain that the optimized Asymmetric logistic and 
Gumbel logistic surfaces conform at a level sufficient to affirm that they actually coincide.

Thereby, the optimized Gumbel logistic model with asymptotically consistent marginals obtained from the optimized univariate ACER, can be used as the parametric representative of the bivariate ACER surface estimated from the given data set.

Also, comparison of the contour lines of $\mathcal{G}_{1}$ and $\mathcal{G}_{96}$ that correspond to the same return period levels shows fairly good agreement in the tail considering the high uncertainty for the case $k = 96$. To highlight the results that would be obtained by adopting the common approach of assuming Gumbel marginal extreme value distributions combined with a suitable copula model, in Figure~\ref{fig:ws8} are plotted the 50 and 100 year return period levels obtained by using the Asymmetric logistic model with asymptotically consistent marginals obtained from the optimized univariate ACER marginals together with the corresponding return levels obtained by using the Gumbel logistic model with Gumbel marginals fitted by the method of moments. It is clear that the discrepancy is significant, which is primarily caused by the use of asymptotic Gumbel marginals. Finally, it is noted that the bivariate ACER methodology has been studied in more detail in an initial study on synthetic data with known extreme value distribution, and therefore a known 100-year return period level $\left(\xi^{100yr},\,\eta^{100yr}\right)$. To get an idea about the performance of the ACER method and the existing Gumbel logistic model with Gumbel marginals fitted by the method of moments, Monte Carlo simulations were carried out to produce 100 bivariate data samples. It was observed that the predicted 100-year return period levels were consistently better for the ACER method. 

\begin{figure}[htb!]
	\centering
	\includegraphics[width=2.75in]{Figures/T_levels_k1_50_100_bw.pdf}  % T_levels
	\caption{Contour plot of the return period levels for the optimized Asymmetric logistic $\mathcal{A}_{1}$ surfaces (solid line \textbf{---})
		and the Gumbel logistic model with the Gumbel marginals $\mathcal{G}_{MM}$ (dash-dot line \textbf{--}$\,\cdot\,$\textbf{--}). Boxes indicate return period levels in years.}
	\label{fig:ws8}
\end{figure}




\subsection{Wind speed and wave height measured at a North Sea weather station}\label{WSWH}

This example presents results obtained by \cite{naka15jomae}. Wind speed (WS - 3 hours mean [m/s]) and significant wave height (Hs - total sea [m]) data measured\index{wind speed measurements!North Sea weather station}\index{wave height measurements!North Sea weather station} in the Norwegian sea
at location N 65.29, E 7.32 were analyzed to obtain numerical estimates
of the bivariate extreme value distribution. Figure~\ref{fig01} shows the geographical position of the measurement site. The data were recorded during 54 years (1957 - 2011), 8 times per day (every three hours).
%%
\begin{figure}[htb!]
	\centerline{\includegraphics[width=2.5in]{Figures/mapka.pdf}}
	\caption{Map of a part of Norway with the marked location.}
	\label{fig01}
\end{figure}
%%

Figure~\ref{fig022} demonstrates the plot of the observed data. As it is seen from the plot, there is a rather strong dependence between the two time series.
%%
%%
\begin{figure}[htb!]
	\centerline{\includegraphics[width=2.75in]{Figures/TS.pdf}} % scale=0.5
	\caption{Coupled observations of wind speed data ($\xi$ axis) and significant wave height data (total sea, $\eta$ axis).}
	\label{fig022}
\end{figure}
%%


The Pearson product-moment correlation coefficient is found to be $\rho = 0.79$.
The Kendall's rank correlation coefficient is $\tau = 0.56$ and Spearman's $\rho$ is equal to $0.7$, which also indicates a nonlinear connection between WS and Hs.

 

It should be noted that the available bivariate observations have low accuracy.
This especially concerns the significant wave height data, where the graduating mark is 0.1 meter, and on average, there are 98 unique numerical values of the Hs data per year.
Obviously, for fairly accurate estimation of the bivariate ACER functions, more data wouuld be required. In order to get a good sample size of records (realizations), 
it was decided to divide the data series into 18 three-years records for the analysis. 
By this procedure, also the standard deviation of the ACER function estimates can be calculated fairly accurately. 


As in the previous example, the univariate ACER functions $\hat{{\varepsilon }}_k$ were estimated first, 
using the Matlab-based standalone downloadable application \citep{ACERmanual}.
In Figures~\ref{fig031}~-~\ref{fig032},  $\hat{{\varepsilon }}_k$ is plotted versus different levels of wind speeds and wave heights, respectively, for different values of $k$. 
From both figures, it is clearly seen that there is significant time dependence between WS observations, as well as between Hs data. 
It is also understood that this dependence effect is largely accounted for by $k = 16$ since there is a marked degree of convergence in the tail of $\hat{{\varepsilon }}_k$ for $k \ge 16$ in both cases.
% of wind speed and wave height data analysis.
Obviously, $k = 16$ corresponds to exceedances separated by at least two days of non-exceedances for 
three hours observations. For $k \ge 32$, which corresponds to four days declustered data, full convergence has been achieved. Figures~\ref{fig031}~-~\ref{fig032} 
also demonstrate that for extreme value estimation, $\hat{{\varepsilon }}_2$ 
can be used since the ACER functions for $k \ge 2$ all converge in the far tail. Again, this clearly demonstrates the power of an ACER function plot as a diagnostic tool to 
decide on the value of $k$ needed for extreme value estimation in a particular case. 
In spite of significant dependence effects for the WS and Hs data with lower magnitudes, for the extreme values this is largely absent. This makes it possible to choose $k = 2$, 
which makes much more data available for estimation, with a possible 
reduction of uncertainty in estimation as a result. 
%%

\begin{figure}[htb!]
	\centerline{\includegraphics[width=2.75in]{Figures/ACER_Ws.pdf}}
	\caption{ACER estimates for different degrees of conditioning. Wind speed data.}
	\label{fig031}
\end{figure}
%%
%%
\begin{figure}[htb!]
	\centerline{\includegraphics[width=2.75in]{Figures/ACER_Hs.pdf}}
	\caption{ACER estimates for different degrees of conditioning. Significant wave height (total sea).}
	\label{fig032}
\end{figure}
%%
%%

Figures~\ref{fig041}~-~\ref{fig042} show the plots of the optimized parametric fit to the 
data for $\hat{{\varepsilon }}_{k}$ for $k = 2$ for both time series.
In particular, 100-year return level value and its 95\% CI are estimated parametrically and plotted.
For the wind speed data, the optimal parameters are: $q = 0.05$, $b = 0.1$, $a = 1.9\cdot10^{-4}$, $c = 3.14$, 
while in case of Hs data, the parameters of the optimal curve are: $q = 0.04$, $b = -2.27$, $a = 0.02$, $c = 2.23$.
%%
\begin{figure}[htb!]
	\centerline{\includegraphics[width=2.75in]{Figures/Ws_fit.pdf}}
	\caption{Plot of $\hat{{\varepsilon }}_{2}(\xi)$ versus wind speeds $\xi$ on a logarithmic scale for the optimized parameter values; $\xi_{1} = 14.5$.}
	\label{fig041}
\end{figure}
%%
%%
\begin{figure}[htb!]
	\centerline{\includegraphics[width=2.75in]{Figures/Hs_fit.pdf}}
	\caption{Plot of $\hat{{\varepsilon }}_{2}(\eta)$ against wave heights $\eta$ on a logarithmic scale for the optimized parameter values; $\eta_{1} = 4.5$. }
	\label{fig042}
\end{figure}
%%

Figure~\ref{fig05} demonstrates the empirically estimated bivariate ACER surfaces $\hat{\mathcal{E}}_k (\xi, \eta)$
for different values of $k$ on a logarithmic scale. $\hat{\mathcal{E}}_k (\xi, \eta)$ with $k = 1$ is the uppermost. As it is seen from the figure, the cross-section of the surfaces at a high level of the wave height $\eta$ gives the univariate ACER functions of the wind speed data, while the cross-section at a high wind speed level represents the univariate ACER of the Hs time series, respectively.

The same arguments as in the univariate case are applied to make the decision about the bivariate ACER surface to be used in the analyses. That is, as long as the surfaces for $k \geq 2$ all converge in the tail and estimation of $\hat{\mathcal{E}}_2 (\xi, \eta)$ is more accurate due to availability of more data, we would choose the surface with the degree of conditioning $k = 2$. 
%%
\begin{figure}[htb!]
	\centering
	\includegraphics[width=3in]{Figures/surfaces_wh.pdf}
	\caption{Bivariate ACER surface estimates for different degrees of conditioning. $\hat{\mathcal{E}}_k (\xi, \eta)$ surfaces are plotted on a logarithmic scale.}
	\label{fig05}
\end{figure}
%%

The optimal parameters of the Asymmetric logistic fit were found to be: $m_{AL} = 7$, $\theta = 1$ and $\phi = 0.91$. 
Since the additional parameters $\theta$ and $\phi$ are close to one, it would seem reasonable to consider a Gumbel logistic model.
The same optimization procedure as in the Asymmetric logistic case, was applied to obtain the optimal dependence parameter of the Gumbel logistic copula $m_{GL} = 4.78$.

Figures~\ref{fig07}~-~\ref{fig08} show the contour plots of the optimized Asymmetric logistic fit $\mathcal{AL}_{2}(\xi, \eta)$ to the data for the $\hat{\mathcal{E}}_{2}(\xi, \eta)$ surface and also the contour plots of the optimized Gumbel logistic surface $\mathcal{GL}_{2}(\xi, \eta)$.

In Figure~\ref{fig07}, contour lines of the three surfaces are plotted for those levels of $\xi$ and $\eta$, 
where the bivariate ACER surface $\hat{\mathcal{E}}_{2}(\xi, \eta)$  have been empirically estimated. Contour lines that correspond to the return period levels are presented in Figure~\ref{fig08}.

The figures reveal that the empirical bivariate ACER surface $\hat{\mathcal{E}}_{2}$ captures high correlation between the data and so do the optimally fitted $\mathcal{GL}_{2}$ and $\mathcal{AL}_{2}$ surfaces. Note that the contour lines of the bivariate ACER surface of fully correlated data would show up as lines that consist of only horizontal and vertical line segments. This happens because for such data, $H_{xy}(\xi,\,\eta) = \min\big(F_{x}(\xi),\ G_{y}(\eta)\big)$, which implies that ${\mathcal{E}}_k (\xi, \eta) = \max\big(\varepsilon^{x}_{k}(\xi), \varepsilon^{y}_{k}(\eta)\big)$, cf. Eq.~\eqref{Marginals}.  It is seen that the level of agreement between the estimated bivariate ACER and optimized Asymmetric logistic and Gumbel logistic surfaces is equally significant. Thereby, the optimized Gumbel logistic model with asymptotically consistent marginals obtained from the optimized univariate ACER, can be used as the parametric representative of the bivariate ACER surface estimated from the given data set. The Matlab programs for ACER 2D analyses from \cite{ACER2D14} were used to obtain the results presented in this section.


%%
\begin{figure}[htb!]
	\centering
	\includegraphics[width=2.75in]{Figures/V_levels_AL_GL.pdf}
	\caption{Contour plot of the empirically estimated $\hat{\mathcal{E}}_{2}(\xi, \eta)$ surface ($\bullet$),
		optimized Gumbel logistic $\mathcal{GL}_{2}(\xi, \eta)$ ($\circ$) and optimized Asymmetric logistic $\mathcal{AL}_{2}(\xi, \eta)$ (\textbf{---}) surfaces. 
		Boxes indicate levels on a logarithmic scale.}
	\label{fig07}
\end{figure}
%%
\begin{figure}[htb!]
	\centering
	\includegraphics[width=2.75in]{Figures/T_return_AL_GL.pdf}
	\caption{Contour plot of the return period levels for $\hat{\mathcal{E}}_{2}(\xi, \eta)$ surface ($\bullet$),
		optimized Asymmetric logistic $\mathcal{AL}_{2}(\xi, \eta)$ (\textbf{---}) and optimized Gumbel logistic
		$\mathcal{GL}_{2}(\xi, \eta)$ ($\circ$) surfaces. Boxes indicate return period levels in years.}
	\label{fig08}
\end{figure}
%%



\section{Appendix 1: The Sequence of Conditioning Approximations}\label{ApA}

Consider a bivariate stochastic process $Z(t) = \big(X(t), Y(t)\big)$ with dependent component processes, which has been observed over
a time interval, $(0,T)$ say. Assume that the sampled values $Z_1 = (X_1, Y_1), \ldots,
Z_N = (X_N, Y_N)$ are allocated to the (usually equidistant) discrete times $t_1, \ldots, t_N$ in $(0,T)$.
The goal is to determine the joint distribution function of
the extreme value vector $M_N = \big(M_{x,N}, M_{y,N}\big)$, where $M_{x,N} = \max\left\{ X_j\,; j = 1, \ldots, N \right\}$, and with a similar definition of $M_{y,N}$.
Specifically, the goal is to estimate $P(\xi, \eta) = \Prob\big(M_{x,N}\le \xi, M_{y,N} \le \eta \big)$ accurately for large values of $\xi$ and $\eta$. 

Following the pattern of Chapter~\ref{chap:acer}, the implementation of a sequence of approximations based on
conditioning is now outlined, where the first is a one-step memory approximation. This
approximation concept is described by \cite{naes85aor,naes90jsv}.

Hereafter, whenever expedient to ease the notation, $\zeta = (\xi,\,\eta)$ is used with a componentwise ordering relationship for $Z_i$, e.g. $Z_i\leq \zeta$ means $X_i \leq \xi$ and $Y_i \leq \eta$.
Also, the event $\mathcal{C}_{kj}(\zeta) = \mathcal{C}_{kj}(\xi, \eta) = \big\{ Z_{j-1} \leq \zeta, \ldots, Z_{j-k+1} \leq \zeta \big\}$  of $k-1$ consecutive componentwise non-exceedances $(k\geq2)$ is introduced. Then, from the definition of $P(\xi, \eta)$, it emerges that,
%%
\begin{align}\label{eq:aM2}
\begin{aligned}
& P(\xi, \eta) = P(\zeta) = \Prob \big( \mathcal{C}_{N+1,N+1}(\zeta) \big)  \\
& = \, \Prob \big( Z_N \leq \zeta \, |\,  \mathcal{C}_{NN}(\zeta) \big) \cdot \, \Prob \big(\mathcal{C}_{NN}(\zeta) \big)  \\
& = \prod_{j=2}^N \, \Prob \big( Z_j \leq \zeta \,|\, \mathcal{C}_{jj}(\zeta)\big) \cdot \, \Prob \big( \mathcal{C}_{22}(\zeta)\big).
\end{aligned}
\end{align}
%%
The first approximation of the sequence is obtained by assuming that the observed data pairs are independent, that is, data points $Z_i$ and $Z_j$
are statistically independent for all $i, j$, $i \ne j$, so that all conditioning in Eq.~\eqref{eq:aM2} can be neglected.

In this special case, it is obtained that,
%%
\begin{align}\label{eq:aM21}
\begin{aligned}
& P(\zeta) \approx {P}_1(\zeta) = \prod_{j=1}^N \, \Prob \big( Z_j \leq \zeta \big)  \\
& = \prod_{j=1}^N \, \big( 1 - \Prob(X_j > \xi) - \Prob(Y_j > \eta) + \Prob(X_j > \xi, Y_j > \eta) \big).
\end{aligned}
\end{align}
%%
%
Now, the designations $\alpha_{1j}(\xi; \eta)\,=\,\Prob(X_j > \xi)$, $\beta_{1j}(\eta; \xi)\,=\,\Prob(Y_j > \eta)$ and $\gamma_{1j}(\xi, \eta) \,=\,\gamma_{1j}(\zeta)\,=\,\Prob(Z_j > \zeta)$ for $1 \le j \le N$ are introduced. It should be noted that although neither $\alpha_{1j}(\xi; \eta)$ depends on $\eta$ nor $\beta_{1j}(\eta; \xi)$ 
depends on $\xi$, yet this notation is kept for the correct further derivations.

Eq.~\eqref{eq:aM21} can now be rewritten as,
\begin{align}\label{eq:aM22}
\begin{aligned}
P(\zeta) & \approx {P}_1(\zeta)  = \prod_{j=1}^N \big\{ 1 - \alpha_{1j}(\xi; \eta) - \beta_{1j}(\eta; \xi) + \gamma_{1j}(\zeta) \big\} \\
\approx & \exp \Big\{ - \sum_{j=1}^N \big( \alpha_{1j}(\xi; \eta) + \beta_{1j}(\eta; \xi) - \gamma_{1j}(\zeta) \big) \Big\} \,;
\end{aligned}
\end{align}
for large values of the components of $\zeta$, where the approximation $1-x \approx \exp(-x)$ for small values of $x$, has been applied, cf. Chapter~\ref{chap:acer}. 

In general, the variables $Z_{j}$ are statistically dependent in the componentwise sense.
In this case, clearly the first genuine conditioning approximation is obtained by neglecting all previous data except the immediate predecessor in Eq.~\eqref{eq:aM2}. While this approximation sometimes captures the effect of dependence in the time series on the extreme value distribution, it is in general not sufficient. Consequently, the following $k-1$-step memory approximation is adopted,
%%
\begin{align}\label{eq:aM23}
\begin{aligned}
P(\zeta) & \approx {P}_k(\zeta)  = \prod_{j=k}^N \, \Prob \big( Z_j \leq \zeta \,|\, \mathcal{C}_{kj}(\zeta) \big) \cdot \Prob \big( \mathcal{C}_{kk}(\zeta) \big) \\
& = \prod_{j=k}^N \,\Big\{ 1 - \Big( 1 - \Prob \big( Z_j \leq \zeta \,|\, \mathcal{C}_{kj}(\zeta) \big) \Big) \Big\} \cdot \Prob \big( \mathcal{C}_{kk}(\zeta) \big)\, \\
& \approx \exp \Big\{ - \sum_{j=k}^N \Big( 1 - \Prob \big( Z_j \leq \zeta \,|\, \mathcal{C}_{kj}(\zeta) \big) \Big) \Big\} \cdot \Prob \big( \mathcal{C}_{kk}(\zeta) \big),
\end{aligned}
\end{align}
%
for large values of the components of $\zeta$. Then, $k = 2$ means conditioning on only the previous observation.

By introducing the notations $\alpha_{kj}(\xi; \eta) = \Prob \big( X_j > \xi \,|\,\mathcal{C}_{kj}(\zeta) \big)$, $\beta_{kj}(\eta; \xi) = \Prob \big( Y_j > \eta \,|\,\mathcal{C}_{kj}(\zeta) \big)$ and 
$\gamma_{kj}(\xi, \eta) = \gamma_{kj}(\zeta) = \Prob \big( Z_j > \zeta\,|\,\mathcal{C}_{kj}(\zeta) \big)$,
for $k \le j \le N$, it can now be shown that, 
%%
\begin{align}\label{eq:aM36}
\exp \Big\{ - \sum_{j=k}^N \Big( 1 - \Prob \big( Z_j \leq \zeta \,|\, \mathcal{C}_{kj}(\zeta) \big) \Big) \Big\} = \exp \Big\{ - \sum_{j=k}^N \big( \alpha_{kj}(\xi; \eta) 
+ \beta_{kj}(\eta; \xi) - \gamma_{kj}(\zeta) \big) \Big\}\,. 
\end{align}

Similarly, it is found that, 
%%  
\begin{align}\label{eq:aM37}
\Prob \big(\mathcal{C}_{kk}(\zeta) \big) \approx \exp \Big\{ - \sum_{j=1}^{k-1} \big( \alpha_{jj}(\xi; \eta) + \beta_{jj}(\eta; \xi)
- \gamma_{jj}(\zeta) \big) \Big\}\,;
\end{align}
%%
for large values of the components of $\zeta$. Under these conditions, this leads to the result,
\begin{align}\label{eq:M38}
{P}_k(\zeta) \approx  &\exp \Big\{ - \sum_{j=k}^N \big( \alpha_{kj}(\xi; \eta) + \beta_{kj}(\eta; \xi) - \gamma_{kj}(\zeta) \big) \nonumber \\
&- \sum_{j=1}^{k-1} \big( \alpha_{jj}(\xi; \eta) + \beta_{jj}(\eta; \xi) - \gamma_{jj}(\zeta) \big) \Big\} \,. 
\end{align}


Thereby, based on the definition of the extreme value distribution $P(\zeta)$ and the properties of conditional probability, a set $\big\{P_k(\zeta)\big\}_{k=1}^{N}$ of conditional probability distributions has been constructed, which converges to the target 
distribution $P(\zeta)$ of the extreme value $M_N$ in the limit as $k$ increases for large values of the components of $\zeta$. 


For most applications, and for practical significance, the following assumption on this sequence of approximations is made: there is an effective $k_e$ satisfying $k_e \ll N$ such that $P(\zeta) = P_{k_e}(\zeta)$. Then, $P_{1}(\zeta) \le P_{2}(\zeta) \le \ldots \le P_{k_e}(\zeta) = P(\zeta)$. It may be noted that for a $k$-dependent stationary bivariate stochastic process $Z(t)$, that is, for data where $Z_i$ and $Z_j$ are independent componentwise whenever $|j-i| > k$, then $P(\zeta) = P_{k+1}(\zeta)$ exactly, as in the univariate case. 

It will be verified that the property $k_e \ll N$ is indeed satisfied for the wind speed data analysed in the present paper. Also, under this assumption, $\sum_{j=1}^{k-1}\big( \alpha_{jj}(\xi; \eta) + \beta_{jj}(\eta; \xi) - \gamma_{jj}(\zeta) \big)$ is generally negligible compared to $\sum_{j=k}^N \big( \alpha_{kj}(\xi; \eta) + \beta_{kj}(\eta; \xi) - \gamma_{kj}(\zeta) \big)$. This leads to the approximation ($k = k_e$), 
%%
\begin{align}\label{eq:M39}
P(\zeta) \approx P_k(\zeta) \approx \exp \Big\{ - \sum_{j=k}^N \big( \alpha_{kj}(\xi; \eta) 
+ \beta_{kj}(\eta; \xi) - \gamma_{kj}(\zeta) \big)\Big\}\,;
\end{align}
%%
for large values of $\xi$ and $\eta$, from which it emerges that for the estimation of the bivariate extreme value distribution, it is sufficient to estimate the sequence of functions $\big\{\alpha_{kj}(\xi; \eta) + \beta_{kj}(\eta; \xi) - \gamma_{kj}(\zeta)\big\}_{j=k}^{N}$.





\section{Appendix 2: Empirical Estimation of the Bivariate ACER Functions}\label{ApB}\index{bivariate ACER functions!empirical estimation}


To get a more compact representation, it is expedient to introduce the concept of $k$'th order bivariate average conditional exceedance rate (ACER) function as follows,
%%
\begin{align}\label{eq:bA0}
\mathcal{E}_k (\zeta)\,= \,\frac{1}{N-k+1}\,\sum\limits_{j=k}^N \big(\alpha_{kj}(\xi; \eta) + \beta_{kj}(\eta; \xi) - \gamma_{kj}(\zeta)\big)\: ; \: k=1,\,2,\ldots
\end{align}
%
Hence, when $N \gg k$, and for large values of the components of $\zeta$,
\begin{align}\label{eq:bA4}
P_k(\zeta) \approx \exp \left\{\, - \,(N-k+1)\, \mathcal{E}_k (\zeta) \right\} \,.
\end{align}
%%
Note that for higher dimensional cases, it is expedient to use an alternative formulation of the ACER functions, which is related to the problem of numerical calculation of the ACER functions, viz.,
\begin{align}\label{eq:bA1}
	\mathcal{E}_k (\zeta)\,= \,\frac{1}{N-k+1}\,\sum_{j=k}^N \Big( 1 - \Prob \big( Z_j \leq \zeta \,|\, \mathcal{C}_{kj}(\zeta) \big) \Big) \: ; \: k=\,2, 3, \ldots
\end{align}


A few more details on the numerical estimation of
the bivariate ACER functions are worth noting. It is useful to start by introducing a set of random functions. For $k = 2, \ldots, N$, and $k \le j \le N$, let,
%%
\begin{align}\label{eq:p0} % \hspace{-5mm}
\begin{aligned}
& A_{kj}(\xi; \eta)\  =\  \textbf{1} \Big\{{X_j > \xi}\,\cap\, \mathcal{C}_{kj}(\zeta)\Big\},  \\
& B_{kj}(\eta; \xi)\  =\  \textbf{1} \Big\{{Y_j > \eta}\,\cap\, \mathcal{C}_{kj}(\zeta)\Big\}, \\
& G_{kj}(\xi, \eta)\ =\ G_{kj}(\zeta) =\  \textbf{1} \Big\{{Z_j > \zeta}\,\cap\, \mathcal{C}_{kj}(\zeta) \Big\}, \\
& C_{kj}(\xi, \eta)\ =\ C_{kj}(\zeta) =\  \textbf{1} \Big\{\mathcal{C}_{kj}(\zeta)\Big\}, 
\end{aligned}
\end{align}
%%
where $\textbf{1}\{\mathcal{A}\}$ denotes the indicator function
of some event $\mathcal{A}$. 

From these definitions it follows that, for instance, 
%%
\begin{equation}\label{eq:bp4}
\alpha_{kj}(\xi; \eta) = \frac{\EE[A_{kj}(\xi; \eta)]}{\EE[C_{kj}(\zeta)]} \,,
\end{equation}
%%
where $\EE[\cdot]$ denotes the expectation operator. 
A similar equation holds for $\beta_{kj}(\eta; \xi)$ with $B_{kj}(\eta; \xi)$ instead of $A_{kj}(\xi; \eta)$ in the numerator, and for $\gamma_{kj}(\zeta)$ with $G_{kj}(\zeta)$ instead of $A_{kj}(\xi; \eta)$,  by analogy.

Assuming ergodicity of the process $Z(t) = \big(X(t), Y(t)\big)$, then obviously $\mathcal{E}_k (\zeta) =
\big(\alpha_{kk}(\xi; \eta) + \beta_{kk}(\eta; \xi) - \gamma_{kk}(\zeta)\big) = \ldots = \big(\alpha_{kN}(\xi; \eta) + \beta_{kN}(\eta; \xi) - \gamma_{kN}(\zeta)\big)$, and it may be assumed that for the bivariate time series at hand,
%
\begin{align}\label{eq:bp7}
\mathcal{E}_k (\zeta) 
= \lim_{N \to \infty}
\frac{\sum_{j=k}^N \big(\text{a}_{kj}(\xi; \eta) + \text{b}_{kj} (\eta; \xi) - \text{g}_{kj}(\zeta)\big)}{\sum_{j=k}^N {\text{c}_{kj}(\zeta)}} \,,
\end{align}
where $\text{a}_{kj}(\xi; \eta)$, $\text{b}_{kj} (\eta; \xi)$, $\text{g}_{kj}(\zeta)$ and $\text{c}_{kj}(\zeta)$ are the realized values of $A_{kj}(\xi; \eta)$, $B_{kj}(\eta; \xi)$, $G_{kj}(\zeta)$ and $C_{kj}(\zeta)$, respectively, for the observed time series.
%

While Eq.~\eqref{eq:bp7} clearly displays how the conditional exceedances are counted for the estimation of the bivariate ACER functions, this approach becomes very cumbersome for higher dimensional problems. Fortunately, this can be circumvented by invoking Eq.~\eqref{eq:bA1}. It is observed that,
%
\begin{equation}\label{eq:bP4}
	\Prob \big( Z_j \leq \zeta \,|\, \mathcal{C}_{kj}(\zeta) \big) = \frac{\Prob \big( Z_j \leq \zeta \,\cap\, \mathcal{C}_{kj}(\zeta) \big)}{\Prob \big( \mathcal{C}_{kj}(\zeta) \big)} = \frac{\EE[C_{k+1,j+1}(\zeta)]}{\EE[C_{kj}(\zeta)]}  \,,
\end{equation}
%%
since clearly $\{Z_j \leq \zeta \} \,\cap\, \mathcal{C}_{kj}(\zeta) = \mathcal{C}_{k+1,j+1}(\zeta)$. A consequence of this, is that Eq.~\eqref{eq:bp7} may be replaced by the following equation,
%
\begin{align}\label{eq:bP7}
	\mathcal{E}_k (\zeta) 
	= 1 - \lim_{N \to \infty}
	\frac{\sum_{j=k}^{N-1} \text{c}_{k+1,j+1}(\zeta)}{\sum_{j=k}^{N} {\text{c}_{kj}(\zeta)}} \,,
\end{align}
which is a form particularly suitable for estimating $\mathcal{E}_k (\zeta)$ for high dimensional distributions.
%




Clearly, $\lim\limits_{\xi,\eta \to \infty} \EE[C_{kj}(\zeta)] = 1$. Hence, 
%%
\begin{equation}\label{eq:bp7_1}
\lim\limits_{\xi,\eta \to \infty} {\tilde{\mathcal{E}}_k(\zeta)}\,/\,{\mathcal{E}_k(\zeta)} \,= 1\,,
\end{equation}
%%
where
%%
\begin{align}\label{eq:bp8}
\tilde{\mathcal{E}}_k(\zeta) = \lim_{N \to \infty} \frac{\sum_{j=k}^N \Big( \EE[A_{kj}(\xi; \eta)] + \EE[B_{kj}(\xi; \eta)] - \EE[G_{kj}(\zeta)] \Big)}{ N - k +1 } \,. % (\xi; \eta) (\eta; \xi) (\xi, \eta)
\end{align}
%%
The advantage of using the modified bivariate ACER function
$\tilde{\mathcal{E}}_k(\zeta)$ for $k \ge 2$ is that it is somewhat easier
to use for non-stationary or long-term statistics than
$\mathcal{E}_k(\zeta)$. This aspect is discussed in Chapter~\ref{chap:acer}, see also \cite{kana13jweia}. Since our focus is on the values of the ACER functions at the extreme levels, any function may be used that
provides correct predictions of the appropriate ACER function at these extreme
levels.
%

Now, let us look at the problem of estimating confidence
intervals for the bivariate ACER function. 
If several realizations of the time series  $Z(t) = \big(X(t), Y(t)\big)$ are provided or the time series can be appropriately sectioned into 
several records, e.g. several annual, or other time span records, the sample
estimate of $\mathcal{E}_k(\zeta)$ would be,
%%
\begin{equation}\label{eq:o5}
\hat{\mathcal{E}}_k(\zeta) =\frac{1}{R}\sum_{r=1}^R
\hat{\mathcal{E}}_k^{(r)}(\zeta)\,,
\end{equation}
%%
where $R$ is the number of realizations (samples). $\hat{\mathcal{E}}_k^{(r)}(\zeta)$ is estimated using the result from Eq.~\eqref{eq:bp7}
for the stationary time series, or using the result from Eq.~\eqref{eq:bp8}
for non-stationary time series, where the index $(r)$ refers to realization no. $r$.
The sample standard deviation $\hat{s}_{k}(\zeta)$ then can be estimated by the standard
formula,
%%
\begin{equation}\label{eq:bp12}
\hat{s}_{k}(\zeta)^2 = \frac{1}{R-1} \sum_{r=1}^R  \, \left(\hat{\mathcal{E}}_k^{(r)}(\zeta) - \hat{\mathcal{E}}_k(\zeta)\right)^2\,.
\end{equation}
%%
Assuming that realizations are independent, Eq.~(\ref{eq:bp12}) leads to a good  
approximation of the 95\% confidence interval CI =
$\big(\CI^-(\zeta), \CI^+(\zeta)\big)$ for the value
$\mathcal{E}_k(\zeta)$, where,
%%
\begin{equation}\label{eq:bcb}
\CI^\pm(\zeta)=\hat{\mathcal{E}}_k(\zeta) \pm \tau \cdot \frac{\hat{s}_{k}(\zeta)}{\sqrt R } \,,
\end{equation}
%%
and $\tau = t^{-1}\Big((1-0.95)/2,\ R - 1\Big)$  is the corresponding quantile of the Student's $t$-distribution with $R - 1$ degrees of freedom.





