

\section{Introduction}

Extreme value statistics, even in applications, is generally
based on asymptotic results. This is done either by assuming
that the epochal extremes, for example yearly extreme wind speeds at a given location,
are distributed according to the so called generalized (asymptotic) extreme
value distribution with unknown parameters to be estimated on the
basis of the observed data, cf. Chapter~\ref{chap:aevt}, or \cite{cole01,begosete04}. Or, it is assumed that the exceedances above high thresholds follow a generalized (asymptotic) Pareto
distribution with parameters that are estimated from the data, cf. Chapter~\ref{chap:pot}, or \cite{cole01,begosete04,dasm90,reth07}. 
As was discussed in Chapter~\ref{chap:challenges}, the major problem with both of these approaches is
that the asymptotic extreme value theory itself cannot be used in practice to decide to what extent it is applicable for the
observed data. And since statistical tests to decide this issue are rarely precise enough to settle this problem, 
the assumption that a specific asymptotic extreme value distribution is the appropriate distribution for the
observed data, is based more or less on faith or convenience. 

On the other hand, one can reasonably assume that in many cases long time series obtained from practical measurements do contain values that are large enough to provide useful information about extreme events that are truly asymptotic. This cannot be strictly proved in general, of course, but the accumulated experience indicates that asymptotic extreme value distributions do provide reasonable, if not always very accurate, predictions when based on measured data. This is amply documented in the vast literature on the subject, and good references to this literature are \cite{begosete04,emklmi97,fahure04}. 

However, even if the situation might be tolerable, it is clearly not satisfactory. In an effort to improve on the described state of affairs, an approach to the extreme value prediction problem has been developed that is less restrictive and more flexible than the ones using only asymptotic theory \citep{naga09ss,nagaka13jps}. The approach is based on two separate components which are designed to improve on two important aspects of extreme value prediction based on observed data. 

\begin{list}{$\Rightarrow$}{}
\item The first component provides a nonparametric empirical representation of the extreme value distribution inherent in the data. It also has the capability to accurately capture and display the effect of statistical dependence in the data on the extreme value distribution, which opens for the opportunity of using all the available data in the analysis. 
\item The second component is then constructed based on the realization that all, or at least almost all, the available data are sub-asymptotic. Hence, it is important to set up the method of analysis to also incorporate sub-asymptotic data into the estimation of extreme value distributions, which will be shown to have some importance for accurate prediction. 
\end{list}

The proposed method has been used on a wide variety of estimation problems, and the experience is that it represents a very powerful addition to the toolbox of methods for extreme value estimation and prediction. Several such examples will be presented in this book, cf. Chapters~\ref{chap:paeva}, \ref{chap:mc} and \ref{chap:ewl}. A recent and interesting example of its practical use on the distribution of defects in metallic materials is given by \cite{cena12ijf}. Fatigue failures in such materials is intimately related to the distribution of defects, cf. \cite{mura02}. Another recent example, is its application to sea ice dynamics proposed by \cite{sinalech22coe}. Needless to say, what is presented in this chapter is by no means considered the end of the development of the ACER method. It is a novel method, and it is to be expected that several aspects of the proposed approach will see improvements.  





\section{A Sequence of Conditioning Approximations}\label{M}

In this section a sequence of nonparametric distribution functions will be constructed that converges to the exact extreme value distribution for the time series considered. This constitutes the core of the proposed approach. 

Consider a stochastic process $Z(t)$, which has been observed over
a time interval, $(0,T)$ say. Assume that values $X_1, \ldots,
X_N$, which have been derived from the observed process, are
allocated to the discrete times $t_1, \ldots, t_N$ in $(0,T)$.
This could be simply the observed values of $Z(t)$ at each $t_j$,
$j = 1, \ldots, N$, or it could be average values or peak values
over smaller time intervals centered at the $t_j$'s. Our goal is to accurately determine the distribution function of the extreme value $M_N = \max\{ X_j\,; j = 1, \ldots, N \}$.
Specifically, we want to estimate $P(\eta) = \Prob(M_N \le \eta)$
accurately for large values of $\eta$ without asymptotics. Clearly, $P(\eta) = \Prob(X_1 \le \eta, X_2 \le \eta, \ldots, X_N \le \eta)$. Since $N$ typically will be a large number, direct estimation of this joint distribution function is not a practical option. Hence, we need to develop another approach. 

An underlying premise for the development in this chapter is that a rational approach to the
study of the extreme values of the sampled time series is to
consider exceedances of the individual random variables $X_j$
above given thresholds, as in classical extreme value theory. The
alternative approach of considering the exceedances by upcrossing
of given thresholds by a continuous stochastic process has already been discussed in Chapter~\ref{chap:cr}, see also \cite{naga08jem, nagaha07oe}. 
The approach taken in the present chapter seems to be the appropriate way to deal with the
recorded data time series of, for example, the hourly or daily largest wind
speeds observed at a given location, just to cite a concrete example. 


The following basic rule from probability theory: $\Prob( \mathcal{A} \cap \mathcal{B}) = \Prob( \mathcal{A}| \mathcal{B})\Prob( \mathcal{B})$ for two events $\mathcal{A}$ and $ \mathcal{B}$, turns out to be an important key to estimating $P(\eta)$. Using this basic rule repeatedly, it is obtained that,
\begin{align}\label{eq:ca1}
&P(\eta) = \Prob(M_N \le \eta) = \Prob\{X_N \leq \eta,\ldots, X_1 \leq \eta\} \nonumber \\
&= \Prob\{X_N \leq \eta |\, X_{N-1} \leq \eta,\ldots, X_1 \leq \eta\}  \Prob\{X_{N-1} \leq \eta, \ldots, X_1 \leq \eta \} \nonumber \\
&= \Prob\{X_N \leq \eta |\, X_{N-1} \leq \eta,\ldots, X_1 \leq \eta\}  \nonumber \\
&\:\: \cdot \:\: \Prob\{X_{N-1} \leq \eta | \,X_{N-2} \leq \eta, \ldots, X_1 \leq \eta \}  \Prob\{X_{N-2} \leq \eta, \ldots, X_1 \leq \eta \} \nonumber \\   
&\vdots \nonumber  \\
&= \prod_{j=2}^N \Prob\{X_j \leq \eta |\, X_{j-1} \leq \eta,\ldots,
X_1 \le \eta\} \cdot \Prob(X_1 \leq \eta).
\end{align}

