[1,2]Justus Contzen

1]Section Paleoclimate Dynamics, Alfred Wegener Institute Helmholtz Center for
Polar and Marine Research, Bremerhaven, Germany

2]Department of Environmental Physics, University of Bremen, Bremen, Germany

3]Institute for Statistics, University of Bremen, Bremen, Germany

# Regionalization approaches for the spatial analysis of extremal dependence

<EMAIL_ADDRESS>Thorsten Dickhaus  Gerrit Lohmann  [  [  [

###### Abstract

The impact of an extreme climate event depends strongly on its geographical
scale. Max-stable processes can be used for the statistical investigation of
climate extremes and their spatial dependencies on a continuous area. Most
existing parametric models of max-stable processes assume spatial stationarity
and are therefore not suitable for the application to data that cover a large
and heterogeneous area. For this reason, it has recently been proposed to use
a clustering algorithm to divide the area of investigation into smaller
regions and to fit parametric max-stable processes to the data within those
regions. We investigate this clustering algorithm further and point out that
there are cases in which it results in regions on which spatial stationarity
is not a reasonable assumption. We propose an alternative clustering algorithm
and demonstrate in a simulation study that it can lead to improved results.

###### keywords:

Max-stable Processes, Clustering, Extremal Dependence, Multivariate Extremes

##  1 Introduction

Extreme climate events like droughts, heatwaves or heavy rainfall events
usually do not occur in isolation at a certain observed location, but in
combination with more or less extreme events in the vicinity. Having
information on how strongly extremes at different locations influence each
other is of high interest as the consequences of large-scale extreme weather
are often especially severe. Resources to mitigate the impacts of such events
might be limited, while their effects on the economy, the environment or the
healthcare system can be drastic. When analyzing climate data in a certain
region with regard to extreme events, it is therefore important to focus not
only on univariate distributions (i.e. distributions at one location alone),
but also on their interdependence and on the multivariate distribution of the
data.

One common approach in extreme value analysis is the investigation of the
block-wise maxima of given time series (for example, annual maxima of daily
data). For the spatial analysis of data in a certain region, stochastic
processes on a compact subset of $\operatorname{{\mathbb{R}}}^{2}$ with almost
surely continuous sample paths are often employed, which makes it possible to
investigate spatial dependencies not only for the locations of the spatial
data at hand, but also for unobserved locations in the area. Under mild
regularity conditions, spatial block-wise maxima can be described using max-
stable processes. Using a spectral representation, parametric sub-classes of
max-stable processes can be derived which allow for parametric inference. A
composite maximum likelihood estimator can be used to fit such processes to
given data. Different sub-classes of max-stable processes have been proposed
[1, 2, 3, 4], but they all describe processes that are spatially stationary,
i.e. their finite-dimensional marginal distributions are invariant with
respect to translations in space. Such an assumption may be reasonable if the
investigated area is small and rather homogeneous, and indeed such models have
successfully been applied to climate extremes in Switzerland [5]. For the
application to data from a larger area of investigation they are not well
suited, and non-stationary extensions of them are an active research topic.
[6] propose such a max-stable model based on non-stationary Gaussian
processes. Their model flexibly captures non-stationarity, but it relies on
covariates which are not given for every data set. Another approach to
investigate large-scale data using max-stable processes is due to [7]. They
employ a clustering algorithm by [8] to split the investigated area into
several small regions. To each of the smaller regions, a stationary max-stable
process is then fitted.

In this work, we will build up on the aforementioned clustering approach.
While clustering is a promising idea to make stationary max-stable processes
applicable, we will see that the algorithm by [8] does not necessarily produce
regions in which spatial stationarity is a reasonable assumption. We will
illustrate this using a simple example of a non-stationary max-stable process,
and we will propose an alternative clustering algorithm and compare the two in
a simulation study.

The rest of this paper is structured as follows: In the next section, we will
give a brief introduction to max-stable processes. The clustering algorithm by
[7] and our suggested algorithm are discussed in Section 3. Using a simulation
study, we investigate the performance of the two algorithms in Section 4.
Section 5 on conclusions and a discussion finalize the article.

##  2 Theoretical foundations

This section gives a brief overview of the theory of max-stable processes, for
a more detailed treatment see [9] and [5]. Let $\operatorname{{\mathcal{S}}}$
be a compact subset of $\operatorname{{\mathbb{R}}}^{2}$. Throughout the
paper, we will investigate stochastic processes on
$\operatorname{{\mathcal{S}}}$ with almost surely continuous sample paths and
we use the short notation $Y$ for
$\\{Y(s)\\}_{s\in\operatorname{{\mathcal{S}}}}$. The distribution of such a
process is completely determined by its finite-dimensional marginal
distributions.

A stochastic process $Y$ is called max-stable if for all
$n\in\operatorname{{\mathbb{N}}}$ there exist continuous functions
$\\{c_{n}(s)\\}_{s\in\operatorname{{\mathcal{S}}}}$ and
$\\{d_{n}(s)\\}_{s\in\operatorname{{\mathcal{S}}}}$ with $c_{n}(s)>0$ such
that for independent copies $Y_{1},\dots,Y_{n}$ of $Y$

$\displaystyle\Big{\\{}\frac{\max_{i=1,\dots,n}Y_{i}(s)-d_{n}(s)}{c_{n}(s)}\Big{\\}}_{s\in\operatorname{{\mathcal{S}}}}\overset{\mathcal{D}}{=}\\{Y(s)\\}_{s\in\operatorname{{\mathcal{S}}}}$
(1)

holds, with $\overset{\mathcal{D}}{=}$ denoting equality in distribution. For
$n\in\operatorname{{\mathbb{N}}}$, let $X_{1},\dots,X_{n}$ be independent and
identically distributed stochastic processes and define
$M_{(n)}(s):=\max_{i=1,\dots,n}(X_{i}(s))$. We assume the existence of
functions $a_{n}(\cdot)>0$ and $b_{n}(\cdot)$ on
$\operatorname{{\mathcal{S}}}$ such that $M_{(n)}$ converges in distribution
($\xrightarrow{\mathcal{D}}$) to some process $Y$ as $n$ tends to infinity:

$\displaystyle\Big{\\{}\frac{\\{M_{(n)}(s)-b_{n}(s)}{a_{n}(s)}\Big{\\}}_{s\in\operatorname{{\mathcal{S}}}}\xrightarrow{\mathcal{D}}\\{Y(s)\\}_{s\in\operatorname{{\mathcal{S}}}}.$
(2)

Then, $Y$ must be a max-stable process if it is not degenerate. In practical
applications, the data are grouped into blocks of a fixed block size $n$, and
it is assumed that a max-stable process is a reasonable approximation to the
resulting block-wise maxima. The margins of a max-stable process are
necessarily max-stable random variables, and it is well known that as such,
they must follow a generalized extreme value (GEV) distribution. The GEV
distributions form a parametric family for which statistical inference is
well-established and several methods of parameter estimation exist [see e.g.
10, Chapter 7]. Using marginal transformations, it can be assumed without loss
of generality that all marginal distributions of the investigated max-stable
processes are identical (note that the GEV distributions are absolutely
continuous), and it is common to consider margins that are unit Fréchet
distributed (cdf $\Phi$ given by $\Phi(z)=\exp(-z^{-1})$, $z>0$). A max-stable
process with unit Fréchet margins is called a simple max-stable process.

The finite-dimensional marginal distribution of a simple max-stable process
$Y$ at locations $s_{1},\dots,s_{n}\in\operatorname{{\mathcal{S}}}$,
$n\in\operatorname{{\mathbb{N}}}$, can be written as

$\displaystyle\Pr(Y(s_{1})\leq y_{1},...,Y(s_{n})\leq
y_{n})=\exp(-V_{s_{1},\dots,s_{n}}(y_{1},\dots,y_{n}))$ (3)

for a function $V_{s_{1},\dots,s_{n}}$ that fulfills
$V_{s_{1},\dots,s_{n}}(ay_{1},\dots,ay_{n})=a^{-1}V_{s_{1},\dots,s_{n}}(y_{1},\dots,y_{n})$
for all $a>0$ [9]. In practical applications, a max-stable process is fitted
to data that are given on a finite set of locations
$\operatorname{{\mathcal{X}}}\subseteq\operatorname{{\mathcal{S}}}$ (for
example the locations of weather stations). The data at each location have
been transformed to a unit Fréchet distribution. Denote the data at location
$x\in\operatorname{{\mathcal{X}}}$ by $y_{x}^{(1)},\dots,y_{x}^{(m)}$ with
$m\in\operatorname{{\mathbb{N}}}$ the sample size. Below, we will present
several parametric models for max-stable processes, but before doing so, we
will discuss how to fit parametric models to data. Assume a parametric
subclass of max-stable processes with a parameter space $\Psi$, which implies
in particular that the functions $V_{s_{1},\dots,s_{n}}$ are parametrized as
$V_{s_{1},\dots,s_{n}}^{\psi}$ with $\psi\in\Psi$ for all
$s_{1},\dots,s_{n}\in\operatorname{{\mathcal{S}}}$. The multivariate density
of a max-stable process $Y$ at $s_{1},\dots,s_{n}\in S$ can be derived from
Eq. (3) by taking the mixed derivative $\partial s_{1}\dots\partial s_{n}$. It
contains a sum of products of partial derivatives of
$V_{s_{1},\dots,s_{n}}^{\psi}$. Even for a moderate value of $n$ like $n=10$,
the number of summands in the density is so large that numerical maximization
is not computationally feasible. The classical maximum-likelihood approach is
therefore usually not applicable to the multivariate density of all data
points in $\operatorname{{\mathcal{X}}}$ [5]. As a remedy, it is common to use
a composite likelihood approach instead: Maximum-likelihood estimation is not
carried out by maximizing the full log-likelihood, but instead by maximizing
the sum of bivariate marginal log-likelihoods

$\displaystyle\hat{\psi}=\operatorname*{arg\,max}_{\psi\in\Psi}L(\psi)$ (4)

with

$\displaystyle L(\psi)=\sum_{i=1}^{m}\sum_{x_{1}\neq
x_{2}\in\operatorname{{\mathcal{X}}}}\log
f^{\psi}_{x_{1},x_{2}}(y_{x_{1}}^{(i)},y_{x_{2}}^{(i)}).$ (5)

The bivariate densities can be calculated as

$f^{\psi}_{x_{1},x_{2}}(y_{1},y_{2})=\exp\Big{(}\\!-V^{\psi}_{x_{1},x_{2}}(y_{1},y_{2})\Big{)}\,\cdot\\\
\qquad\Big{(}\frac{\partial}{\partial
y_{1}}V^{\psi}_{x_{1},x_{2}}(y_{1},y_{2})\frac{\partial}{\partial
y_{2}}V^{\psi}_{x_{1},x_{2}}(y_{1},y_{2})-\frac{\partial^{2}}{\partial
y_{1}\partial y_{2}}V^{\psi}_{x_{1},x_{2}}(y_{1},y_{2})\Big{)}.$ (6)

To reduce numerical complexity, it is also common practice to include only
those pairs $(x_{1},x_{2})$ in Eq. (5) for which $\|x_{1}-x_{2}\|$ does not
exceed a certain threshold (throughout the paper, $\|\cdot\|$ denotes the
Euclidean norm).

As mentioned in the introduction, we will use clustering algorithms to split
the set $\operatorname{{\mathcal{X}}}$ into different subsets and and we will
fit max-stable processes to the data on each subset. Introducing the notation
for this, let $M\subseteq\operatorname{{\mathcal{X}}}$ be a subset of
$\operatorname{{\mathcal{X}}}$ with at least two elements, then we write the
corresponding composite likelihood as

$\displaystyle L_{M}(\psi)=\sum_{i=1}^{m}\sum_{x_{1}\neq x_{2}\in M}\log
f^{\psi}_{x_{1},x_{2}}(y_{x_{1}}^{(i)},y_{x_{2}}^{(i)})$ (7)

and the corresponding estimator as

$\displaystyle\hat{\psi}_{M}=\operatorname*{arg\,max}_{\psi\in\Psi}L_{M}(\psi).$
(8)

Parametric subfamilies of max-stable processes are constructed using the
following stochastic representation [11, 12]: Let $Y$ be a simple max-stable
process. Then,

$\displaystyle\Big{\\{}Y(s)\Big{\\}}_{s\in\operatorname{{\mathcal{S}}}}\overset{\mathcal{D}}{=}\Big{\\{}\max_{i\geq
1}\zeta_{i}Z^{(i)}(s)\Big{\\}}_{s\in\operatorname{{\mathcal{S}}}},$ (9)

where $\\{\zeta_{i}\\}_{i\in\operatorname{{\mathbb{N}}}}$ are the points of a
Poisson point process with intensity measure
$d\Lambda(\zeta)=\zeta^{-2}d\zeta$ and
$\\{Z^{(i)}\\}_{i\in\operatorname{{\mathbb{N}}}}$ are independent copies of a
non-negative stochastic process $Z$ fulfilling
$\operatorname{{\mathbb{E}}}[Z(s)]=1$ for all
$s\in\operatorname{{\mathcal{S}}}$. The process $Z$ is called spectral process
of $Y$.

Thus, the distribution of $Y$ can be modeled by choosing a model for $Z$. One
of the first parametric subfamilies proposed was the Schlather model [1],
using the underlying process $Z$ given by

$\displaystyle Z(s)=\sqrt{2\pi}\max(0,G(s))$ (10)

with $G$ a standard Gaussian process with a spatially stationary covariance
function $\rho(h)$ (depending only on
$h=s_{1}-s_{2}\in\operatorname{{\mathbb{R}}}^{2}$,
$s_{1},s_{2}\in\operatorname{{\mathcal{S}}}$). [4] showed that the Schlather
process has a limited scope of applicability because its bivariate
distributions always exhibit extremal dependence, even if the underlying
Gaussian variables are uncorrelated. To avoid this limitation, he extended
Schlather’s model to the extremal-t model, using the spectral process $Z$
given by

$\displaystyle Z(s)=c_{\nu}\max(0,G(s))^{\nu}$ (11)

with $\nu\geq 1$ fixed, $G$ as in the Schlather model and $c_{\nu}$ a constant
ensuring that $\operatorname{{\mathbb{E}}}[Z(s)]=1$. For this process, it
holds

$V_{s,s+h}(y_{1},y_{2})=\frac{1}{y}_{1}T_{\nu+1}\Big{(}\frac{(y_{2}/y_{1})^{1/\nu}-\rho(h)}{(\nu+1)^{-1/2}(1-\rho(h)^{2})^{1/2}}\Big{)}+\\\
\qquad\qquad\frac{1}{y}_{2}T_{\nu+1}\Big{(}\frac{(y_{1}/y_{2})^{1/\nu}-\rho(h)}{(\nu+1)^{-1/2}(1-\rho(h)^{2})^{1/2}}\Big{)}$
(12)

with $T_{\nu+1}$ the cdf of the Student-t distribution with $(\nu+1)$ degrees
of freedom. Bivariate densities can be calculated based on this equation.

One possible choice for $\rho$ is $\rho(h)=\exp(-\|h/\lambda\|^{\alpha})$ with
$\lambda>0$ and $0<\alpha\leq 2$ [5]. To model anisotropic behavior, we follow
[13] and use

$\displaystyle\rho(h)=\exp(-\|Ah\|^{\alpha})$ (13)

instead, with $A$ a $2\times 2$ transformation matrix, meaning that

$\displaystyle A=\begin{pmatrix}\sin(\gamma)/a&\cos(\gamma)/(a+b)\\\
-\cos(\gamma)/(a+b)&\sin(\gamma)/b\end{pmatrix}$ (14)

with parameters $a>0$, $b\geq 0$ and $\gamma\in[0,\pi)$. The level sets of
$\rho$ are then ellipses with $a$ and $a+b$ proportional to the lengths of the
minor and major axes and $\gamma$ the angle between the horizontal axis and
the major axis of the ellipsis. In climate data, elliptical level sets are
preferable over the circular ones of isotropic covariance functions because
the presence of ocean currents, a predominant wind direction or topographical
boundaries frequently cause the extremal dependence structure of two points to
depend not only on their distance, but also on the direction.

##  3 Clustering algorithms

###  3.1 Extremal Dependence Clustering

As mentioned in the introduction, one approach to investigate data for which
the assumption of spatial stationarity is not reasonable is using a clustering
algorithm that splits the area of investigation into smaller regions.
Stationary max-stable processes like the extremal-t process can then be fitted
to the data within those regions. This approach is due to [7], and they
perform the regionalization using a hierarchical clustering algorithm (see
their paper for an introduction to hierarchical clustering). To apply a
hierarchical clustering algorithm, a dissimilarity measure for all pairs of
elements in $\operatorname{{\mathcal{X}}}$ is required. A dissimilarity
measure is a non-negative symmetric function
$D:\operatorname{{\mathcal{X}}}\times\operatorname{{\mathcal{X}}}\rightarrow\operatorname{{\mathbb{R}}}$
that fulfills $D(x,x)=0$ for all $x\in\operatorname{{\mathcal{X}}}$. [7]
propose to use a dissimilarity measure first developed by [8] that is based on
pairwise extremal coefficients, which in turn are defined as

$\theta_{x_{1},x_{2}}=V_{x_{1},x_{2}}(1,1).$

The extremal coefficient is a useful summary measure for the dependency of
$Y(x_{1})$ and $Y(x_{2})$. It takes values between one and two, with a value
of one corresponding to the variables being comonotonic and a value of two
corresponding to them being stochastically independent, and it can be
estimated from a data sample by using the madogram estimator $\hat{\theta}$ by
[14] and [15]. Defining

$D_{1}(x_{1},x_{2}):=\hat{\theta}_{x_{1},x_{2}}-1$ (15)

yields a dissimilarity function as required. In the following, we will use the
term Extremal Dependence Clustering (EDC) for the clustering based on this
dissimilarity function.

###  3.2 Applicability of the Extremal Dependence Clustering

The EDC algorithm is performed with the goal of defining regions that are
suitable for fitting spatially stationary max-stable processes to the data
within them. It should therefore group points together in such a way that
within the clusters stationarity can be assumed, or is at least a reasonable
approximation. The dissimilarity measure by [8] is based on the comparison of
extremal coefficients and it therefore groups together points with a tendency
for concurrent extremes. Within the resulting regions, pairwise extremal
dependencies tend to be high in general, which might reduce the possible
extent of spatial non-stationarity. Nevertheless, spatial stationarity is not
a justified assumption within the clusters defined that way, and the
dissimilarity measure by [8] was not designed with the intention of finding
such regions.

We illustrate this using a concrete example of a non-stationary max-stable
process. We can construct such a process using an approach by [6]. They
extended the stationary extremal-t process by using a non-stationary Gaussian
process as underlying spectral process: Instead of one $2\times 2$
transformation matrix $A$ for the whole space, they use for each point
$s\in\operatorname{{\mathcal{S}}}$ a matrix $A_{s}$ such that the map
$s\mapsto A_{s}$ is continuous. Let as before $0<\alpha\leq 2$ be fixed. Using
the notations $\Omega_{s}=(A_{s}^{T}A_{s})^{-1}$ with $A^{T}$ denoting the
transpose of $A$ and $R_{\alpha}(x)=\exp(-x^{\alpha})$ for
$x\in\operatorname{{\mathbb{R}}}_{\geq 0}$, they show how to construct a non-
stationary Gaussian process using kernel convolution [16]. The covariance
structure of the resulting process is given by

$\rho(s_{1},s_{2})=|\Omega_{s_{1}}|^{\frac{1}{4}}|\Omega_{s_{2}}|^{\frac{1}{4}}\Big{|}\frac{\Omega_{s_{1}}+\Omega_{s_{2}}}{2}\Big{|}^{-\frac{1}{2}}\cdot\\\
R_{\alpha}\Big{(}\sqrt{({s_{1}}-{s_{2}})^{T}\Big{(}\frac{\Omega_{s_{1}}+\Omega_{s_{2}}}{2}\Big{)}^{-1}({s_{1}}-{s_{2}})}\Big{)},$
(16)

which reduces to the stationary extremal-t process from the previous section
if $A_{s}$ is constant on the whole space.

As a simple example of a non-stationary process we use a Huser-Genton process
on the set $\operatorname{{\mathcal{S}}}=[-5,5]\times[-5,5]$ with matrix
parameters $a_{s}=2$ constant, $b_{s}=(x+5)/2$,
$s=(x,y)\in\operatorname{{\mathcal{S}}}$ and $\gamma_{s}=0$ constant. The
global model parameters are $\nu=5$ and $\alpha=1$. This process is obviously
stationary on the sets $\\{x\\}\times[-5,5]$ for all $x\in[-5,5]$, and if we
investigate a vertical stripe of the form
$[x-\epsilon,x+\epsilon]\times[-5,5]$, $\epsilon>0$ small, the values of
$a_{s}$, $b_{s}$ and $\gamma_{s}$ in that region are very similar and
stationarity is a reasonable approximation. The clustering based on extremal
coefficients, however, does not result in clusters of such a form. In Fig. 1
we depict the pairwise extremal coefficients $\theta_{s,t}$,
$t\in\operatorname{{\mathcal{S}}}$ for four selected values of $s$:
$s_{1}=(-3,2)$, $s_{2}=(3,2)$, $s_{3}=(-3,-2)$ and $s_{4}=(3,-2)$. It can be
observed also from Fig. 1 that the dependence structures for the points with
the same value of $x$ are identical and that it would therefore be reasonable
to group them into the same cluster. However, the extremal coefficient
$\theta_{s_{1},s_{2}}$ is close to $2$, so the points $s_{1}$ and $s_{2}$ will
likely not be grouped into the same cluster by the EDC clustering. The same
holds for the points $s_{3}$ and $s_{4}$. Instead, pairs of points with a low
extremal coefficient, like for example $(2,2)$ and $(3,2)$ will be grouped
together even though the dependency structures around these points differ.
Indeed, if we apply the EDC clustering to the true values of the extremal
coefficients, we obtain the clusters shown in Fig. 2, confirming the
theoretical considerations we just made.

Figure 1: Illustration of the non-stationary max-stable process on
$\operatorname{{\mathcal{S}}}=[-5,5]\times[-5,5]$ given by the matrix
parameters $a_{s}=2$, $b_{s}=(x+5)/2$,
$s=(x,y)\in\operatorname{{\mathcal{S}}}$ and $\gamma_{s}=0$ and the global
parameters $\nu=5,\alpha=1$. Displayed are the pairwise extremal coefficients
for the points in $\operatorname{{\mathcal{S}}}$ relative to (a)
$s_{1}=(-3,2)$, (b) $s_{2}=(3,2)$ , (c) $s_{3}=(-3,-2)$ and (d)
$s_{4}=(3,-2).$ Figure 2: The results of the EDC clustering algorithm, applied
to the true extremal coefficients of the non-stationary max-stable process
from Fig. 1 using a number of clusters of five.

###  3.3 Clustering based on local estimates

An approach that is expected to be more suitable to finding spatially
stationary regions is based on the direct comparison of the structures of
extremal dependence around two points. We assume that for each point
$s\in\operatorname{{\mathcal{S}}}$, the extremal coefficients form a locally
elliptic structure, that is, for all points $t\in\operatorname{{\mathcal{S}}}$
in the vicinity of $s$ the extremal coefficient $\theta_{s,t}$ can be
approximated by $R_{\alpha}(\|A_{s}(s-t)\|)$ for some transformation matrix
$A_{s}$ and some $\alpha\in(0,2]$ (fixed on the whole process), with
$R_{\alpha}(x)=\exp(-x^{\alpha})$ as in Section 3.2. This is obviously true
for stationary extremal-t processes and [6] show that it is also true for the
non-stationary processes they designed, so we do not consider this to be a too
severe restriction. Fixing a small $\epsilon>0$ and choosing values for the
global parameters $\alpha$ and $\nu$, we obtain for each of the locations
$x\in\operatorname{{\mathcal{X}}}$ at which data are given an estimate for
$A_{x}$ by using composite maximum likelihood on the sum of the pairwise log-
likelihoods for all pairs $(x,u)$, $u\in
U(x):=\\{u\in\operatorname{{\mathcal{X}}}\big{|}\,\|x-u\|<\epsilon\\}$.
Depending on the spatial structure of $\operatorname{{\mathcal{X}}}$, some
isolated points may have to be excluded beforehand to ensure that $U(x)$ is
always of a sufficient size.

In order to identify structures in the estimates better and to reduce the
influence of outliers, we apply spatial smoothing to the estimated parameters
$a,b$, and $\gamma$ at every point. In our application, to smooth the values
at $x\in\operatorname{{\mathcal{X}}}$ we used local averages on a vicinity of
$x$. Other spatial smoothing algorithms, for example kernel smoothing, are
also possible [see for example 17]. This results in an estimated matrix
$\hat{A}_{x}$ for each location $x\in\operatorname{{\mathcal{X}}}$. If we
group points with similar values of $\hat{A}_{x}$ into one cluster, it is
reasonable to assume that within this cluster spatial stationarity holds
approximately.

To determine which of the matrices are considered ’similar’, we compare the
shapes of the elliptical level sets corresponding to them using the Jaccard
index. Let
$L_{s}:=\\{h\in\operatorname{{\mathbb{R}}}^{2}\big{|}\,||\hat{A}_{s}h||\geq
0.5\\}$. We define a dissimilarity measure as

$\displaystyle D_{2}(x_{1},x_{2})=1-\frac{|L_{x_{1}}\cap
L_{x_{2}}|}{|L_{x_{1}}\cup L_{x_{2}}|}$ (17)

with $|\cdot|$ denoting the area. The value of $0.5$ in the definition of
$L_{x}$ is arbitrary, any other value in $(0,1)$ would yield the same result.
We will use the term Local Estimates Clustering (LEC) for this clustering from
now on.

###  3.4 Comparison of clusterings

To compare the two algorithms, we investigate which of the resulting
clusterings is better suited to fitting stationary max-stable processes to it.
Assume that we have fitted both clustering algorithms to some data and that we
have fitted a stationary max-stable process to each of the resulting clusters.
Denote the clusters of the EDC algorithm by
$\\{\sigma_{1},\dots,\sigma_{p}\\}$ and the clusters of the LEC algorithm by
$\\{\tau_{1},\dots,\tau_{q}\\}$. Each of these elements is a subset of
$\operatorname{{\mathcal{X}}}$, the different clusters in one clustering are
disjoint and their union is the whole set $\operatorname{{\mathcal{X}}}$.
Define $\mathcal{M}=\\{\sigma_{1},\dots,\sigma_{p},\tau_{1},\dots,\tau_{q}\\}$
as the set of all clusters. Remember that for each cluster $M\in\mathcal{M}$,
the estimated parameters of the corresponding max-stable process
$\hat{\psi}_{M}$ have been calculated by maximizing the composite likelihood
$L_{M}$ from Eq. (7).

Note that for either clustering algorithm, max-stable processes are fitted
only to the data within the same cluster, so if two points fall into two
different clusters, a statistical model for their dependency is not provided.
For this reason, it is not possible to compute or compare the composite
likelihoods on the whole set $\operatorname{{\mathcal{X}}}$ from Eq. (5),
which would otherwise be the standard approach for evaluating the goodness of
fit. The likelihoods of the fitted processes of the different clusters are not
comparable to each other because each one of them is based on different
underlying data. However, it is possible to calculate composite likelihoods on
the intersections of clusters of the two clusterings, that is, on the sets
$\upsilon_{ij}=\sigma_{i}\cap\tau_{j}$, provided they contain two or more
elements. On each of these intersections a stationary max-stable process has
been fitted for both algorithms, and by comparing the likelihoods the goodness
of fit of the processes on this area can be compared. This leads to a measure
for the goodness of fit on $\upsilon_{ij}$ for both algorithms:

$\displaystyle\hat{L}_{i,j}^{EDC}=L_{\upsilon_{ij}}(\hat{\psi}_{\sigma_{i}}),\qquad\qquad\hat{L}_{i,j}^{LEC}=L_{\upsilon_{ij}}(\hat{\psi}_{\tau_{j}}).$
(18)

Note that both models have the same number of parameters, so we can compare
the likelihoods directly and do not need a penalty term as in the Akaike or
Bayesian Information Criterion.

##  4 Simulation study

In this section, we compare the two clustering algorithms by means of a
simulation study. To do this, we simulate data from the Huser-Genton model we
already investigated in Section 3.2. Remember that we use for this model as
global parameters $\alpha=1$ and $\nu=5$ and as parameters for the local
dependencies $a_{s}=2$, $b_{s}=(x+5)/2$ and $\gamma_{s}=0$ for
$s=(x,y)\in\operatorname{{\mathcal{S}}}$ (see Fig. 1). We choose a horizontal
and vertical resolution of the space
$\operatorname{{\mathcal{S}}}=[-5,5]\times[-5,5]$ equal to $0.2$ and simulate
data from processes with $250$ independent observations.

For the clustering algorithms, we choose a number of clusters equal to five in
both algorithms. In a first investigation, we apply the algorithms using as
global parameters the true values of $\nu$ and $\alpha$. In Fig. 3a and Fig.
3b, the clusters produced by the two algorithms are displayed. Stationary max-
stable processes are fitted to the data in the clusters, and the color inside
each cluster in Fig. 3a and Fig. 3b depicts the value of the corresponding
estimate for the parameter $b$. For reference, the true values of the
parameter $b_{s}$ are depicted in Fig. 3c. The true values of the other two
parameters $a_{s}$ and $\gamma_{s}$ are constant over the whole space; their
estimates are also similar for all clusters and are not shown. It can be
observed that the clusters of the EDC algorithm are similar to those derived
when applying the EDC algorithm using the true values (Fig. 2). In particular,
as in the theoretical case, there is considerable variation in the true values
of $b_{s}$ within some of the clusters. A fitted stationary process cannot
account for that variation. The LEC algorithm results in clusters that form
vertical stripes, and on these clusters there is less variation in the true
values of $b_{s}$. The fitted values on the clusters are therefore often
closer to the true values than for the EDC algorithm (compare Fig. 3b and Fig.
3c).

Figure 3: Results of the EDC algorithm (a) and the LEC algorithm (b) applied
to simulated data of a Huser-Genton process as in Fig. 1 with $250$
observations. The number of clusters used is five. The colors within the
clusters show the estimated value for parameter $b$ on the stationary-max-
stable process that has been fitted to the data on the cluster. For reference,
the true values of the parameter $b_{s}$ used to simulate the data are shown
in (c).

The two clustering algorithms are compared using the method described in
Section 3.4. In Fig. 4, we depict the intersections of the clusters of the two
algorithms. The color of each region indicates which algorithm has the better
goodness of fit there (darker color — EDC, lighter color — LEC). It can be
observed that the LEC algorithm results in a better goodness of fit on most
regions.

Figure 4: Comparison of the goodness of fit of the two clustering algorithms.
For each intersection of the clusters of the LEC and the EDC algorithm from
Fig. 3, the goodness of fit of the two algorithms is compared. Dark blue color
indicates that the EDC algorithm has a better goodness of fit, light blue
colors indicates a better goodness of fit of the LEC algorithm.

The above analyses present the results for just one simulation and using the
true values for $\nu$ and $\alpha$. For a more general investigation, the
simulation of a non-stationary process and its investigation is repeated
another $24$ times. As the true values of $\nu$ and $\alpha$ are not known in
practical applications, we conduct the analyses also for other parameter
values, using for $\nu$ the values $3,5,$ and $7$ and for $\alpha$ the values
$0.7,1.0,1.3$. The resulting clusters are of course slightly different each
time, but the general structures that can be identified in Fig. 3 stay the
same (not shown). In Fig. 5, for each combination of the values for $\nu$ and
$\alpha$, we depict for each point $s\in\operatorname{{\mathcal{S}}}$ the
percentage of the $24$ simulations for which the LEC algorithms has a better
goodness of fit on the cluster the point is in. It can be observed that for
each choice of the global parameter values and throughout the area of
investigation, the LEC algorithm exhibits at each point a better goodness of
fit for more than $75\%$ of the realizations.

Figure 5: The analysis shown in Fig. 4 is repeated for $24$ independent
simulations and the clustering algorithms are applied using different values
for the global parameters $\nu$ and $\alpha$ (true values: $5$ and $1.0$,
respectively). For each point $s\in\S$, the percentage of results for which
the LEC algorithm has a better goodness of fit than EDC algorithm on the
region the point is in is depicted.

As a second example for a non-stationary process we use a process with
parameters $a_{s}=1$ constant, $b_{s}=3$ constant and
$g_{s}=(-x+5)\cdot\pi/2$. The dependence structure around each point is
locally an ellipse with the angle of the major axis rotating clockwise with
increasing $x$ (vertical for $x=-5$ , horizontal for $x=0$, vertical again for
$x=5$). In Fig. 6a and Fig. 6b we depict for one simulation the clusters
produced by the two algorithms, with the colors of the clusters indicating the
value of the estimate of parameter $g$ of the fitted processes. We depict the
true values of parameter $g_{s}$ in Fig. 6c. Again, the clusters of the EDC
algorithm group points with a high interdependency together, so their shape
follows the rotation of the ellipses (this is visible especially well in the
cluster in the middle of Fig. 6a and the two clusters above it), while the LEC
algorithm results in clusters that form vertical stripes and reconstruct the
spatial structure in the parameter values $g_{s}$,
$s\in\operatorname{{\mathcal{S}}}$ better. The analysis of the goodness of fit
yields similar results as for the first algorithm, with the LEC algorithm
being preferred in at least $75\%$ of the cases for every point.

As a third example, we use values of $a_{s}=(7.5-\|s\|)/2+1$, $b_{s}=0$,
$g_{s}=0$. This time, the true spatial structure is a bit different, it does
not feature vertical stripes, but instead a circular structure with the values
of $g_{s}$ depending on the distance of $s$ to the center. This circular
structure is visible in the clusters of the LEC algorithm (Fig. 7b). It is not
reproduced by the EDC algorithm, which results in clusters that are quite
uninformative this time, as the fitted values for $a$ are very similar for
each cluster (Fig. 7a). The results for the goodness of fit are similar to
those for the two examples before.

Figure 6: Results of the EDC algorithm (a) and the LEC algorithm (b) applied
to simulated data of a Huser-Genton process with parameters $a_{s}=1$
constant, $b_{s}=3$ constant and $g_{s}=(-x+5)\cdot\pi/2$ with $250$
observations. The number of clusters used is five. The colors within the
clusters show the estimated value for parameter $g$ on the stationary-max-
stable process that has been fitted to the data on the cluster. For reference,
the true values of the parameter $g_{s}$ are shown in (c). Figure 7: Results
of the EDC algorithm (a) and the LEC algorithm (b) applied to simulated data
of a Huser-Genton process with parameters $a_{s}=(7.5-\|s\|)/2+1$, $b_{s}=0$,
$g_{s}=0$ with $250$ observations. The number of clusters used is five. The
colors within the clusters show the estimated value for parameter $a$ on the
stationary-max-stable process that has been fitted to the data on the cluster.
For reference, the true values of the parameter $a_{s}$ are shown in (c).

##  5 Conclusions and Discussion

We have discussed a clustering algorithm by [7] that is used in multivariate
extreme value theory to group an area of investigation into smaller regions.
The obtained regions are then used to fit parametric stationary max-stable
processes to the data. This allows the application of such processes also to
data for which stationarity on the whole area cannot be assumed. While
clustering is in general a sensible and valid approach to the problem, the
dissimilarity measure used previously is not necessarily suitable to find
regions in which stationarity can be assumed, as we have discussed for one
concrete example of a non-stationary max-stable process. We propose a
different dissimilarity measure based on local estimates and demonstrate in a
simulation study that for three different examples of non-stationary data we
used, it indeed reconstructs the dependency structures of the data more
accurately and that the processes that were fitted to the data on the clusters
mostly have a better goodness of fit. For more general results regarding the
performance of the two methods, a more formal and mathematically more strict
analysis is required.

The clusters derived using the method by [7] do have a useful and meaningful
interpretation in another context: They show regions in which there is a high
extremal dependency between pairs of points. Therefore, a large cluster in a
certain area indicates a tendency for more large-scale extreme events in that
region. Finding such clusters is of relevance for example in the context of
insurances or risk management.

The clustering approach in this work and the subsequent fitting of max-stable
processes requires choosing two global parameters, and so far, we have not
found a systematic method to do so. For exactly two different choices of the
values, the clusters could be calculated and the goodness of fit of the
processes could be compared using the measure from Section 3.4. Unfortunately,
this measure does not allow for the comparison of more than two clusterings at
the same time. Carrying out a lot of pairwise comparisons is time-consuming
and does not seem to be a very convenient approach. Besides, there is no
guarantee that the results of these pairwise comparisons do not contradict
each other. A more sophisticated approach is definitely desirable here.

A more general limitation with the approach of using cluster algorithms and
then fitting regional max-stable models is that while those models can be used
to describe the data within one cluster and also give meaningful information
about how dependence structures vary spatially, they do not enable us to model
directly the dependency between two points that are in different clusters. In
this regard, the clustering approach is inferior to non-stationary max-stable
processes like the ones presented in [6] (which are on the other hand more
difficult to apply and require the availability of suitable covariates). It is
an interesting future research direction to use the max-stable processes that
have been fitted on the clusters and try to find a way to combine and extend
them to a parsimonious process covering the whole area of investigation.

Another possible application of the clustering algorithms is the spatio-
temporal investigation of extremes. In order to model changes in climate
extremes over time at one specific location, a GEV distribution with time-
dependent parameters can be used [18]. After fitting time-dependent models for
data at different locations, the data can then be transformed to have time-
stationary unit Fréchet distributions. By applying clustering algorithms to
the transformed data at different time windows, it is possible to combine the
temporal investigation of changes in extremes with an investigation of changes
in the spatial dependence structure.

## Statements and Declarations

  * •

Funding: Justus Contzen is funded through the Helmholtz School for Marine Data
Science (https://www.mardata.de; grant no. HIDSS-0005). Gerrit Lohmann
receives funding through "Ocean and Cryosphere under climate change" in the
Program "Changing Earth – Sustaining our Future" of the Helmholtz Society
(https://www.helmholtz.de/en/about-us/structure-and-governance/program-
oriented-funding/) and through PalMod by the Bundesministerium für Bildung und
Forschung (https://www.palmod.de/; grant no. 01LP1917A). The funders had no
role in study design, data collection and analysis, decision to publish, or
preparation of the manuscript.

  * •

Conflict of interest: The authors declare that they have no conflict of
interest.

  * •

Authors’ contributions: All authors contributed to the study conception and
design. The computer simulations were performed by Justus Contzen, who also
wrote the first draft of the manuscript. Gerrit Lohmann and Thorsten Dickhaus
supervised the study and commented on previous versions of the manuscript. All
authors contributed to reviewing and editing the manuscript. The final
manuscript was read and approved by all authors.

## References

  *   * Schlather [2002] Schlather, M.: Models for stationary max-stable random fields. Extremes 5, 33–44 (2002) https://doi.org/10.1023/A:1020977924878
  * Brown and Resnick [1977] Brown, B.M., Resnick, S.I.: Extreme values of independent stochastic processes. J. of Applied Probability 14(4), 732–739 (1977) https://doi.org/10.2307/3213346
  * Kabluchko et al. [2009] Kabluchko, Z., Schlather, M., de Haan, L.: Stationary max-stable fields associated to negative definite functions. Annals of Probability 37(5), 2042–2065 (2009) https://doi.org/10.1214/09-AOP455
  * Opitz [2013] Opitz, T.: Extremal t-processes: Elliptical domain of attraction and a spectral representation. J. of Multivariate Analysis 122, 409–413 (2013) https://doi.org/10.1016/j.jmva.2013.08.008
  * Ribatet [2017] Ribatet, M.: Modelling spatial extremes using max-stable processes. In: Franzke, C.L.E., O’Kane, T.J. (eds.) Nonlinear and Stochastic Climate Dynamics, pp. 369–391. Cambridge University Press, Cambridge (2017). https://doi.org/10.1017/9781316339251.014
  * Huser and Genton [2016] Huser, R., Genton, M.G.: Non-stationary dependence structures for spatial extremes. J. of Agricultural, Biological, and Environmental Statistics 21, 470–491 (2016) https://doi.org/10.1007/s13253-016-0247-4
  * Saunders et al. [2021] Saunders, K., Stephenson, A.G., Karoly, D.J.: A regionalisation approach for rainfall based on extremal dependence. Extremes 24, 215–240 (2021) https://doi.org/10.1007/s10687-020-00395-y
  * Bernard et al. [2013] Bernard, E., Naveau, P., Vrac, M., Mestre, O.: Clustering of maxima: Spatial dependencies among heavy rainfall in france. J. of Climate 26(20), 7929–7937 (2013) https://doi.org/10.1175/JCLI-D-12-00836.1
  * de Haan and Ferreira [2006] de Haan, L., Ferreira, A.: Extreme Value Theory: An Introduction. Springer, New York (2006)
  * McNeil et al. [2015] McNeil, A.J., Frey, R., Embrechts, P.: Quantitative Risk Management: Concepts, Techniques and Tools. Revised Edition. Economics Books. Princeton University Press, Princeton (2015)
  * de Haan [1984] de Haan, L.: A spectral representation for max-stable processes. Annals of Probability 12(4), 1194–1204 (1984) https://doi.org/10.1214/aop/1176993148
  * Penrose [1992] Penrose, M.D.: Semi-min-stable processes. Annals of Probability 20(3), 1450–1463 (1992) https://doi.org/10.1214/aop/1176989700
  * Davis et al. [2013] Davis, R.A., Klüppelberg, C., Steinkohl, C.: Max-stable processes for modeling extremes observed in space and time. J. of the Korean Statistical Society 42(3), 399–414 (2013) https://doi.org/10.1016/j.jkss.2013.01.002
  * Ribatet et al. [2015] Ribatet, M., Dombry, C., Oesting, M.: Spatial extremes and max-stable processes. In: Dey, D., Yan, J. (eds.) Extreme Value Modeling and Risk Analysis: Methods and Applications, pp. 179–194. Chapman and Hall/CRC, New York, ??? (2015). https://doi.org/10.1201/b19721-10
  * Cooley et al. [2006] Cooley, D., Naveau, P., Poncet, P.: Variograms for spatial max-stable random fields. In: Bertail, P., Soulier, P., Doukhan, P. (eds.) Dependence in Probability and Statistics. Springer, New York (2006)
  * Paciorek and Schervish [2006] Paciorek, C.J., Schervish, M.J.: Spatial modelling using a new class of nonstationary covariance functions. Environmetrics 17(5), 483–506 (2006) https://doi.org/10.1002/env.785
  * Wand and Jones [1994] Wand, M.P., Jones, M.C.: Kernel Smoothing (1st Ed.). Chapman and Hall/CRC, New York (1994). https://doi.org/10.1201/b14876
  * Contzen et al. [2023] Contzen, J., Dickhaus, T., Lohmann, G.: Long-term temporal evolution of extreme temperature in a warming earth. PLOS ONE 18(2), 1–30 (2023) https://doi.org/10.1371/journal.pone.0280503