In general, the variables $X_j$ are statistically dependent.
Hence, instead of assuming that all the $X_j$ are statistically
independent, which leads to the classical approximation,
\begin{equation}\label{eq:ca2}
P(\eta) \approx P_1(\eta) := \prod_{j=1}^N \Prob(X_j \leq \eta) \,,
\end{equation}
where $:=$ means 'by definition', the following one-step memory approximation will to
some extent account for the dependence between the $X_j$'s,
\begin{equation}\label{eq:ca3}
\Prob\{X_j \leq \eta |\, X_{j-1} \leq \eta, \ldots, X_1 \leq \eta \}
\approx \Prob\{ X_j \leq \eta |\, X_{j-1} \leq \eta\},
\end{equation}
for $2 \le j \le N$. With this approximation, it is obtained that,
\begin{equation}\label{eq:ca4}
P(\eta) \approx  P_2(\eta) := \prod_{j=2}^N \Prob\{ X_j \leq \eta |\, X_{j-1} \leq \eta\} \Prob(X_1 \leq \eta)\,.
\end{equation}
By conditioning on one more data point, the one-step memory approximation is extended to,
\begin{equation}\label{eq:ca5}
\Prob\{X_j \leq \eta |\, X_{j-1} \leq \eta,\ldots, X_1 \le \eta\}
\approx \Prob\{X_j \leq \eta |\, X_{j-1} \leq \eta, X_{j-2} \leq
\eta \} \,,
\end{equation}
where $3 \le j \le N$, which leads to the approximation, 
\begin{align}\label{eq:ca6}
P(\eta) \approx P_3(\eta) :=  \prod_{j=3}^N \, &\Prob\{X_j \leq \eta |\, X_{j-1} \leq \eta, X_{j-2} \leq
\eta \}\,\nonumber \\ \cdot \, &\Prob\{ X_2 \leq \eta |\, X_1 \leq \eta\} \,\Prob(X_1 \leq \eta)\,.
\end{align}
For a general $k$, $2 \le k \le N$, it is obtained that,
\begin{align}\label{eq:ca7}
P(\eta) \approx P_k(\eta) :=  \prod_{j=k}^N \, &\Prob\{X_j \leq \eta |\, X_{j-1} \leq \eta, \ldots, X_{j-k+1} \leq
\eta \} \,\nonumber \\ 
\cdot \, \prod_{j=2}^{k-1} &\Prob\{ X_j \leq \eta |\, X_{j-1} \leq \eta\, \ldots, X_{1} \leq \eta \} \cdot \,\Prob(X_1 \leq \eta)\,,
\end{align}
where $P(\eta) = P_N(\eta)$. It follows that the sequence of approximations $P_1(\eta), P_2(\eta), \ldots$ constitute a sequence of increasingly accurate representations of the exact extreme value distribution $P(\eta)$.


It should be noted that the one-step memory approximation adopted above is not a Markov chain approximation, as being discussed in \cite{smit92,cole94,smtaco97}, nor do the $k$-step memory approximations lead to $k$th-order Markov chains, which are proposed in \cite{yun98,yun00}. An effort to relinquish the Markov chain assumption to obtain an approximate distribution of clusters of extremes is reported by \cite{sege05}.   


It is now necessary to have a closer look at the values for $P(\eta)$ obtained by
using Eq.~\eqref{eq:ca7} as compared to Eq.~\eqref{eq:ca2}. Eq.~\eqref{eq:ca2} can be rewritten in the form,
\begin{equation}\label{eq:ca8}
P(\eta) \approx P_1(\eta) = \prod_{j=1}^N \big(1 - \alpha_{1j}(\eta) \big),
\end{equation}
where $\alpha_{1j}(\eta) = \mbox{Prob}\{X_j > \eta\}$, $j = 1, \ldots, N$. Then the approximation based on assuming independent data can be written as,
\begin{equation}\label{eq:ca9}
P(\eta) \approx F_1(\eta) := \exp \Big(- \sum_{j=1}^N
\alpha_{1j}(\eta) \Big)\,, \:\: \eta \to \infty \,,
\end{equation}
noting that $1-x \approx \exp(-x)$ with high accuracy for small $x$. The relative error of this approximation is less than 0.5\% for values $|x|<0.1$ and it decreases rapidly for decreasing values of $|x|$. In this chapter, the notation $\eta \to \infty$ is adopted to denote large values of $\eta$ relative to the typical values. Hence, it should not be construed as a true limiting operation.

Similarly, Eq.~\eqref{eq:ca7} can be expressed as,
\begin{equation}\label{eq:ca10}
P(\eta) \approx P_k(\eta) = \prod_{j=k}^N \big( 1 - \alpha_{kj}(\eta)
\big)\, \prod_{j=1}^{k-1} \big( 1 - \alpha_{jj}(\eta) \big)  \,,
\end{equation}
where $\alpha_{kj}(\eta) = \Prob\{X_j> \eta \,|\, X_{j-1} \leq \eta, \ldots, X_{j-k+1} \leq \eta \}$, for
$j \ge k \ge 2$, denotes the exceedance probability conditional on $k-1$ previous
non-exceedances. From Eq.~\eqref{eq:ca10} it is obtained that,
\begin{equation}\label{eq:ca11}
P(\eta) \approx F_k(\eta) := \exp\Big( - \sum_{j=k}^N
\alpha_{kj}(\eta) - \sum_{j=1}^{k-1}\alpha_{jj}(\eta) \Big)\,, \:\: \eta \to \infty \,,
\end{equation}
and $F_k(\eta) \to P(\eta)$ as $k \to N$ with $F_N(\eta) = P(\eta)$ for $\eta \to \infty$. 

For the sequence of approximations $F_k(\eta)$ to have practical significance, it is implicitly assumed that there is a cut-off value $k_c$ satisfying $k_c \ll N$ such that effectively $F_{k_c}(\eta) = F_N(\eta)$. It may be noted that for $k$-dependent stationary data sequences, that is, for data where $X_i$ and $X_j$ are independent whenever $|j-i| > k$, then $P(\eta) = P_{k+1}(\eta)$ exactly, and, under rather mild conditions on the joint distributions of the data, $\lim_{N \to \infty} P_{1}(\eta) = \lim_{N \to \infty} P(\eta)$ \citep{watson54}. In fact, it can be shown that $\lim_{N \to \infty} P_{1}(\eta) = \lim_{N \to \infty} P(\eta)$ is true for weaker conditions than $k$-dependence \citep{leliro83}. However, for finite values of $N$, the picture is much more complex, and purely asymptotic results should be used with some caution. 


Returning to Eq.~\eqref{eq:ca11}, extreme value prediction by the conditioning
approach described above reduces to estimation of (combinations)
of the $\alpha_{kj}(\eta)$ functions. In accordance with the previous assumption about a cut-off value $k_c$, for all $k$-values of interest, $k \ll N$, so that $\sum_{j=1}^{k-1}\alpha_{jj}(\eta)$ is effectively negligible compared to $\sum_{j=k}^N
\alpha_{kj}(\eta)$. Hence, for simplicity, the following approximation is adopted, which is applicable to both stationary and nonstationary data,
\begin{equation}\label{eq:ca12}
F_k(\eta) = \exp \Big( - \sum_{j=k}^N \alpha_{kj}(\eta)  \Big)\, \: , \:\: k \ge 1\,.
\end{equation}



Going back to the definition of $\alpha_{1j}(\eta)$, it follows that $\sum_{j=1}^N
\alpha_{1j}(\eta)$ is equal to the expected number of exceedances
of the threshold $\eta$ during the time interval $(0,T)$.
Eq.~(\ref{eq:ca9}) therefore expresses the approximation that the
stream of exceedance events constitute a (non-stationary) Poisson
process. This opens for an understanding of Eq.~(\ref{eq:ca12}) by interpreting the expressions
$\sum_{j=k}^N \alpha_{kj}(\eta)$ as the
expected effective number of (assumed) independent exceedance events provided by
conditioning on $k-1$ previous observations.





\section{Empirical Estimation of the Average Conditional Exceedance Rates}\label{Emp}


The concept of average conditional exceedance rate\index{average conditional exceedance rate} (ACER) of order $k$ is now introduced as follows,
\begin{align}\label{eq:ee1}
{\varepsilon}_k(\eta) = \frac{1}{N-k+1} \sum_{j=k}^N \,
\alpha_{kj}(\eta)\,, \: k = 1, 2, \ldots
\end{align}
It is noted that this empirical ACER function\index{ACER function} also depends on the number of data points $N$. In contrast to the average upcrossing rate of Chapter~\ref{chap:cr}, which would typically express the average number of upcrossings per time unit, the ACER functions are exceedance rates per data point. 

The behaviour and diagnostic power of the ACER functions will be demonstrated for several examples in Chapter~\ref{chap:paeva}. In terms of the ACER function, we may now write,
\begin{equation}\label{eq:ca13}
F_k(\eta) = \exp \Big( - (N-k+1) {\varepsilon}_k(\eta) \Big)\, \: , \:\: k \ge 1\,.
\end{equation}

In practice, there are typically two scenarios for the underlying
process $Z(t)$. Either we may consider it to be a stationary
process, or, in fact, even an ergodic process, which allows the replacement of ensemble averages with time averages, cf. \cite{doob53,crle67,woha85}. The alternative is
to view $Z(t)$ as a process that depends on certain parameters
whose variation in time may be modelled as an ergodic process in
its own right. For each set of values of the parameters, the
premise is that $Z(t)$ can then be modelled as an ergodic process. This
would be the scenario that can be used to model long-term
statistics \citep{naes84baor,scfara91}.

For both these scenarios, the empirical estimation of the 
ACER function\index{ACER function!empirical estimation} ${\varepsilon}_k(\eta)$ proceeds in a
completely analogous way by counting the total number of
favourable incidents, that is, exceedances combined with the
requisite number of preceding non-exceedances, for the total data
time series and then finally dividing by $N-k+1 \approx N$. This
can be shown to apply for the long-term situation, as briefly discussed below.


A few more details on the numerical estimation of
${\varepsilon}_k(\eta)$ for $k \ge 2$ may be appropriate. Initially,
the following random functions are introduced,
\begin{equation}\label{eq:ee2}
A_{kj}(\eta) = \mathbf{1} \{X_j > \eta, X_{j-1} \leq \eta, \ldots,
X_{j-k+1} \leq \eta \}\,, \: j = k, \ldots, N, \: k = 2, 3, \ldots
\end{equation}
and
\begin{equation}\label{eq:ee3}
B_{kj}(\eta) = \mathbf{1} \{X_{j-1} \leq \eta, \ldots, X_{j-k+1}
\leq \eta \}\,, \: j = k, \ldots, N, \: k = 2, \ldots\,,
\end{equation}
where $\mathbf{1} \{ \mathcal{A} \}$ denotes the indicator function
of some event $\mathcal{A}$, that is, $\mathbf{1} \{ \mathcal{A} \} = 1$ if the event occurs, $\mathbf{1} \{ \mathcal{A} \}$ = 0 if not. Then, since $\EE[\mathbf{1} \{ \mathcal{A} \}] = \Prob\{\mathcal{A}\}$, where $\EE[\cdot]$ denotes the expectation operator, it follows that,
\begin{equation}\label{eq:ee4}
\alpha_{kj}(\eta) = \frac{\Prob\{X_j> \eta , X_{j-1} \leq \eta, \ldots, X_{j-k+1} \leq \eta \}}{\Prob\{ X_{j-1} \leq \eta, \ldots, X_{j-k+1} \leq \eta \}} = \frac{\EE[A_{kj}(\eta)]}{\EE[B_{kj}(\eta)]} \,,
\: j = k, \ldots, N, \: k = 2, \ldots\,,
\end{equation}
Assuming an ergodic process, which implies stationarity, then obviously ${\varepsilon}_k(\eta) =
\alpha_{kk}(\eta) = \ldots = \alpha_{kN}(\eta)$, and by ergodicity, replacing ensemble means with corresponding time averages, it may be assumed that for the time series at hand,
\begin{equation}\label{eq:ee5}
{\varepsilon}_k(\eta) = \lim_{N \to \infty}
\frac{\sum_{j=k}^N a_{kj}(\eta)}{\sum_{j=k}^N b_{kj}(\eta)} \,,
\end{equation}
where $a_{kj}(\eta)$ and $b_{kj}(\eta)$ are the realized values of $A_{kj}(\eta)$ and $B_{kj}(\eta)$, respectively, for the observed time series. Denoting a realization of the time series $X_1, X_2, \ldots$ by $x_1, x_2, \ldots$, then $a_{kj}(\eta) = 1$ if $x_j> \eta , x_{j-1} \leq \eta, \ldots, x_{j-k+1} \leq \eta$, otherwise $a_{kj}(\eta) = 0$. Similarly, $b_{kj}(\eta) = 1$ if $x_{j-1} \leq \eta, \ldots, x_{j-k+1} \leq \eta$, otherwise $b_{kj}(\eta) = 0$.


Clearly, $\lim_{\eta \to \infty} \EE[B_{kj}(\eta)] = 1$. Hence, $\lim_{\eta \to \infty}
\tilde{\varepsilon}_k(\eta) / {\varepsilon}_k(\eta)$ $= 1$,
where,
\begin{equation}\label{eq:ee6}
\tilde{\varepsilon}_k(\eta) = \frac{\sum_{j=k}^N \EE[A_{kj}(\eta)]}{ N -k +1 } \,.
\end{equation}
The advantage of using the modified ACER function\index{ACER function!modified}
$\tilde{\varepsilon}_k(\eta)$ for $k \ge 2$, is that it is easier
to use for non-stationary or long-term statistics than
${\varepsilon}_k(\eta)$. Since our focus is on the values
of the ACER functions at the extreme levels, we may use any function that
provides correct predictions of the appropriate ACER function at these extreme
levels.

To see why Eq.~\eqref{eq:ee6} may be applicable for nonstationary time series\index{nonstationary time series}, it is recognized that,
\begin{align}\label{eq:ee7}
P(\eta) &\approx  \exp \big( - \sum_{j=k}^N \alpha_{kj}(\eta)
\big) = \exp \Big( - \sum_{j=k}^N  \frac{\EE[A_{kj}(\eta)]}{\EE[B_{kj}(\eta)]} \Big) \nonumber \\
&\underset{\eta \to \infty}{\simeq} \exp \big( - \sum_{j=k}^N  \EE[ A_{kj}(\eta)] \big) \,. 
\end{align}
If the time series can be segmented into $K$ blocks such that $\EE[A_{kj}(\eta)]$ remains approximately constant within each block and such that $\sum_{j \in C_i}  \EE[ A_{kj}(\eta)] \approx \sum_{j \in C_i} a_{kj}(\eta)$ for a sufficient range of $\eta$-values, where $C_i$ denotes the set of indices for block no. $i$, $i = 1, \ldots, K$, then $\sum_{j=k}^N  \EE[ A_{kj}(\eta)] \approx \sum_{j=k}^N a_{kj}(\eta)$. Hence,
\begin{equation}\label{eq:ee8}
P(\eta) \approx  \exp \big( - (N - k + 1) \hat{\varepsilon}_k(\eta) \big)\,, \,\, \eta \to \infty,
\end{equation}
where, in this case, 
\begin{equation}\label{eq:ee9}
\hat{\varepsilon}_k(\eta) =  \frac{1}{N - k + 1} \sum_{j=k}^N a_{kj}(\eta) \,,
\end{equation}
which is the empirical counterpart of Eq.~\eqref{eq:ee6}.
   

It is of interest to note what events are actually counted for the
estimation of the various ${\varepsilon}_k(\eta)$, $k \ge 2$.
Let us start with ${\varepsilon}_2(\eta)$. It follows from the
definition of ${\varepsilon}_2(\eta)$ that
${\varepsilon}_2(\eta)\,(N-1)$ can be interpreted as the
expected number of exceedances above the level $\eta$ satisfying
the condition that an exceedance is counted only if it is
immediately preceded by a non-exceedance. Hence, if a realization of a stochastic process is sampled with a sufficiently high frequency to approximate the realization with good accuracy, then ${\varepsilon}_2(\eta)\,(N-1)$ is approximately the expected number of upcrossings of $\eta$ among the $N-1$ data points. This implies that ${\varepsilon}_2(\eta)$ is approximately the average upcrossing rate per data point. Therefore, if the $N-1$ data points correspond to a time interval of length $T$, then ${\varepsilon}_2(\eta)\,(N-1)/T$ expresses the average upcrossing rate per time unit.  

Assuming that we are studying a narrow band stochastic response process, sampling a realization by extracting the local peak values would be a typical procedure for studying the extremes. A reinterpretation of ${\varepsilon}_2(\eta)\,(N-1)$ for this case would lead us to the conclusion that it equals the average
number of clumps of exceedances above $\eta$ for the realizations
considered, where a clump of exceedances is defined as a maximum
number of consecutive exceedances of the peak values above $\eta$. Figure~\ref{fig:acer3.1} illustrates how clumps of exceedances would be counted for the estimation of ${\varepsilon}_2(\eta)$. A scrutiny of the figure shows that four clumps can be identified.  


\begin{figure}[!hbt]
	\centerline{\includegraphics[width=4in]{Figures/ts_exceedances.pdf}}
	\caption{An illustration of clumps of exceedances counted for the estimation of ${\varepsilon}_2(\eta)$. }
	\label{fig:acer3.1}
\end{figure}



In general, ${\varepsilon}_k(\eta)\,(N-k+1)$ then equals the average number
of clumps of exceedances\index{clumps of exceedances} above $\eta$ separated by at least $k-1$
non-exceedances. Hence, for the illustrating example of Figure~\ref{fig:acer3.1}, it follows that for the estimation of ${\varepsilon}_3(\eta)$, the two arrows show that what were three clumps for the estimation of ${\varepsilon}_2(\eta)$, now become one clump for the estimation of ${\varepsilon}_3(\eta)$. Hence, in Figure~\ref{fig:acer3.1}, we can now identify only two clumps for the estimation of ${\varepsilon}_3(\eta)$.  

If the time series analysed is obtained by extracting local peak values from a narrow band response process, it is interesting to note that there is a certain similarity between the ACER approximations and the envelope approximations for extreme value prediction \citep{naga08jem,vanm75}. For alternative statistical approaches to account for the effect of clustering on the extreme value distribution, the reader may consult \cite{lead83,hsin87,hsin91,lead95,fese03,robe09}. In these works the emphasis is on the notion of the extremal index, which characterizes the clumping or clustering tendency of the data and its effect on the extreme value distribution, cf. Section~\ref{sec:dss}. In the ACER functions, these effects are automatically accounted for. 


Now, let us look at the problem of estimating a confidence
interval\index{ACER function!confidence interval} for ${\varepsilon}_k(\eta)$, assuming a stationary time series. If $R$ realizations of the requisite length of the time series is available, or, if one long realization can be segmented into $R$ subseries, then the sample standard
deviation $\hat{s}_k(\eta)$ can be estimated by the standard formula,
%
\begin{equation}\label{eq:ee10}
\hat{s}_k(\eta)^2 = \frac{1}{R-1} \sum_{r=1}^R  \, \Big(
         \hat{\varepsilon}_k^{(r)}(\eta) - \hat{\varepsilon}_k(\eta) \,
         \Big)^2\,,
\end{equation}
%
where $\hat{\varepsilon}_k^{(r)}(\eta)$ denotes the ACER function estimate from realization no. $r$, and $\hat{\varepsilon}_k(\eta) = \sum_{r=1}^R      \hat{\varepsilon}_k^{(r)}(\eta) / R$.

Assuming that realizations are independent, for a suitable number
$R$, e.g. $R \ge 20$, Eq.~(\ref{eq:ee10}) leads to a good
approximation of the 95 \% confidence interval CI =
$\big(C^-(\eta), C^+(\eta)\big)$ for the value
${\varepsilon}_k(\eta)$, where,
%
\begin{equation}\label{eq:ee11}
C^\pm(\eta)=\hat{\varepsilon}_k(\eta) \pm 1.96 \,
\hat{s}_k(\eta) / \sqrt{R}\,.
\end{equation}


Alternatively, and which also applies to the non-stationary case, it is consistent with the adopted approach to assume that the stream of conditional exceedances over a threshold $\eta$ constitute a Poisson process, possibly non-homogeneous. Hence, the variance of the estimator $\hat{E}_k(\eta)$ of $\tilde{\varepsilon}_k(\eta)$, where,
%
\begin{equation}\label{eq:ee12}
\hat{E}_k(\eta) = \frac{\sum_{j=k}^N A_{kj}(\eta)}{ N -k +1 } \,,
\end{equation}
%
is $\Var[\hat{E}_k(\eta)] = \tilde{\varepsilon}_k(\eta)$. Therefore, for high levels $\eta$, the approximate limits of a 95~\% confidence interval of ${\varepsilon}_k(\eta)$, can be written as,
%
\begin{equation}\label{eq:ee13}
C^\pm(\eta)=\hat{\varepsilon}_k(\eta) \Big( 1 \pm \frac{1.96 }{\sqrt{(N-k+1)\hat{\varepsilon}_k(\eta)}} \Big)\,.
\end{equation}
  



\section{Long-Term Extreme Value Analysis by the ACER Method}\label{sec:ltevacer}\index{ACER method!long-term analysis}


In Section~\ref{sec:ltev} we studied long-term extreme value distributions for the point process method. This was related to a scatter diagram of short term environmental conditions. If the whole time series over a long-term scenario is available, we have shown in the previous section that the long-term statistics using ACER functions may be estimated directly from the time series, cf. Eq.~\ref{eq:ee8}. However, in many cases it would be more practical to analyze each short term condition separately and combine the obtained ACER functions after that. This would, e.g. be the typical approach in a simulation based long-term statistical analysis where the short term response time series would be simulated and the resulting time series subjected to an ACER analysis. The obtained ACER functions would then be combined to form the long-term ACER function. One advantage of doing a long-term analysis this way, would be the opportunity to apply the ACER function definition of Eq.~\eqref{eq:ee5} to each short term condition. 

Adapting Eq.~\eqref{lev1} to the discrete scenario for the ACER analysis over the scatter diagram of Table~\ref{tab:10.2}, an estimated long-term ACER function would be expressed as follows,
\begin{equation}\label{eq:lev1}
\hat{\varepsilon}_k(\eta) =  \sum_{h_i} \sum_{t_j}
\hat{\varepsilon}_k(\eta|h_i, t_j)\, f_{H_sT_s}(h_i, t_j) \, \Delta h_i \Delta t_j \,,
\end{equation}
where the ACER functions $\hat{\varepsilon}_k(\eta|h_i, t_j)$ have been estimated for each separate short term condition $(h_i, t_j)$, or $(i,j)$ for short.

An alternative, equivalent formulation is obtained as in Subsection~\ref{ssec:ltev}. Assume that the number of sea states in condition $(i,j)$ is $K_{ij}$, $i = 1, \ldots, m$ and $j = 1, \ldots, n$, and $K = \sum_{i=1}^m \sum_{j=1}^n K_{ij}$, cf. Subsection~\ref{ssec:apv}. We may then write,

\begin{equation}\label{eq:lev2}
\hat{\varepsilon}_k(\eta) =  \sum_{i=1}^{m} \sum_{j=1}^{n}\hat{\varepsilon}_k^{(ij)}(\eta)\,\: \frac{K_{ij}}{K}\,,
\end{equation}
where the ACER function $\hat{\varepsilon}_k^{(ij)}(\eta)$ is estimated for condition $(i,j)$. 
So, again we obtain the long-term extreme value distribution as,

\begin{equation}\label{eq:lev3}
P(\eta) \approx  \exp \big( - (N - k + 1) \hat{\varepsilon}_k(\eta) \big)\,,
\end{equation}
using either Eq.\eqref{eq:lev1} or \eqref{eq:lev2}, where $N$ is the number of data over the long-term period. However, please note the cautionary comments on the use of a scatter diagram for long-term analysis in Subsection~\ref{ssec:ltev}.




\section{Estimation of Extremes for the Asymptotic Gumbel Case}\label{P}\index{ACER function!asymptotic Gumbel case}


The second component of the approach to extreme value estimation presented in this chapter, was originally derived for a time series with an asymptotic extreme value distribution of the Gumbel type, cf. \cite{naga09ss}. This case is therefore presented first, also because the extension of the asymptotic distribution to a parametric class of extreme value distribution tails that are capable of capturing to some extent subasymptotic behaviour is more transparent, and perhaps more obvious, for the Gumbel case. The reason behind the efforts to extend the extreme value distributions to the subasymptotic range is the fact that the ACER functions allow us to use not only asymptotic data, which is clearly an advantage since proving that observed extremes are asymptotic is not possible.      


The effect of the asymptotic distribution being of the Gumbel type on the possible subasymptotic functional forms of ${\varepsilon}_k(\eta)$ cannot easily be decided in any detail. However, using the asymptotic form as a guide, it is assumed that the behaviour of the ACER in the tail
is dominated by a function of the form $\exp\{ - a(\eta - b)^c\}$
($\eta \ge \eta_1 \ge b$) where $a$, $b$ and $c$ are suitable
constants, and $\eta_1$ is an appropriately chosen tail marker.
Hence, it will be assumed that,
\begin{align}\label{eq:ag1}
{\varepsilon}_k(\eta) = q_k(\eta) \, \exp\{- a_k (\eta
- b_k)^{c_k} \} \, , \:\: \eta \ge \eta_1,
\end{align}
where the function $q_k(\eta)$ is slowly varying compared with the
exponential function $\exp\{- a_k (\eta - b_k)^{c_k} \}$ and $a_k,
b_k$, and $c_k$ are suitable constants, that in general will be
dependent on $k$. Note that the value $c_k = q_k(\eta) = 1$ corresponds to the
asymptotic Gumbel distribution, which is then a special case of the assumed tail behaviour. And, of course, any extreme value distribution with an ACER function of the form given by Eq.~\eqref{eq:ag1}, is asymptotically Gumbel, cf. comments at the end of Section~\ref{sec:evgp}.

Note that under the assumptions made, a plot of $-\,\log \big|
\log \big({\varepsilon}_k(\eta)/q_k(\eta) \big)\big|$ versus
$\log (\eta - b_k)$ will exhibit a perfectly linear tail
behaviour. This will be illustrated in Chapter~\ref{chap:mc}.

It is realized that if the function $q_k(\eta)$ could be replaced
by a constant value, $q_k$ say, one would immediately be in a
position to apply a linear extrapolation strategy for deep tail prediction problems. In general,
$q_k(\eta)$ is not constant, but its variation in the tail region
is often sufficiently slow to allow for its replacement by a
constant, possibly by adjusting the tail marker $\eta_1$. The proposed statistical approach to the prediction of extreme values is therefore based on the assumption that we can write, 
\begin{align}\label{eq:ag3}
{\varepsilon}_k(\eta) = q_k \, \exp\{- a_k (\eta
- b_k)^{c_k} \} \, , \:\: \eta \ge \eta_1,
\end{align}
where $a_k, b_k, c_k$ and $q_k$ are appropriately chosen constants. In a certain sense, this is a minimal class of parametric functions that can be used for this purpose, which makes it possible to achieve three important goals. Firstly, the parametric class contains the asymptotic form given by $c_k = q_k = 1$ as a special case. Secondly, the class is flexible enough to capture to a certain extent subasymptotic behaviour of any extreme value distribution that is asymptotically Gumbel. Thirdly, the parametric functions agree with a wide range of known special cases, of which a very important example is the extreme value distribution for a regular stationary Gaussian process, which has $c_k=2$.    


The viability of this approach has been successfully demonstrated for extreme value statistics of the response processes related to a wide range of different dynamical systems, cf. \cite{naga08jem,nagaha07oe}. In Chapter~\ref{chap:paeva}, the ACER method will be applied to the problem of predicting extreme wind speeds. The performance of the annual maxima method and the POT method for this purpose, will also be discussed. 

As to the question of finding the parameters $a, b, c, q$ \index{ACER function!parameter estimation}\index{parameter estimation!ACER function} (the subscript $k$, if it applies, is
suppressed), the adopted approach is to determine these parameters by minimizing the following mean square error function with respect to all four arguments,
\begin{equation}\label{eq:ag4}
F(a, b, c, q) = \sum_{j = 1}^J w_j\, \big| \log
\hat{\varepsilon}(\eta_j) - \log q + a (\eta_j - b)^{c} \big|^2\,,
\end{equation}
where $\eta_1 < \ldots < \eta_J$ denotes the levels where the ACER function has been estimated, $w_j$ denotes a weight factor that puts more emphasis on the more reliably estimated $\hat{\varepsilon}(\eta_j)$, cf. the discussion about the BLUE below. However, the choice of weight factor is to some
extent arbitrary. With $\big(C^-(\eta), C^+(\eta)\big)$ denoting the 95\% confidence interval for the value
${\varepsilon}_k(\eta)$, we have previously used $w_j = \big(\log
C^+(\eta_j) - \log C^-(\eta_j) \big)^{-\theta}$ with $\theta = 1$
and 2, combined with a Levenberg-Marquardt least squares
optimization method \citep{gimuwr81}. This has usually worked well
provided reasonable, initial values for the parameters were
chosen. Note that the form of $w_j$ puts some
restriction on the use of the data. Usually, there is a level
$\eta_j$ beyond which $w_j$ is no longer defined, that is, $C^-(\eta_j)$ becomes negative. Hence, the
summation in Eq.~(\ref{eq:ag4}) has to stop before that happens.
Also, the data should be preconditioned by establishing the tail
marker $\eta_1$ based on inspection of the empirical ACER functions.
 

In general, to improve robustness of results, it is recommended to apply a nonlinearly constrained optimization \citep{foho10}. 
The set of constraints is written as,
\begin{equation}\label{eq:ag5}
  \left\{ {\begin{array}{l}
  \log q-a(\eta_i -b)^c < 0 \,, \\ 
  0<q<+\infty \,, \\ 
  \min\limits_j X_j < b \le \eta_1 \,, \\ 
  0 < a < +\infty \,, \\ 
  0 < c < 5 \,. \\ 
 \end{array}} \right.
\end{equation}
Here, the first nonlinear inequality constraint is evident, since under our assumption we 
have $\hat{{\varepsilon}} (\eta_i ) = q\exp \{-a(\eta_i -b)^c\}$, 
and $\hat{{\varepsilon}} (\eta_i ) < 1$ by definition.


A note of caution: When the parameter $c$ is equal to 1.0 or close to it, that is, the distribution is close to the Gumbel distribution, the optimization problem becomes ill-defined or close to ill-defined.
It is seen that when $c=1.0$, there is an infinity of $(b,q)$ values that gives exactly the same value of $F( a, b, c, q)$. Hence, there is no well defined optimum in parameter space. There are simply too many parameters. This problem is alleviated by fixing the $q$-value, and the obvious choice is $q=1$. The restriction $c < 5$ is a practical one and has no real significance beyond limiting the range of $c$-values. A $c$-value larger than $5$ should cause an inspection of the choice of tail marker. For $c$-values larger than 1.0, the $c$-value would typically decrease with increasing tail marker, indicating a deep tail approach to the Gumbel distribution, which is the asymptotic limit. 


Although the Levenberg-Marquardt method generally works well with four or, when appropriate, three parameters, a more direct and transparent optimization method has also been developed
for the problem at hand. It is realized by scrutinizing
Eq.~(\ref{eq:ag4}) that if $b$ and $c$ are fixed, the optimization
problem reduces to a standard weighted linear regression problem.
That is, with both $b$ and $c$ fixed, the optimal values of $a$
and $\log q$ are found using closed form weighted linear
regression formulas in terms of $w_j$, $y_j = \log
{\varepsilon}(\eta_j)$ and $x_j = (\eta_j - b)^{c}$. In that light, it
can also be concluded that the best linear unbiased estimators
(BLUE) are obtained for $w_j = \sigma_{yj}^{-2}$, where
$\sigma_{yj}^2 = \Var[y_j]$ (empirical) \citep{drsm98,mopevi02}.
Unfortunately, this is not a very practical weight factor for the
kind of problem we have here because the summation
in Eq.~(\ref{eq:ag4}) then typically would have to stop at
undesirably small values of $\eta_j$.


It is obtained that the optimal values of $a$ and $q$ are given by
the relations,
\begin{equation}\label{eq:ag6}
a^*(b, c) = - \frac{\sum_{j = 1}^N w_j (x_j - \overline{x})(y_j -
\overline{y})}{\sum_{j = 1}^N w_j (x_j - \overline{x})^2 } \,,
\end{equation}
and
\begin{equation}\label{eq:ag7}
\log q^*(b, c) = \overline{y} + a^*(b,c) \overline{x} \,,
\end{equation}
where $\overline{x} = \sum_{j=1}^N w_j x_j / \sum_{j=1}^N w_j$, with a similar definition of $\overline{y}$.

To calculate the final optimal set of parameters, one may use the
Levenberg-Marquardt method on the function $\tilde{F}(b,c) =
F(a^*(b,c), b, c, q^*(b,c))$ to find the optimal values $b^*$ and
$c^*$, and then use Eqs.~(\ref{eq:ag6}) and (\ref{eq:ag7}) to
calculate the corresponding $a^*$ and $q^*$.


For a simple construction of a confidence interval for the target deep tail extreme value provided by a particular ACER function as given by the fitted parametric curve, the empirical confidence band is reanchored to the fitted curve by centering the individual confidence intervals CI$_{0.95}$ for the point estimates of the ACER function on the fitted curve. Under the premise that the specified class of parametric curves fully describes the behaviour of the ACER functions in the tail, parametric curves are fitted as described above to the boundaries of the reanchored confidence band. These curves are used to determine a first estimate of a 95\% confidence interval for the  target extreme value. To obtain a more precise estimate of the confidence interval, a bootstrapping method would be recommended. A comparison of estimated confidence intervals by both these methods will be presented in Section~\ref{sec:SD} on extreme value prediction for synthetic data. 

As a final point, it has been observed that the predicted value is not very sensitive to the choice of $\eta_1$, provided it is chosen with some care. This property is easily recognized by looking at the way the optimized fitting is done. If the tail marker is in the appropriate domain of the ACER function, the optimal fitted curve does not change appreciably by moving the tail marker.     




\section{Estimation of Extremes for the General Case}\label{sec:Q}\index{ACER function!asymptotic GEV case}


For independent data in the general case, the ACER function ${\varepsilon}_1(\eta)$ can be expressed asymptotically as, 
\begin{align}\label{eq:gev1}
{\varepsilon}_1(\eta) \underset{\eta \to \infty}{\simeq}  \left[ 1 + \gamma \big( a(\eta - b) \big) \right]^{- \frac{1}{\gamma}} \, , 
\end{align}
where $a > 0$, $b$, $\gamma$ are constants. This follows from the explicit form of the generalized extreme value (GEV) distribution, which has been discussed in Chapter~\ref{chap:aevt}. 



Again, the implication of this assumption on the possible subasymptotic
functional forms of ${\varepsilon}_k(\eta)$ in the general case is not a trivial matter. 
The approach we have chosen is to assume that the class of parametric functions needed for the prediction of extreme values for the general case can be modelled on the relation between the Gumbel distribution and the GEV distribution. While the extension of the asymptotic Gumbel case to the proposed class of subasymptotic distributions was fairly transparent, this is not equally so for the general case. However, using a similar kind of approximation, the behaviour of the mean exceedance rate in the subasymptotic part of the tail
is assumed to follow a function largely of the form $\left[ 1  +  \gamma \big( a(\eta - b)^c \big) \right]^{- \frac{1}{\gamma}}$
($\eta \ge \eta_1 \ge b$) where $a > 0$, $b$, $c > 0$ and $\gamma > 0$ are suitable
constants, and $\eta_1$ is an appropriately chosen tail level.
Hence, it will be assumed that \citep{naes10p14},
\begin{align}\label{eq:gev2}
{\varepsilon}_k(\eta) = q_k(\eta) \, \left[ 1 + \gamma_k \big( a_k(\eta - b_k)^{c_k} \big) \right]^{- \frac{1}{\gamma_k}} \, , \:\: \eta \ge \eta_1,
\end{align}
where the function $q_k(\eta)$ is weakly varying compared with the
function \linebreak $\left[ 1 + \gamma_k \big( a_k(\eta - b_k)^{c_k} \big)\right]^{- \frac{1}{\gamma_k}}$ and $a_k > 0$, $b_k$, $c_k > 0$ and $\gamma_k > 0$ are suitable constants, that in general will be
dependent on $k$. Note that the values $c_k = 1$ and $q_k(\eta) = 1$ correspond to the
asymptotic limit, which is then a special case of the general expression given in Eq.~\eqref{eq:gev2}. Another method to account for subasymptotic effects has been proposed by \cite{eata12}, building on ideas developed by \cite{tawn90}, \cite{leta96} and \cite{heta04}. In this approach, the asymptotic form of the marginal distribution of exceedances is kept, but it is modified by a multiplicative factor accounting for the dependence structure of exceedances within a cluster.



An alternative form to Eq.~\eqref{eq:gev2} would be to assume that,  
\begin{align}\label{eq:emcr5a}
{\varepsilon}_k(\eta) =  \, \left[ 1 + \gamma_k \big( a_k(\eta - b_k)^{c_k} + d_k(\eta)  \big) \right]^{- \frac{1}{\gamma_k}} \, , \:\: \eta \ge \eta_1,
\end{align}
where the function $d_k(\eta)$ is weakly varying compared with the
function $a_k(\eta - b_k)^{c_k}$. However, for estimation purposes, the form given by Eq.~\eqref{eq:gev2} appears to be preferable as it leads to somewhat simpler estimation procedures.



For practical identification of the ACER functions given by Eq.~\eqref{eq:gev2}, it is expedient to assume that the unknown function $q_k(\eta)$ varies sufficiently slowly to be replaced by a constant. In general,
$q_k(\eta)$ is not constant, but its variation in the tail region
is assumed to be sufficiently slow to allow for its replacement by a
constant. Hence, as in the Gumbel case, it is in effect assumed that $q_k(\eta)$ can be replaced by a constant for $\eta \ge \eta_1$, for an appropriate choice of tail marker $\eta_1$. For simplicity of notation, in the following we shall suppress the index $k$ on the ACER functions, which will then be written as,
\begin{align}\label{eq:gev3}
{\varepsilon}(\eta) = q \, \left[ 1 + \tilde{a} \, (\eta - b)^{c} \right]^{- \xi} \, , \:\: \eta \ge \eta_1,
\end{align} 
where $\xi = 1/ \gamma$, $\tilde{a} = a \gamma$.

For the analysis of data, first the tail marker $\eta_1$ is provisionally identified from visual
inspection of the log plot $(\eta, \ln\hat{\varepsilon}(\eta))$. The value chosen for $\eta_1$
corresponds to the beginning of regular tail behaviour in a sense
to be discussed below. 

The optimization process to estimate the parameters\index{parameter estimation!ACER function}\index{ACER function!parameter estimation} is done relative to the log plot, as for the Gumbel case.
The mean square error function to be minimized is in the general case written as,
\begin{equation}\label{eq:lm2}
 F( \tilde{a}, b, c, q, \xi) =  \sum_{j = 1}^N w_j\, \big| \log \hat{\varepsilon}(\eta_j) - \log q + 
\xi \log \left[ 1  +  \tilde{a}(\eta_j - b)^c \right]\big|^2\,,
\end{equation}
where $w_j$ is a weight factor as previously defined. 

An option for estimating the five parameters $\tilde{a}, b, c, q, \xi$ is again to use the Levenberg-Marquardt least squares optimization method, which can be simplified also in this case by observing
that if $\tilde{a}$, $b$ and $c$ are fixed in Eq.~(\ref{eq:lm2}), the optimization
problem reduces to a standard weighted linear regression problem.
That is, with $\tilde{a}$, $b$ and $c$ fixed, the optimal values of $\xi$
and $\log q$ are found using closed form weighted linear
regression formulas in terms of $w_j$, $y_j = \log
\hat{\varepsilon}(\eta_j)$ and $x_j = \log(1 + \tilde{a}(\eta_j - b)^{c})$. 

It is obtained that the optimal values of $\xi$ and $\log q$ are given by
relations similar to Eqs.~\eqref{eq:ag6} and \eqref{eq:ag7}.
To calculate the final optimal set of parameters, the
Levenberg-Marquardt method may then be used on the function $\tilde{F}(\tilde{a}, b, c) =
F( \tilde{a}, b, c, q^*(\tilde{a}, b, c), \xi^*(\tilde{a}, b, c))$ to find the optimal values $\tilde{a}^*$, $b^*$ and $c^*$, and then the corresponding $\xi^*$ and $q^*$ can be calculated. The optimal values of the parameters may e.g also be found by a sequential quadratic programming (SQP) method~\citep{nag}. This general case formulation of the ACER method used on financial risk issues, which often requires heavy tailed distributions, is discussed in Chapter~\ref{chap:eevra}. 

For the construction of confidence intervals, the procedure discussed at the end of the previous section can easily be adapted to the present case.





