# Comparative statistical study of two local clustering coefficient
formulations as tropical cyclone markers for climate networks

Mikhail Krivonosov  Lobachevsky State University of Nizhny Novgorod, Nizhny
Novgorod, Russia  Olga Vershinina  Lobachevsky State University of Nizhny
Novgorod, Nizhny Novgorod, Russia  Anna Pirova  Lobachevsky State University
of Nizhny Novgorod, Nizhny Novgorod, Russia  Shraddha Gupta  Potsdam Institute
for Climate Impact Research, PO Box 601203, 14412 Potsdam, Germany  Department
of Physics, Humboldt University, Berlin, Germany  Oleg Kanakov  Lobachevsky
State University of Nizhny Novgorod, Nizhny Novgorod, Russia  Jürgen Kurths
Lobachevsky State University of Nizhny Novgorod, Nizhny Novgorod, Russia
Potsdam Institute for Climate Impact Research, PO Box 601203, 14412 Potsdam,
Germany  Department of Physics, Humboldt University, Berlin, Germany

###### Abstract

We introduce a new formulation of local clustering coefficient for weighted
correlation networks. This new formulation is based upon a definition
introduced previously in the neuroscience context and aimed at compensating
for spurious correlations caused by indirect interactions. We modify this
definition further by replacing Pearson’s pairwise correlation coefficients
and three-way partial correlation coefficients by the respective Kendall’s
rank correlations. This reduces statistical sample size requirements to
compute the correlations, which translates into the possibility of using
shorter time windows and hence into a shorter response time of the real-time
climate network analysis. We construct evolving climate networks of mean sea
level pressure fluctuations and analyze anomalies of local clustering
coefficient in these networks. We develop a broadly applicable statistical
methodology to study association between spatially inhomogeneous georeferenced
multivariate time series and binary-valued spatiotemporal data (or other data
reducible to this representation) and use it to compare the newly proposed
formulation of local clustering coefficient (for weighted correlation
networks) to the conventional one (for unweighted graphs) in terms of the
association of these measures in climate networks to tropical cyclones. Thus
we substantiate the previously made observation that tropical cyclones are
associated with anomalously high values of local clustering coefficient, and
confirm that the new formulation shows a stronger association.

##  I Introduction

Network analysis is recognized nowadays as a powerful tool in complex system
studies, including climate science. In essence, it is a specialization of the
general correlation networks approach Langfelder and Horvath (2008), which is
also well established in systems biology Bruggeman and Westerhoff (2007); Van
der Greef  _et al._ (2007); Friedman and Alm (2012); Rosato  _et al._ (2018)
and neuroscience Fransson and Marrelec (2008); Bassett and Sporns (2017).
Correlation networks in climate science are constructed from cross-
correlations in multivariate time series of chosen climate variables between
nodes of a spatial grid (e.g. a geographic coordinate grid on Earth surface).
Such networks are then used to compute various graph metrics, which in turn
become subject to subsequent analysis aimed at identifying patterns of climate
dynamics Tsonis  _et al._ (2006); Donges  _et al._ (2009); Gozolchiani  _et
al._ (2011).

Tropical cyclones (TCs) are extreme weather events that can cause tremendous
destruction to life and property. Accurate forecasting and analysis of TCs and
their tracks remains a challenge and a major area of research in climate
science, even today. Climate networks, which have been widely used to study
different weather dynamics, were recently employed in the study of extreme
rainfall events associated with TCs Ozturk  _et al._ (2018, 2019); Traxl  _et
al._ (2016). In Gupta  _et al._ (2021), the authors have undertaken a time-
evolving functional network approach to capture the rapid spatio-temporal
evolution of the network topology of the regional weather system associated
with a TC. They propose suitable network metrics to detect TC tracks from mean
sea level pressure (MSLP) data and study their dynamics.

It is essential for the efficiency of this approach that as much as possible
of information contained in the initial data are retained in the feature set
produced by network analysis. Climate networks are typically (including those
studied in Gupta  _et al._ (2021)) unweighted graphs obtained by thresholding
the correlation matrix of the input data. Network nodes are associated to the
variables in the input data; in case of a single climatic observable, e.g.
MSLP, and, hence, in a single-layer network, the network nodes are in one-to-
one correspondence with the spatial grid nodes on the Earth surface. A link
between two specific nodes in the network is assumed to exist if the
corresponding component of the correlation matrix (i.e. the correlation
coefficient in this specific pair of variables) exceeds a chosen threshold.
The threshold value is a free parameter of the method, and generally allows
tuning to maximize performance (in whatever quantitative sense), but it is
also common to fix the threshold value by specifying the edge density (defined
as the fraction which linked pairs of nodes constitute among the total number
of node pairs), e.g. at 5% or 10%, without any further optimization.
Regardless of whether the threshold value is optimized or not, thresholding
inevitably leads to information loss.

A known way to mitigate the information loss introduced by thresholding is to
construct and analyze an ensemble of several unweighted graphs at once, by
using a series of threshold values; this approach has been successfully
applied to genetic network analysis Bockmayr  _et al._ (2013). It is also
possible to eliminate the thresholding operation at all, by analyzing a full
weighted graph whose edge weights are determined by (in the simplest case,
taken equal to) the respective pairwise correlation coefficients. The
transition from unweighted to weighted networks calls for a proper extension
of graph metrics definitions. For many metrics such extensions are available
Newman (2016), but may be not unique, thus giving rise to an additional
problem of choosing the best weighted-network modification of a particular
graph metric. An important consideration to guide this choice is that the
networks of interest are, in essence, correlation networks.

The weighted-network approach has been successfully implemented in the
neuroscience context Masuda  _et al._ (2018), where an improved formulation of
local clustering coefficient (LCC) specially focused on weighted correlation
networks was proposed and demonstrated to outperform the conventional
thresholding approach in terms of an illustrative neuroscience problem
(revealing the age dependence of human brain network structure based on
functional magnetic resonance data). The key distinctive feature of the
special definition of LCC for weighted correlation networks Masuda  _et al._
(2018) consists in accounting for three-way partial correlations in order to
compensate for spurious correlations caused by indirect interactions, which
otherwise disguise the true interaction structure in the correlation matrix.

No literature is currently available to implement this approach in climate
science. The present study aims at filling this gap by adapting the
correlations-focused weighted-network LCC definition from Masuda  _et al._
(2018) to climate network analysis. We further modify the LCC definition of
Masuda  _et al._ (2018) by switching from (parametric) Pearson’s to (non-
parametric) Kendall’s correlations (which moves from capturing linear to more
general monotonic dependencies between variables and imposes weaker
requirements to sample size as compared to both Pearson’s and Spearman’s
correlation Bonett and Wright (2000)) and compute this new LCC formulation for
MSLP correlation networks. We find and quantify the relation of thus computed
LCC to available data on the trajectories of tropical cyclones, thereby
supporting and extending the results of Gupta  _et al._ (2021), where LCC
anomalies of MSLP correlation networks were discovered in the vicinity of TCs.
We show that the LCC modification based on three-way partial correlations
outperforms the unweighted thresholded network LCC as a marker of tropical
cyclone. We expect this new measure to be especially effective for the
analysis of strongly time-dependent climate phenomena.

##  II Methods

###  II.1 Terminology, conventions and notation

Attribute will mean here any physical quantity which is measurable directly or
indirectly at a specific time and location (e.g. MSLP, see Sec. II.2), or a
computable quantity which is also attributed to a specific time and location
(e.g. LCC of an evolving correlation network, see Sec. II.5).

Observation is an individual value of an attribute (obtained by measurement or
computation).

Time grid is a set of moments in time at which observations are made. Here the
time grid $\\{t_{m}\\}$ is equidistant (periodic), and its period is referred
to as sampling period.

Spatial grid is a set of fixed locations in space at which observations are
made. We use a spatial grid of ERA5 single-level climate data Hersbach  _et
al._ (2020) (see details in Sec. II.2), which has no altitude dimension (MSLP
data are attributed to the mean sea level by definition), and the grid is
periodic in geographic coordinates (lattitude and longitude), but not
equidistant on the Earth surface. Since we use a single attribute
(specifically, MSLP) to construct correlation networks (see Sec. II.4), the
nodes of these networks are by construction in one-to-one correspondence to
the spatial grid nodes; such networks are referred to as single-layer
networks. The quantity of spatial grid nodes is denoted as $N$.

Space-time grid is the Cartesian product of the spatial grid and the time
grid.

Time series is a sequence of observations made at consecutive points of the
time grid. Univariate time series contain a single observation per each moment
of time, and multivariate consist here of simultaneous observations of a
single attribute (e.g. MSLP or LCC) at multiple locations (space grid nodes).
Generally, multivariate time series may as well encompass several attributes,
giving rise to multi-layer correlation networks.

Subscript indices $i$, $j$, $l$ will label the spatial grid nodes (and,
equivalently, the correlation network nodes); superscript bracketed indices
$(m)$, $(n)$ will label the time grid points. For example, the observation of
an attribute $v$ at the spatial node $i$ at the time $t_{m}$ will be denoted
as $v_{i}(t_{m})$, or, equivalently, as $v_{i}^{(m)}$; notations with braces
(e.g. $\\{v_{i}^{(m)}\\}$) mean the entire time series.

Variate generally means a random quantity, but here it will also refer to a
component of multivariate time series, when described in terms of probability
theory. A notation for such a variate (e.g. $v_{i}$) is obtained by omitting
the temporal index from the notation of the observations in the series.

Sample of a variate is generally a set of its observations (not necessarily
time-stamped or even ordered); here, a sample is always represented by a time
series. It is common in statistics to assume that the observations in a sample
are independent and identically distributed, but this is not the case within
our problem statement. Here we treat all time series under study as
realizations of discrete-time stochastic processes (Papoulis and Pillai, 2002,
Chap. 9). Stationarity and ergodicity of these underlying processes are
important issues to consider in order to treat a time series as a statistical
sample (see Sec. II.6, II.7).

Georeferenced time series means multivariate time series, where each component
(variate) is attributed to a specific location in space, which does not change
over time; these locations are nodes of the spatial grid.

Trajectory is another form of spatiotemporal data, where the location of a
moving object (here, TC center) is recorded over a time grid; these locations
are not assumed to coincide with spatial grid nodes.

Our study involves several statistical and graph-theoretic methods, which are
applied in sequence. Further, we provide a general description for each
method, and specify how the method is applied within our computation chain.
General input data in a method’s description mean the type of data to which
the method is generically applicable; actual input data mean the data to which
the method is actually applied in our computation.

We study association (dependence) between processed georeferenced time series
(here, that of LCC in two different formulations) on the one hand, and TCs
(represented by a set of available TC center trajectories) on the other. For
generality, we refer to a particular attribute under this association study
(here, to each particular formulation of LCC) as indicator, and to a class of
events (here, TCs) which occur in space and time and may be represented as
time series of a binary valued (taking on values ‘yes’ or ‘no’, see Sec. II.7)
attribute over the space-time grid as event of interest. Moreover, prior to a
statistical analysis for association, the indicator time series are
transformed (regularized) so as to provide that all variates within the
transformed series (each variate pertaining to a specific node of the spatial
grid) are characterized by the same universal probability distribution
(namely, uniform distribution on $[0,1]$, see Sec. II.6). The outcome of this
regularization will be referred to as transformed indicator. The association
analysis methodology described in Sec. II.6–II.7 may be applied to any
indicator and event of interest, provided they meet the assumptions outlined
in the respective methods descriptions.

###  II.2 Cimate data

As a source of initial data (in the form of georeferenced multivariate time
series) for constructing climate correlation networks, we use the ERA5
reanalysis data Hersbach  _et al._ (2020), which is essentially a model-based
interpolation of available climatic observational data to a regular grid over
geographic coordinates and time. The attribute under study is mean sea level
pressure (MSLP), i.e. the calculated value of local atmospheric pressure at
the standard (long-term average) local sea level. We use time series of MSLP
on a coordinate grid with step of 0.75° over latitude and longitude, and time
grid with sampling period of 3 hours. The data are taken for the sea surface
(land excluded) within the range 4.5°N to 30.75°N latitude and 49.5°E to
100.5°E longitude (northern part of the Indian Ocean), over the time range
00:00:00 01.01.1982 to 21:00:00 31.12.2020, thus the time grid consists of
$M_{\text{source}}$=113960 points.

As long as correlation networks are meant to characterize fluctuations of the
variates under study, it is beneficial to suppress the deterministic component
in the time series prior to the correlation analysis. For that sake, following
the approach of Gupta  _et al._ (2021), we preprocess the original ERA5 data
by subtracting the respective local daily climate normals for each location
(spatial grid node) and each date of year, which are computed by averaging the
attribute (here, MSLP) at the specific spatial grid node over all observations
for the particular date of year and over all years of observation (for
simplicity, data for 29 February of leap years are discarded). This eliminates
the annual cycle from each component of the time series. The result of this
subtraction, further referred to as MSLP anomaly and denoted
$\\{v_{i}^{(m)}\\}$, constitutes georeferenced time series on the same space
and time grids as the original ERA5 data. The MSLP anomaly data are subject to
the subsequent correlation analysis (see Sec. II.3, II.4).

The other type of spatiotemporal data in our study are the Best Track data
Mohapatra  _et al._ (2012) on tropical cyclones in the region, which are
produced and published by the Regional Specialized Meteorological Centre for
Tropical Cyclones Over North Indian Ocean (India) and contain information on
the position and strength of all registered tropical cyclones in the region
over time. Currently, we ignore the TC strength information and use only TC
center trajectories (see Sec. II.7).

###  II.3 Kendall’s correlation

To quantify correlation, in line with Gupta  _et al._ (2021), we use the
Kendall’s rank correlation coefficient $\tau$ Kendall (1938) (more precisely,
the Goodman-Kruskal modification thereof, see below). General input data for
computing a correlation coefficient of two random variates is their joint
bivariate sample. We compute pairwise correlation coefficients $\tau_{ij}$ for
the respective pairs of MSLP anomalies $v_{i}$ and $v_{j}$ (pertaining to the
spatial grid nodes $i$ and $j$). The actual input data for obtaining a
coefficient $\tau_{ij}$ is a joint sample of the variates $v_{i}$ and $v_{j}$
constructed from a section of the multivariate time series of MSLP anomalies
taken within a sliding time window (see Sec. II.4).

Given a bivariate sample consisting of $w$ joint observations
$\\{(v_{i}^{(k)},v_{j}^{(k)})\\}_{k=1..w}$, the Kendall’s rank correlation
coefficient $\tau_{ij}$ is defined as a normalized difference between the
counts of concordant and discordant pairs, i.e. pairwise combinations of
bivariate observations from the sample. Namely, a pair of bivariate
observations $\big{(}(v_{i}^{(m)}$, $v_{j}^{(m)})$, $(v_{i}^{(n)}$,
$v_{j}^{(n)})\big{)}$ is called concordant (discordant), if the product
$(v_{i}^{(m)}-v_{i}^{(n)})\cdot(v_{j}^{(m)}-v_{j}^{(n)})$ is positive
(negative), implying that these particular observations show an increasing
(decreasing) dependence between $v_{i}$ and $v_{j}$. Normalization ensures
that $\tau_{ij}\in[-1,1]$; in particular, $\tau_{ij}=1$ ($\tau_{ij}=-1$)
implies that all observations within the sample are compliant with some
deterministic increasing (decreasing) functional dependence between $v_{i}$
and $v_{j}$. Different formulations for the normalizing denominator are
available in the literature, depending on the chosen way to account for ties
(i.e. equality cases with $v_{i}^{(m)}=v_{i}^{(n)}$ or
$v_{j}^{(m)}=v_{j}^{(n)}$), which is still a matter of research Amerise and
Tarsitano (2015). As long as ties in our problem have a negligible impact due
to the continuous nature of the climate variables under study, the method of
their resolution is not actually important. For definiteness, we discard tied
pairs if such occur, and define the normalizing denominator as the total of
concordant and discordant pairs (ties excluded), as suggested in Adler (1957),
(Kruskal, 1958, Eq. (15.2) and below therein), the resultant quantity also
known as the Goodman-Kruskal gamma coefficient (Gibbons and Chakraborti, 2003,
Eq. (2.39)), which is named after the authors who applied it to variates
taking on finite sets of values (Goodman and Kruskal, 1954, Eq. (21)).

Our choice of the Kendall’s rank correlation coefficient over other available
measures of association is due to the following considerations: (i) it imposes
weaker requirements on the sample size $w$ as compared to both Pearson’s and
Spearman’s correlation coefficients Bonett and Wright (2000), which translates
into a shorter length of the sliding time window (see Sec. II.4), and hence a
better time resolution; (ii) like any rank-order statistics (including
Spearman’s correlation), it is invariant to any monotone nonlinear
transformations of variables (Gibbons and Chakraborti, 2003, Sec. 5.5), which
is not the case with Pearson’s correlation being only invariant to linear
transformations; (iii) it is more robust to outliers (large-amplitude noise)
than both Pearson’s and Spearman’s correlation coefficients Xu  _et al._
(2013).

###  II.4 Evolving correlation networks

Here we describe a method of constructing a time sequence of correlation
networks, which is referred to as an evolving network. This is a modification
of the approach taken in Gupta  _et al._ (2021). General input data can be any
multivariate time series; the actual input data are the georeferenced time
series of MSLP anomaly $\\{v_{i}^{(m)}\\}$ (see Sec. II.2).

The full symmetric matrix of Kendall’s pairwise correlation coefficients
$(\tau_{ij})$ is computed for input multivariate time series over a sliding
time window of fixed length $w$ points of the time grid, and attributed to the
ending time of the window. We take $w=16$; given the sampling period of 3
hours, the window length amounts to 2 days. We performed computations with
greater window sizes, but they produce weaker association of LCC anomalies to
TCs, which is explained by a worse time resolution, i.e. slower response time
of the correlation network to the processes associated with TCs.

Namely, for each point of the time grid of the input data, starting from the
$w$th and onwards (further referred to as the “current time”), we use a short
section (“window”) of the input time series consisting of $w$ consecutive
joint observations, of which the latest is at the current time, as a sample to
compute all pairwise correlation coefficients $\tau_{ij}$ between the variates
(as described in Sec. II.3).

We consider both weighted and unweighted correlation networks, the former
defined by the link weight matrix taken equal to the Kendall’s correlation
matrix $(\tau_{ij})$, and the latter constructed by thresholding the
correlation matrix: a pair of nodes are considered connected when their
correlation coefficient exceeds a threshold, which in turn is chosen so that
the fraction of connected nodes among the total number of node pairs (edge
density) is a specified quantity (taken equal to 5% and 10% in the
computations). Thus, the adjacency matrix $(a_{ij})$ of an unweighted network
(by definition, $a_{ij}=1$ if a link is present between the nodes $i$ and $j$,
and $a_{ij}=0$ otherwise) is expressed as

$a_{ij}=1(\tau_{ij}>\theta),$ (1)

where $\theta$ is the threshold value; the notation $1(\tau_{ij}>\theta)$ is
defined to be unity whenever the condition in the brackets holds, and zero
otherwise. As soon as the correlation matrices are symmetric, all correlation
networks under consideration are undirected graphs.

The procedure of constructing correlation networks is repeated with the
current time advanced sequentially to each further time grid point, and the
time window shifted accordingly 111We developed an optimized algorithm to
compute Kendall’s rank correlation coefficient over a sliding time window,
which benefits from reusing the computation results obtained at the previous
positions of the window.. This way, we obtain a sequence of networks, each
attributed to the respective current time, which can be seen as instances of a
time-evolving network. The time interval between consecutive instances of this
evolving network equals the sampling period of the input time series. The
chosen positioning of the current time within the time window (at the ending
point of the window) ensures causality of our analysis, which means that
computing any quantity attributed to the current time never requires the
knowledge of any data which are in the future relative to this current time.
The chosen time window length of $w=16$ observations (or 2 days) is a
reasonable trade-off between being too short for the Kendall’s coefficient to
be meaningful, and being too long to respond quickly enough to the TC-
associated changes in the climate system.

The resultant evolving networks are represented by a time sequence of weight
matrices $\\{(\tau_{ij})^{(m)}\\}$ for the weighted network, and adjacency
matrices $\\{(a_{ij})^{(m)}\\}$ for the unweighted network. Since the earliest
instance of a correlation network is attributed to the $w$th time grid point
of the source time series, the time grid of the evolving network lacks the
$w-1$ initial points of the source time grid; the total number of network
instances $M$ is accordingly smaller than the length of the source ERA5 time
series $M_{\text{source}}$:

$M=M_{\text{source}}-w+1.$ (2)

Since the time series of network measures (see Sec. II.5) obtained from an
evolving network inherit its time grid, this grid (of length $M$, starting at
the $w$th point of the source grid) is always implied below.

###  II.5 Local clustering coefficient

This Section describes the procedure of computing LCC in correlation networks,
which can be unweighted or weighted. General input data for defining and
computing LCC in an unweighted (weighted) undirected graph can consist of any
representation of the graph, e.g. the quantity of nodes $N$ and the adjacency
matrix $(a_{ij})$ (link weight matrix $(\tau_{ij})$). The actual input data
are the instances of the evolving correlation networks $\\{(a_{ij})^{(m)}\\}$,
$\\{(\tau_{ij})^{(m)}\\}$ obtained from the procedure described in Sec. II.4.
Further we omit the temporal index $(m)$, assuming that LCC computation is
performed in turn for each instance of an evolving network. Since the graphs
under study are essentially correlation networks, they admit the use of
specially defined LCC for weighted correlation networks Masuda  _et al._
(2018), as described below.

$A$ $B$

Figure 1: Sample network fragments illustrating LCC values of 1 and 0. Labeled
circles and solid lines shown in magenta are the nodes and the edges of the
subgraph which determines LCC of the node $i$. Green lines with dashed endings
denote insignificant connections to the rest of the network. Panel $A$: node
$i$ together with its direct neighbors constitute a clique, and thus LCC of
this node equals 1 (all neighbors of $i$ are also neighbors to each other).
Panel $B$: node $i$ with its neighbors constitute a star subgraph, and LCC of
$i$ equals 0 (none of the neighbors of $i$ are neighbors of each other).

The local clustering coefficient for a given node $i$ in a general unweighted
undirected graph is defined as the ratio of the number of closed triangles
(i.e. triplets of nodes with all three pairwise connections present)
containing the node $i$ to the total number of node triplets constructed from
the node $i$ and its neighbors (adjacent nodes) Watts and Strogatz (1998).
Equivalently, LCC is the fraction of directly connected (mutually adjacent)
pairs of nodes among all pairs composed of the direct neighbors of the node
$i$. According to Watts and Strogatz (1998), LCC characterizes the
“cliquishness” of the direct neighborhood of the node $i$, i.e. how close (in
terms of the quantity of mutual connections) this neighborhood is to being a
clique (which is a subgraph where all nodes are connected to each other). In
particular, LCC equals one for any node which actually belongs to a clique
composed of this node and all its neighbors (this is the case for the node $i$
in Fig. 1$A$); conversely, if a node with all its neighbors constitute a star
subgraph, then LCC of the star center is zero (this holds for the node $i$ in
Fig. 1$B$). Another instructive example is the Gilbert’s formulation of Erdős-
Rényi random graphs Gilbert (1959); Erdős and Rényi (1959), where an edge
between any pair of nodes is assumed to exist with the same preassigned
probability $p$ independently on any other edges in the graph. In such random
graphs, expected LCC of any node equals $p$, which becomes a small value in
large networks with sparse connections; hence, the global clustering
coefficient defined as average LCC over all nodes is also small, see
discussion in Watts and Strogatz (1998), (Newman, 2003, Sec. 4.1). In network
science LCC along with the global clustering coefficient are commonly accepted
as markers of so-called small-world networks, which are characterized, given
the total number of nodes and edges, simultaneously by large values of
clustering coefficients (as compared to Erdős-Rényi random graphs) and small
values of mean shortest path length (as compared to regular lattices); these
features have been actually found in a variety of real networks Watts and
Strogatz (1998); Newman (2003).

Following the definition above, LCC may be expressed in terms of the adjacency
matrix $(a_{ij})$ of an unweighted graph as

$LCC_{i}^{\text{unw}}=\frac{\sum_{1\leq j<l\leq N(j,l\neq i)}a_{ij}\cdot
a_{il}\cdot a_{jl}}{\sum_{1\leq j<l\leq N(j,l\neq i)}a_{ij}\cdot a_{il}}.$ (3)

Note that the denominator in (3) is essentially the number of unordered pairs
(2-combinations) among the nodes adjacent to $i$ and can be equivalently
expressed as $k_{i}(k_{i}-1)/2$, where $k_{i}$ is the degree of the node $i$
(the number of the nodes adjacent to $i$). We do not use this replacement in
the expression (3) in order to facilitate establishing its correspondence to
the weighted-network modifications of LCC formulation.

Multiple extensions of the LCC definition to account for link weights are
available in the literature Masuda  _et al._ (2018); Barrat  _et al._ (2004);
Onnela  _et al._ (2005); Zhang and Horvath (2005); Holme  _et al._ (2007);
Opsahl and Panzarasa (2009). Among them, the most straightforward modification
of Eq. (3) is the definition suggested in Zhang and Horvath (2005), which, as
shown in Kalna and Higham (2006), can be obtained from Eq. (3) simply by
replacing all adjacency matrix elements $a_{ij}$ with the respective
normalized link weights. Similar expressions were also derived in Grindrod
(2002); Ahnert  _et al._ (2007) based upon unweighted random graphs with
probabilistic links, where the probability for a link to exist between a pair
of nodes in a random graph is taken proportional to the respective link weight
from the underlying deterministic weighted network. Comparative reviews of
different formulations of LCC for weighted networks can be found in Saramäki
_et al._ (2007); Wang  _et al._ (2017).

The LCC formulation introduced in Masuda  _et al._ (2018) is specifically
tailored to correlation networks, which sets it apart from other LCC
definitions adapted to weighted graphs in general. The approach of Masuda  _et
al._ (2018) implies that a correlation matrix is an empirical manifestation of
some underlying interactions, and seeks to provide an empirical estimate of
LCC adequate to the network of that underlying interactions. The idea behind
the LCC definition in Masuda  _et al._ (2018) is to compensate for the impact
of indirect interaction paths (e.g. node $j$ interacting with node $i$, which
in turn interacts with node $l$, as in both panels of Fig. 1) upon the
correlation matrix. Namely, the nodes $j$ and $l$ in the example above are
expected to exhibit some level of correlation not only in the presence of an
actual direct interaction between them (e.g. when the actual interactions are
as in Fig. 1$A$), but also even in the absence of such interaction (as in Fig.
1$B$), due to indirect interaction paths (in particular, but not limited to,
the path mediated by the node $i$). Consequently, such indirect correlations
make the true underlying interaction structure hard to identify by the
correlation matrix, and moreover, this may dramatically affect the empirical
LCC value; e.g. the actual interaction structure shown in Fig. 1$B$ may
produce, in the worst case scenario, a correlation network with as many
spurious connections as in Fig. 1$A$, which translates into an inflation of
LCC from 0 in the actual underlying interactions network to 1 in the empirical
correlation network. In order to reduce the influence of such indirect
correlations on the resultant LCC value, the definition of LCC in Masuda  _et
al._ (2018) makes use of the three-way partial Pearson’s product-moment
correlation coefficient $\rho^{\text{part}}_{jl|i}$ Whittaker and Whittaker
(1990), which essentially indicates the surplus correlation between the nodes
$j$ and $l$ beyond their indirect correlation through the node $i$.

In our study we modify the LCC definition from Masuda  _et al._ (2018) by
replacing all Pearson’s product-moment correlation coefficients by the
respective Kendall’s rank correlation coefficients, which produces the
expression

$LCC_{i}^{\text{wei}}=\frac{\sum_{1\leq j<l\leq N,(j,l\neq
i)}|\tau_{ij}\cdot\tau_{il}\cdot\tau^{\text{part}}_{jl|i}|}{\sum_{1\leq
j<l\leq N(j,l\neq i)}|\tau_{ij}\cdot\tau_{il}|}$ (4)

differing from the initial unweighted definition (3) in the following: (i) the
adjacency matrix elements $a_{ij}$ and $a_{il}$, which denote the connections
between the node $i$ and its neighbors, are replaced by the respective
Kendall’s correlation coefficients $\tau_{ij}$ and $\tau_{il}$; (ii) the
adjacency matrix element $a_{jl}$, denoting the mutual interaction between the
node $i$’s neighbors $j$ and $l$, is replaced by the three-way partial
Kendall’s rank correlation coefficient $\tau^{\text{part}}_{jl|i}$, which in
turn is defined according to (Gibbons and Chakraborti, 2003, Section 12.6) as

$\tau^{\text{part}}_{jl|i}=\frac{\tau_{jl}-\tau_{ij}\cdot\tau_{il}}{\sqrt{1-\tau_{ij}^{2}}\sqrt{1-\tau_{il}^{2}}};$
(5)

(iii) all correlation coefficients, including the partial one, are taken by
absolute value, which implies that positive and negative correlations are
treated equally. The latter is in line with the approach taken in Masuda  _et
al._ (2018), which is substantiated by the reasoning that mutual correllations
and partial correlations, regardless of their sign, are equally interpreted as
indicators of interaction between the respective nodes. Notably, the
expression (5) is analogous to the definition of the three-way partial
Pearson’s correlation coefficient $\rho^{\text{part}}_{jl|i}$ Masuda  _et al._
(2018); Whittaker and Whittaker (1990) with all pairwise Pearson’s
correlations replaced by the respective Kendall’s correlations. Despite this
resemblance, the expression (5) for $\tau^{\text{part}}_{jl|i}$ is not
obtained merely by a formal substitution, but is properly justified based on
the general concept of partial correlations, see details in (Gibbons and
Chakraborti, 2003, Section 12.6).

As it is pointed out in Masuda  _et al._ (2018), the three-way partial
correlation coefficient is intended to compensate for only one specific
indirect interaction path between the nodes $j$ and $l$, namely, the path
mediated by the node $i$ (for which LCC is being computed), and ignores other
possible indirect interaction paths. This approach is well-suited to the
triangle-based structure of the LCC definition; that said, taking into account
every indirect interaction path could arguably further improve the adequacy of
the reconstructed network topology to the underlying interactions. However,
this task is difficult in terms of both computational burden and poor
tractability for statistical estimation Brier  _et al._ (2015), the latter
being especially important when the sample size (equal to the time window
length) is rather small due to the imposed limitations, as discussed in Sec.
II.4. These difficulties render the mentioned improvement unfeasible.

As long as the expression (5) contains a division by zero whenever
$|\tau_{ij}=1|$ or $|\tau_{il}=1|$, such values of the indices $j$ and $l$ are
excluded from summation both in the numerator and in the denominator of (4),
whenever this occurs (such cases turn out to be rare but not impossible in our
actual computations). Essentially, as mentioned in Sec. II.3, unity absolute
value of Kendall’s $\tau$ implies that the bivariate statistics within such a
pair of variables (between the nodes $i$ and $j$, or $i$ and $l$) is
indistinguishable by the data within the sample from a deterministic monotone
functional dependence. We use this observation to justify the exclusion of
such node pairs from summation, which may be interpreted as such a pair being
temporarily (for a particular position of the time window) lumped together
into a single node, as soon as they turn out to behave as a single node
anyway. That said, excluding unity-correlated pairs of nodes inevitably
introduces a methodical perturbation into the obtained LCC value, and the
final justification for this operation relies mostly on our observation that
such perturbations, being rare and relatively small, do not affect the
capability of thus computed LCC as a climate marker.

We also note that although three-way partial correlations are still not
expected to recover completely the true underlying interaction structure of
whatever network system, we expect this approach, following the reasoning of
Masuda  _et al._ (2018), to produce the most adequate and meaningful empirical
estimate of LCC from correlation analysis of multivariate data among all
available LCC formulations. In particular, it was found in Masuda  _et al._
(2018) experimentally, that LCC defined in terms of three-way partial
correlations better reveals the age dependence of brain connectivity by
functional magnetic resonance data, as compared to other available LCC
formulations for weighted networks Barrat  _et al._ (2004); Onnela  _et al._
(2005); Zhang and Horvath (2005), and to the LCC definition (3) for unweighted
(thresholded) correlation networks. The present study, to the best of our
knowledge, is the first to combine the three-way partial correlations approach
of Masuda  _et al._ (2018) with the Kendall’s formulation of correlation
coefficient in the definition of LCC for correlation networks (4), and to
apply this new LCC definition to climate networks.

As a local network measure, an LCC value is attributed to every network node.
For an evolving network, LCC of all nodes are computed for each network
instance, which can be represented by appending an additional temporal index
$(m)$ to all variables in Eqs. (3–5). The resultant LCC arrays thus constitute
georeferenced multivariate time series $\\{LCC_{i}^{\text{unw}}(t_{m})\\}$ and
$\\{LCC_{i}^{\text{wei}}(t_{m})\\}$ on the same spatial grid as the source
data, and with the time grid of length $M$ (see Eq. (2)) inherited from the
evolving network.

###  II.6 Empirical complementary cumulative distribution function transform

The nodes of the spatial grid are subject to unequal conditions due to the
natural spatial inhomogeneity of the climate system (including the influence
of land and all other objective reasons), and also due to the geometric non-
uniformity of the grid itself (which is periodic in geographic coordinates,
but not in actual space). This inhomogeneity in the source climate data is
partially compensated by subtracting the respective local daily climate
normals (see Sec. II.2). However, fluctuations and their inter-node
correlation properties inevitably retain some inhomogeneity, thus affecting
the constructed correlation networks and their graph-theoretic measures (here,
different LCC formulations). As soon as the georeferenced time series of these
measures (further referred to as indicators) are to be analysed for
association with events of interest (TCs), which are also spatiotemporal
phenomena, the indicator time series have to be transformed (regularized) such
that their spatial inhomogeneity (across the network nodes) is compensated as
far as possible prior to the association analysis. To this end, we use an
approach based on the empirical complementary cumulative distribution function
transform (a modification of rank transform), which we describe below in the
general form, and apply it to the georeferenced LCC time series.

It is common in mathematical statistics to address the problem of non-
identical, non-normal and a priori unknown distributions of variates by means
of the rank transformation Conover and Iman (1981), on which a variety of non-
parametric statistical methods are based Gibbons and Chakraborti (2003),
whereby the actual observations are replaced by their ranks, i.e. sequential
numbers within a sorted sample. The rank transformation can be normalized by
the sample size, so that the result fits into the $[0,1]$ range.

If the ranking order is ascending, the normalized rank transformation
coincides with the empirical cumulative distribution function (ECDF) (Gibbons
and Chakraborti, 2003, Sec. 2.3). Using ECDF transformation to regularize
input data is also a known practice in machine learning Zhuang  _et al._
(2020). ECDF of a variate $x$ is defined as the empirical frequency of
observing a value of $x$ which is same or smaller than the function argument
$\xi$ over the entire available sample (or time series) of the variate:

$ECDF_{x}(\xi)=\frac{1}{M}\sum_{m=1}^{M}1\left(x^{(m)}\leq\xi\right),$ (6)

where $M$ is the sample size (length of the time series), $x^{(m)}$ are the
observations in the sample, notation $1(x^{(m)}\leq\xi)$ has the same meaning
as in Eq. (1).

Under certain conditions (see below) $ECDF_{x}(\xi)$ is an unbiased consistent
empirical estimate of the underlying cumulative distribution function
$CDF_{x}(\xi)$, which means that $ECDF_{x}(\xi)$ has expectation equal to
$CDF_{x}(\xi)$, and converges to it (e.g. in the mean-square sense, meaning
that variance tends to zero) when the sample size increases (Gibbons and
Chakraborti, 2003, p. 18). In turn, CDF is known to transform any continuous
random variate so that the outcome is uniformly distributed on $[0,1]$;
precisely, the transformed variate $y=CDF_{x}(x)$ is distributed uniformly on
$[0,1]$, whenever $CDF_{x}(\cdot)$ is any continuous function (Gibbons and
Chakraborti, 2003, Sec. 2.5), (Embrechts and Hofert, 2013, Proposition 2).
This substantiates the practice of using ECDF (as an empirical estimate of
CDF) to compensate for non-identical and unknown distributions in multivariate
data by reducing all variates to a known universal distribution (namely,
uniform on $[0,1]$) in the limit of large sample size.

In turn, the descending ranking order may be considered similarly. In this
case, the normalized rank transformation is obtained by reversing the
inequality sign in the expression (6). This produces a function that we refer
to as the empirical complementary cumulative distribution function

$ECCDF_{x}(\xi)=\frac{1}{M}\sum_{m=1}^{M}1\left(x^{(m)}\geq\xi\right)$ (7)

with same notations as in Eq. (6). Note that $ECCDF_{x}(\xi)$ defined this way
does not actually complement $ECDF_{x}(\xi)$ to unity
($ECCDF_{x}(\xi)+ECDF_{x}(\xi)\not\equiv 1$); however, given a continuous
variate $x$, $ECCDF_{x}(\xi)$ provides an unbiased consistent estimate of the
underlying complementary cumulative distribution function

$CCDF_{x}(\xi)=1-CDF_{x}(\xi),$ (8)

and, similarly to a mentioned property of CDF, the transformed variate
$y=CCDF_{x}(x)$ is distributed uniformly on $[0,1]$ under the same
assumptions.

In order to choose between the ECDF and ECCDF transformations, we note that a
natural measure to quantify an observed anomaly of a georeferenced variable is
the empirical frequency of observing a same or stronger anomaly in the same
location over the entire observation time, so that the closer is this
frequency to zero, the greater the anomaly. When anomalies of interest are
small values (so that a stronger anomaly implies a smaller value), the
abovementioned frequency of same or stronger anomalies coincides with the
definition of ECDF (6). Conversely, when anomalies of interest are high
values, they are to be characterized by the empirical frequency of observing
same or greater values, which coincides with ECCDF (7). As soon as TCs turn
out to be associated with anomalously high values of LCC Gupta  _et al._
(2021), we use the ECCDF transformation (7). Due to the argument above, this
transformation not only compensates for the non-identical and apriori unknown
distributions within multivariate time series, but is also meaningful as the
frequency of same or stronger anomalies for the particular location.

The required property of ECDF (equivalently, ECCDF) to be an unbiased
consistent estimate of CDF (CCDF) is typically proven under the assumptions of
independence and identical distribution of the observations in the sample
(Gibbons and Chakraborti, 2003, Sec. 2.3). However, if the sample is
constructed from a time series, then typically one or both assumptions are not
fulfilled; in this case, the time series $\\{x^{(m)}\\}_{m=1..M}$ have to be
considered as a realization of an underlying stochastic process
$\\{\tilde{x}^{(m)}\\}$ (Papoulis and Pillai, 2002, Chap. 9), which is a
sequence of (generally, dependent) random variates numbered by the temporal
index $m$. For consistency with cited literature, in the reasoning below we
speak of ECDF and CDF; however, the same applies equivalently to ECCDF and
CCDF. Generally speaking, the variates $\tilde{x}^{(m)}$ may be differently
distributed; in this case there is no single distribution $CDF_{x}(\cdot)$
describing all observations. This problem is eliminated if the stochastic
process is first-order stationary (Papoulis and Pillai, 2002, p. 392) (i.e.
stationary with respect to distribution of an individual observation), meaning
that the distribution functions $CDF_{\tilde{x}^{(m)}}(\cdot)$ of all
individual variates $\tilde{x}^{(m)}$ are identical for all $m$, and thus they
admit a single notation $CDF_{x}(\cdot)$ without a temporal index (see
discussion below on the applicability of this assumption in the context of the
present study). If this holds, then the expectation of each term under
summation in Eq. (6) is

$\left<1\left(x^{(m)}\leq\xi\right)\right>=P\left(x^{(m)}\leq\xi\right)=CDF_{x}(\xi)$
(9)

for any $m$, and thus the expectation $\left<ECDF_{x}(\xi)\right>$ always
equals $CDF_{x}(\xi)$, or, in other words, $ECDF_{x}(\xi)$ is an unbiased
estimate of $CDF_{x}(\xi)$. In turn, if the underlying distributions
$CDF_{\tilde{x}^{(m)}}(\cdot)$ are not identical, then the expectation of
$ECDF_{x}(\xi)$ equals the arithmetic mean of these distribution functions
over the time grid.

Furthermore, the underlying process must be distribution-ergodic (Papoulis and
Pillai, 2002, p. 536) (or ergodic with respect to distribution of an
individual observation), which means that the estimate (6) converges to the
limit (e.g. in the mean-square sense) when the sample size $M$ increases. This
is equivalent to the definition of a consistent estimate (Gibbons and
Chakraborti, 2003, p. 19). A sufficient condition for a stationary process
$\\{\tilde{x}^{(m)}\\}$ to be distribution-ergodic (or, equivalently, for
$ECDF_{x}(\xi)$ to be a consistent estimate of $CDF_{x}(\xi)$) is the
assumption that a pair of individual variates $\tilde{x}^{(m)}$ and
$\tilde{x}^{(m^{\prime})}$ become independent, or their dependence vanishes
quickly enough, whenever they are separated by a sufficiently large time
interval, i.e. when $|m-m^{\prime}|$ is large enough (Papoulis and Pillai,
2002, p. 537). In general, the convergence of a time mean (such as in Eq.
(6))) to the respective underlying expectation is essentially a consequence of
ergodic theorems, which exist in a variety of formulations, e.g. pointwise
(Birkhoff’s) for almost sure convergence Birkhoff (1931), or mean (von
Neumann’s) for mean-square convergence von Neumann (1932); Birkhoff (1939).
Note that a stochastic process is not necessarily required to be stationary in
order to be ergodic. The ergodicity concept has been extended to a more
general class of asymptotically stationary processes, which includes periodic
nonstationarities (see Boyles and Gardner (1983); Anh and Lunney (1991) and
references therein), so that if an asymptotically stationary process is
ergodic, then an arithmetic mean over a finite realization (time series)
converges with increased realization length to a meaningful limit, namely, to
the pooled (time-averaged) expectation Bhagavan (1985), despite the
nonstationarity. To summarize, when the requirements of stationarity (or
asymptotic stationarity) and ergodicity of the underlying process are met, a
finite but sufficiently large time series $\\{x^{(m)}\\}_{m=1..M}$ can be used
as a sample of the variate $x$ in Eqs. (6) and (7), despite a possible mutual
dependence of observations within the series; however, if nonstationarity
(e.g. periodic one) is admitted, then its effect upon the result has to be
addressed in subsequent analysis.

Physically, the assumption of stationarity means that any changes in overall
conditions during the observation time are neglected, or that a possible
effect of such changes is not in the scope of the study. Climate data are
actually subject to nonstationarity, which includes daily and annual cycles,
and long-term trends of climate change. We do not expect daily cycle to
manifest in the LCC time series, as soon as the time window used to construct
evolving correlation networks (see Sec. II.4) is a multiple of 24 hours, and,
thus, always contains an integer number of daily periods. We also do not
expect the relatively small climate change during the available observation
time (which is almost 4 decades, see Sec. II.2) to impact significantly the
phenomena under study (i.e. MSLP correlation network indicators and their
association with TCs). At the same time, a possible annual nonstationarity in
the LCC time series may influence the results of our analysis, especially in
view of considerable seasonal variation in the frequency of tropical cyclones;
this issue is addressed in Sec. II.7.

In turn, ergodicity of time series generally arises from the underlying
dynamics, as described by the ergodic theory Cornfeld  _et al._ (2012);
Coudène (2016); Wiggins and Ottino (2004). Climate dynamics is commonly
assumed to be ergodic, when climate change can be neglected; for discussion of
the interplay between climate change and ergodicity, see Tél  _et al._ (2020).
As discussed above, periodic (including annual) nonstationarity does not
hinder ergodicity, and all other nonstationarities are neglected; thus, we
treat the time series under study as realizations of ergodic processes.

In order to speak of ECDF (ECCDF) as an approximation of the underlying CDF
(CCDF), it is also necessary to assume that the sample size is sufficient,
that is, the observation time interval covered by the available (finite) time
series encompasses a sufficient number of effectively independent observations
(despite their dependence within short time intervals). This is provided by
the observation time greatly exceeding both the annual period and the relevant
correlation time scales in the dynamics of the climate system.

We justify all the assumptions made above by comparing our final results to a
reference obtained from surrogate data (see Sec. II.7, Sec. III).

Note that using the rank transform or, equivalently, its normalized
modifications (ECDF or ECCDF) for regularizing differently distributed
variates in multivariate data implies that the transform is computed
individually for each variate (independently on other variates) Conover and
Iman (1981); Zhuang  _et al._ (2020). In our case this implies that ECCDF
transform is computed locally for each individual spatial grid node,
independently on other nodes, based on the entire available time series of the
indicator under study for this specific node (i.e. for the particular point on
the Earth surface).

The considerations above lead to the following procedure of multivariate time
series transformation (regularization). General input data are time series
denoted below as $\\{x_{i}^{(m)}\\}$. The actual input data are the
georeferenced time series of the two LCC formulations
$\\{LCC_{i}^{\text{unw}}(t_{m})\\}$ and $\\{LCC_{i}^{\text{wei}}(t_{m})\\}$
obtained from the procedure described in Sec. II.5; in actual computations,
they are substituted in turn for $\\{x_{i}^{(m)}\\}$ in the expressions below.

We consider each individual component of the multivariate time series
$\\{x_{i}^{(m)}\\}$ (at any given $i$) as a statistical sample of a random
variate $x_{i}$ in the sense which follows from the assumptions of (at least,
asymptotic) stationarity and ergodicity, as discussed above. Then, writing
down the ECCDF transform (7) for each variate $x_{i}$ instead of the generic
variate $x$, and substituting each individual observation $x_{i}^{(m)}$ in
turn for the function argument $\xi$, we express the transformation as

$f_{i}^{(m)}=ECCDF_{x_{i}}\left(x_{i}^{(m)}\right)=\\\
\frac{1}{M}\sum_{m^{\prime}=1}^{M}1\left(x_{i}^{(m^{\prime})}\geq
x_{i}^{(m)}\right),$ (10)

where the result of the transformation constitutes new multivariate time
series of the transformed indicator $\\{f_{i}^{(m)}\\}$.

Under the assumptions of stationarity and ergodicity of the underlying
stochastic process for any $i$, the transformation (10) converges (in the
mean-square sense) at large sample size $M$ to the complementary CDF

$f_{i}^{(m)}\xrightarrow[M\to\infty]{}1-CDF_{x_{i}}(x_{i}^{(m)}),$ (11)

so that each $f_{i}^{(m)}$ can be seen as an observation of a new
(transformed) variate $f_{i}=CCDF_{x_{i}}(x_{i})$, which is distributed
uniformly on $[0,1]$ for any $i$ by construction. This does not hold if the
stationarity assumption is replaced by asymptotic stationarity; in this case,
the time mean in Eq. (10) converges to the pooled (i.e. time-averaged Bhagavan
(1985)) CCDF

$f_{i}^{(m)}\xrightarrow[M\to\infty]{}1-\overline{CDF}_{x_{i}}(x_{i}^{(m)}),$
(12)

where

$\overline{CDF}_{x_{i}}(\xi)=\lim_{M\to\infty}\frac{1}{M}\sum_{m=1}^{M}CDF_{\tilde{x}_{i}^{(m)}}(\xi).$
(13)

However, if the nonstationarity is periodic and small (which is assumed), then
each individual distribution $CDF_{\tilde{x}_{i}^{(m)}}(\xi)$ weakly differs
from the pooled one, thus the distributions of the transformed variates remain
close to the uniform one, although having a small modulation in time with the
period of the nonstationarity.

Note an instructive case when all observations in the sample
$\\{x_{i}^{(m)}\\}_{i=const,m=1..M}$ are different (there are no ties, i.e. no
such $m$ and $n$ that $m\neq n$ and $x_{i}^{(m)}=x_{i}^{(n)}$); then, the
sample of any variate $f_{i}$, i.e. the set
$\\{f_{i}^{(m)}\\}_{i=const,m=1..M}$ computed according to Eq. (10), is
automatically (without any further assumptions) a permutation of the
equidistant grid $\\{k/M\\}_{k=1..M}$.

Furthermore, according to Eq. (10), each individual transformed observation
$f_{i}^{(m)}$ equals the frequency of observations $x_{i}^{(m^{\prime})}$
satisfying the condition $x_{i}^{(m^{\prime})}\geq x_{i}^{(m)}$ within the
sample $\\{x_{i}^{(m^{\prime})}\\}_{i=const,m^{\prime}=1..M}$, i.e. the
empirical frequency of observing a value of $x_{i}$ greater than or equal to
the respective observation $x_{i}^{(m)}$ at given $i$.

When the georeferenced time series of a chosen indicator
($\\{LCC_{i}^{\text{unw}}(t_{m})\\}$ or $\\{LCC_{i}^{\text{wei}}(t_{m})\\}$)
are substituted instead of $\\{x_{i}^{(m)}\\}$ in the ECCDF transformation
(10), any particular (“current”) observation within the transformed indicator
time series $f_{i}^{(m)}$ (at given $i$ and $m$) means the empirical frequency
of observing a same or stronger anomaly of the indicator in the particular
point on the Earth surface (at the node $i$) over the entire observation time,
as compared to the current (at the same $i$ and $m$) observation of this
indicator. Note that if the described methodology is applied to indicators
expected to demonstrate anomalously low values associated to the event of
interest, then, for the interpretation above to remain valid, the ECDF
transformation (6) should be used instead of ECCDF; this results in reversing
the inequality sign in Eq. (10).

###  II.7 Nonlocal association of georeferenced indicator time series to TC
trajectories

Our goal is to characterize and compare the extent to which the indicators
under study (here, different LCC formulations) demonstrate anomalous behavior
associated with events of interest (TCs). While indicators are represented by
georeferenced time series, TCs are represented by Best Track data Mohapatra
_et al._ (2012), which consist of TC center trajectories and TC strength
information over time. Since we currently ignore the TC strength information
and use only TC trajectories, the problem can be formulated in terms of
establishing association between two different types of spatio-temporal data:
georeferenced transformed indicator time series $\\{f_{i}^{(m)}\\}$ (obtained
from the procedure described in Sec. II.6) on the one hand, and a set of TC
trajectories on the other. Due to the absence of apriori information on the
spatial positioning of a TC-associated indicator anomaly relative to the TC
center (in particular, an indicator anomaly in general may even preferably
occur not at the TC center, but at a distance), the mentioned association
problem is nonlocal in space, in the sense that we speak of association
whenever an indicator anomaly and a TC are observed at the same time within a
certain (apriori unknown) allowable distance in space.

We quantify this association by defining a spatial vicinity of a TC center in
any moment of the time grid, introducing the null hypothesis of independence
between the indicator and TC trajectories, and posing the question whether
(and to what extent) the actual (empirical) statistics of indicator
observations taken within the joint vicinity of all TCs (i.e. conditionally
upon belonging to a TC center vicinity) over the entire time of observation
differs from the result expected under the null hypothesis. Namely, under this
null hypothesis and with stationarity assumed, the conditional distribution of
$f_{i}^{(m)}$ must coincide with the unconditional, which, in turn, is uniform
on $[0,1]$ by construction (see Sec. II.6). Therefore, the divergence of the
empirical conditional distribution of $f_{i}^{(m)}$ from the uniform one will
be interpreted as a measure of association (dependence) between the indicator
under study and the condition (of belonging to a TC center vicinity); at the
same time, the effect of possible annual nonstationarity of the indicator time
series has to be assessed. This reasoning is formalized and detailed below.

The size (and the shape) of the mentioned spatial vicinity is a free parameter
of the method. For the sake of simplicity, we define the vicinity of a TC at a
given time $t_{m}$ as a quadratic-shaped subarray of the spatial grid of a
specified size $R$, centered closest to the current position of the TC center.
This subarray is denoted as $v_{c}^{(m)}(R)$, where $c$ denotes a TC from the
available set. We tried subarray sizes ranging from $2\times 2$ grid points
($R=1$) up to $12\times 12$ points ($R=6$). Note that the spatial grid is
uniformly spaced in terms of geographic coordinates, but not in terms of their
actual distance on the Earth surface, because equal spacing in longitude
translates into different spacings in distance when taken at different
latitudes. The actual spacing of the grid points is proportional to the cosine
of the latitude, which in the region of interest (5°N to 30°N latitude)
produces the greatest difference of approximately 15%, which we neglect for
simplicity, although this effect may need to be taken into account in the
polar regions.

To each space-time grid point $(i,m)$ we ascribe a binary (taking on values
‘yes’ or ‘no’) attribute $V_{R,i}^{(m)}$, which answers the question whether
this grid point belongs to a vicinity of size $R$ of some TC center at the
respective time $t_{m}$, according to the following rule:

$V_{R,i}^{(m)}=\begin{cases}\text{yes},\quad\text{if }\exists c\text{ so that
}i\in v_{c}^{(m)}(R),\\\ \text{no},\quad\text{otherwise}.\end{cases}$ (14)

Consider a subset $S_{V}$ of the space-time grid $\\{(i,m)\\}$ satisfying the
condition $V_{R,i}^{(m)}=\text{yes}$:

$S_{V}=\\{(i,m)|V_{R,i}^{(m)}=\text{yes}\\},$ (15)

which is the joint vicinity of all TC trajectories. Then, consider a
conditional ECDF defined similarly to Eq. (6) over a sample of observations
$\\{f_{i}^{(m)}\\}_{(i,m)\in S_{V}}$, which are taken conditionally upon the
spacetime position of the observation belonging to $S_{V}$, i.e. upon
$V_{R,i}^{(m)}=\text{yes}$, as follows:

$ECDF_{f|S_{V}}(\xi)=\frac{1}{|S_{V}|}\sum_{(i,m)\in
S_{V}}1\left(f_{i}^{(m)}\leq\xi\right),$ (16)

where $|S_{V}|$ is the cardinality (number of elements) of the set $S_{V}$.
Following the lines of Sec. II.6, we consider the multivariate time series
$\\{f_{i}^{(m)}\\}$ and $\\{V_{R,i}^{(m)}\\}$ as realizations of underlying
multivariate stochastic processes $\\{\tilde{f}_{i}^{(m)}\\}$ and
$\\{\tilde{V}_{R,i}^{(m)}\\}$. As soon as the index $i$ maps into space (i.e.
on the Earth surface), these processes are essentially random fields, which
means a random function of point in spacetime (here, discrete).

If the process (field) $\\{\tilde{f}_{i}^{(m)}\\}$ is stationary (which
follows from $\\{\tilde{x}_{i}^{(m)}\\}$ being stationary), then all
individual variates $\tilde{f}_{i}^{(m)}=CDF_{x_{i}}(\tilde{x}_{i}^{(m)})$ are
identically distributed (namely, uniformly on $[0,1]$) in the limit of large
$M$ by construction (see Sec. II.6); hence (similarly to Eq. (9)), Eq. (16) is
an unbiased estimate of the underlying conditional CDF. Similarly to the
argument in Sec. II.6, the assumption of distribution-ergodicity of the
underlying random field implies that the arithmetic mean in Eq. (16) converges
(e.g. in the mean square) with increasing sample size $|S_{V}|$ to the
conditional expectation of the quantity being averaged, and finally reduces to
the conditional CDF, thus providing its consistent estimate:

$ECDF_{f|S_{V}}(\xi)\xrightarrow[|S_{V}|\to\infty]{}\left\langle
1\left(f_{i}^{(m)}\leq\xi\right)\big{|}(i,m)\in S_{V}\right\rangle\\\
=P\left(f_{i}^{(m)}\leq\xi\big{|}V_{R,i}^{(m)}=\text{yes}\right)\\\
=CDF_{f|V}\left(\xi\big{|}V_{R,i}^{(m)}=\text{yes}\right).$ (17)

Moreover, under the null hypothesis, the random fields
$\\{\tilde{f}_{i}^{(m)}\\}$ and $\\{\tilde{V}_{R,i}^{(m)}\\}$ are independent,
and, thus, the conditional distribution equals the unconditional one

$CDF_{f|V}\left(\xi\big{|}V_{R,i}^{(m)}=\text{yes}\right)=CDF_{f}\left(\xi\right),$
(18)

which is uniform on $[0,1]$ by construction under the assumptions taken.

A difference of the argument above as compared to the reasoning behind Eq.
(11) in Sec. II.6 is that the arithmetic mean in Eq. (16) is taken over a
subset $S_{V}$ of the (3-dimensional) space-time grid, in contrast to a mean
over a 1-dimensional time grid in Eq. (6). It implies that the distribution-
ergodicity condition found in (Papoulis and Pillai, 2002, p. 536), which was
used in Sec. II.6, can not be applied here directly, because the consideration
in Papoulis and Pillai (2002) is essentially limited to the case where
stochastic process realizations are averaged over a single dimension (time).
Applying the concept of ergodicity to justify convergence of an arithmetic
mean of a random field over a multi-dimensional space requires an accordingly
generalized ergodic theorem. Such a generalization was first obtained in
Wiener (1939) for a continuous space; alternative formulations and proofs, and
also further generalizations can be found e.g. in Pitt (1942); Becker (1981);
Anh and Lunney (1991); in particular, a generalization for discrete
multidimensional spaces (i.e. for random fields defined over a
multidimensional grid, which is the case in Eq. (16)) is available in (Pitt,
1942, Sec. 4), and for nonstationary (asymptotically stationary) fields (see
Sec. II.6) in Anh and Lunney (1991). Physically, a sufficient condition to
justify the limit transition in Eq. (17) may be formulated similarly to the
one-dimensional case (cf. Sec. II.6) as the requirement that a pair of
variates $\tilde{f}_{i}^{(m)}$ and $\tilde{f}_{i^{\prime}}^{(m^{\prime})}$
become independent, or their dependence decays quickly enough, when the
distance between the points of the space-time grid $(i,m)$ and
$(i^{\prime},m^{\prime})$ is sufficiently large.

Also similarly to the one-dimensional case (cf. Sec. II.6), the finite-sample
empirical conditional distribution (16) is expected to approximate well its
infinite-sample limit (17), if the set $S_{V}$ is large enough that the
(conditional) sample $\\{f_{i}^{(m)}\\}_{(i,m)\in S_{V}}$ encompasses a
sufficient number of effectively independent observations, despite their
dependence within short time intervals and short spatial distances. The
physical justification of this assumption relies on the consideration that the
set $S_{V}$, which is the joint vicinity of all TCs, covers a large number of
TCs scattered over the entire observation time, which greatly exceeds the
relevant correlation time scales of the climate dynamics.

In turn, in case of periodic nonstationarity, the distributions of the
variates $\tilde{f}_{i}^{(m)}$ vary in time in the vicinity of the uniform
distribution (see Sec. II.6) with the period of the nonstationarity. The
ergodic theorem for asymptotically stationary (which include periodic) random
fields (Anh and Lunney, 1991, Theorem 2) requires that the subset over which
the mean is taken (here, $S_{V}$) is rectangular and grows infinitely in all
dimensions. Physically, when the subset encompasses a sufficient number of
nonstationarity periods, the nonstationarity effectively averages out, thus
ensuring ergodicity (convergence of the mean). Note, however, that the actual
set $S_{V}$ in Eq. (16) does not satisfy this requirement, because essentially
it is a union of disjoint TC vicinities in spacetime, while the duration of
each TC is much shorter than the annual cycle, and the frequency of TC
occurrence is also subject to annual nonstationarity. Together with annual
nonstationarity of the indicator time series, this hinders the ergodic limit
transition in Eq. (17), making the outcome of the averaging in Eq. (16)
dependent on the properties of the set $S_{V}$. For example, as TCs in the
region of interest (Northern Indian Ocean) are known to occur most frequently
during the SOND season (September through December) Sugi  _et al._ (2002),
spacetime grid nodes falling within this season will be prevalently present in
$S_{V}$. This is not an issue with stationary indicator time series (because
then all moments of time are equivalent), but an annual nonstationarity (e.g.
if an indicator by itself is prone to anomalies within the SOND season,
regardless of TCs) makes the average over $S_{V}$ in Eq. (16) differ from the
global (unconditional) average even under the null hypothesis.

In order to assess the magnitude of the influence of annual nonstationarity
upon our results, we initially perform computations assuming the LCC time
series to be stationary, and then repeat the processing stages described in
Sec. II.6–II.7 on a subset of the LCC time series consisting only of
observations falling within the SOND season of each year, instead of the full
time series. The deviation between the results obtained from the SOND-only and
full time series are interpreted as a measure of the effect of annual
nonstationarity (see Sec. III).

In turn, a deviation of the empirical conditional distribution (16) from the
uniform one is an empirical evidence in support of rejecting the null
hypothesis, i.e. an evidence of association (dependence) between the indicator
and TC trajectories.

For the sake of visualization, the mentioned empirical distributions are
represented in the form of probability density plots, instead of cumulative
distributions (16). The empirical conditional probability density functions
$w_{f|S_{V}}(\xi)$ are obtained from the conditional sample
$\\{f_{i}^{(m)}\\}_{(i,m)\in S_{V}}$ using the kernel density estimation
method Rosenblatt (1956); Parzen (1962) implemented in the Scikit-Learn
software library Pedregosa  _et al._ (2011).

In order to ensure the correctness of all assumptions, we test the above
methodology against surrogate data, which are obtained by shifting the TC
trajectories from the actual TC Best Track data by a random amount along the
time grid. This shift simulates the fulfillment of the null hypothesis, i.e.
independence between indicators under study and TC tracks. The time series
$\\{{V}_{R,i}^{(m)}\\}$ (and, hence, the conditional sample
$\\{f_{i}^{(m)}\\}_{(i,m)\in S_{V}}$) are constructed again based on these
surrogate data according to Eqs. (14), (15), and conditional empirical
probability densities calculated based on this (surrogate) conditional sample,
all other data same as in the main computation.

##  III Results

We have processed the georeferenced time series of MSLP anomaly (see Sec.
II.2) by applying sequentially the procedures described in Sec. II.4–II.7, as
follows. The weighted evolving correlation network was constructed by
computing all pairwise Kendall’s correlation coefficients in a sliding time
window over the source time series, as described in Sec. II.4; the unweighted
evolving networks were obtained by thresholding the correlation coefficients,
so that the edge density (the proportion of connected node pairs among the
total combinatorial number of node pairs) is a specified value. We constructed
unweighted networks with edge density taken equal to 5% and 10%. Time series
of local clustering coefficients were calculated for the weighted and
unweighted evolving networks using the respective formulations in Sec. II.5.
To provide an assessment of the impact of annual nonstationarity in the LCC
time series upon the final result, additional LCC time series (further
referred to as SOND-only data) were prepared by keeping only the data for
September through December of each calendar year and omitting the rest. The
ECCDF transform, as given by Eq. (10), was applied to each of the available
LCC time series (to the full and SOND-only data, for the weighted and the
unweighted networks). Based on the transformed indicator time series and Best
Track data for the trajectories of tropical cyclone centers (see Sec. II.2),
empirical conditional probability density functions $w_{f|S_{V}}(\xi)$ (i.e.
empirical density of ECCDF-transformed indicator values observed within the
joint set of all TC vicinities) were computed as described in Sec. II.7.
Computations were performed with the sizes of the quadratic-shaped spatial
vicinity $v_{c}^{(m)}(R)$ varied from $2\times 2$ to $12\times 12$ spatial
grid nodes. Surrogate data for TC trajectories were obtained by shifting the
actual trajectories in time, so that independence between TCs and LCC time
series is simulated, and the same empirical conditional probability density
functions were calculated also with these surrogate data.

Recall that the transformed indicator value by itself means the unconditional
frequency of same or stronger anomalies of the indicator for the particular
location, and is distributed uniformly on $[0,1]$ with the unconditional
probability density equal to unity (see Sec. II.6). Hence, the value of the
conditional probability density $w_{f|S_{V}}(\xi)$ (e.g. the hight of a peak)
denotes how much more frequently indicator anomalies with given unconditional
probability $\xi$ are observed in the vicinity of a TC, than unconditionally.
We compare plots of the calculated empirical conditional distributions
visually, and interpret a greater deviation of the conditional distribution
from the uniform one as a stronger dependence between the particular indicator
and the events of interest (TCs).

Examples of the obtained empirical conditional probability density functions
are shown in Fig. 2. The results for the unweighted networks with edge density
equal to 5% and 10% do not differ much; the latter demonstrate slightly better
association with TCs and hence are shown in the plots. The strongest
association (i.e. the greatest deviation of the empirical conditional
distribution from the uniform one) for the weighted LCC is obtained for
spatial vicinity of size $4\times 4$ spatial grid nodes (panel (a) in Fig. 2),
and for the unweighted LCC for spatial vicinity of size $2\times 2$ nodes
(panel (b) in Fig. 2). The weighted LCC formulation shows stronger association
with TCs than the unweighted one, as quantified by an almost twice higher peak
in the conditional probability density, which means an almost twice more
frequent occurrence of the respective (having equal unconditional frequency)
anomalies in the TC vicinity for the weighted as compared to the unweighted
LCC formulation.

(a) (b)

Figure 2: Empirical conditional probability density plots constructed from
ECCDF-transformed indicator observations taken within the joint vicinity of
all TC trajectories. Abscissa axis: value of a transformed indicator, ordinate
axis: probability density. Indicators are local clustering coefficients for
the weighted (red lines) and unweighted (blue lines, edge density is 10%)
correlation networks. Vicinity size is $4\times 4$ spatial grid nodes in panel
(a) and $2\times 2$ in panel (b). Solid (dashed) lines denote the results
obtained from the full (SOND-only) data, and dash-dotted lines — from
surrogate data.

##  IV Conclusion

We have developed a broadly applicable statistical methodology, which is based
upon normalized rank transformation and conditional empirical distributions,
aimed at studying association (dependency) between an indicator represented by
georeferenced multivariate time series and a class of “events of interest”,
which may be represented as an ensemble of spatiotemporal trajectories, or in
any form reducible to a binary-valued (taking on values ‘yes’ or ‘no’) time
series over the same space-time grid. The indicator time series are
interpreted as a realization of a stochastic process (or random field), which
is assumed to be stationary and ergodic with respect to probability
distribution of an individual observation. Stationarity means that this
distribution does not depend on the observation time, and ergodicity is
provided when pairwise dependence between observations taken in different
points of spacetime decays quickly enough (i.e. the pairwise observation
distribution factorizes into a product of individual observation
distributions) with increased distance in spacetime between these
observations. In turn, the indicator probability distribution is not required
to be identical in different points across space; such difference in
distributions is compensated using normalized rank transformation prior to the
association analysis.

We consider two versions of the normalized rank transformation, which differ
only in the direction of the ranking order, and are equivalently applicable in
the context of the present study. We choose the type of the transformation so
that its result may be interpreted as the empirical frequency of observing
same or greater values of the indicator as compared to the current
observation, calculated for a particular location on the Earth surface. The
anomalies of interest (here, anomalously high values) are thus represented by
small values of the transformed indicator. This choice is not essential for a
single indicator, but in studies involving several indicators simultaneously
(e.g. as inputs of a multi-input classifier) it may be beneficial to choose
the type of the transformation for each indicator individually, depending upon
the particular sign of the anomaly (high or low value), which is of interest
for this indicator, so that the anomalies of interest are consistently
transformed into small values for all indicators.

Physically, the stationarity assumption implies the absence of varied external
conditions (or leaving the effects of such variation out of the scope of the
study), and ergodicity is provided by sufficiently strong chaoticity (such as
the mixing property) of the underlying dynamics. The sample size, i.e. the
quantity of effectively independent (despite a possible short-range
dependence) observations, must be sufficient in the sense that empirical
frequencies can be seen as a valid approximation of the underlying
probabilities.

If the underlying process is subject to periodic or quasiperiodic
nonstationarity (generally, if it is asymptotically stationary), then part of
the calculations remain meaningful, but the final result may be biased, thus,
the impact of a possible nonstationarity needs to be assessed. We do that by
repeating the calculation on a subset of the indicator time series obtained by
omitting a part of each time period of possible nonstationarity from the full
time series (here, we keep the data for September through December of each
calendar year and omit the rest), and comparing the result obtained from the
data with omissions to that from full data.

The applicability of other assumptions (including ergodicity of the underlying
process and sufficiency of sample size) is justified by testing the final
result against surrogate data simulating the null hypothesis (i.e.
independence between the indicator and the event of interest), which are
obtained by a time shift of the spatiotemporal data for the events of
interest.

We have introduced a new formulation of local clustering coefficient for
weighted correlation networks, which combines the known approaches based upon
three-way partial correlations and Kendall’s correlations. This modified
formulation is intended to reveal more information on the underlying
interaction structure as compared to the local clustering coefficient for
unweighted networks, which are obtained by thresholding the pairwise
correlation coefficients. Using three-way partial correlations aims at
mitigating the indirect interactions effect upon the correlation network,
while Kendall’s correlation admits a shorter length of the sliding time window
(and hence a better time resolution) when used in the context of evolving
correlation networks, which capture changes in the network topology over time.
This is especially important in the studies of processes occurring on time
scales which are only moderately (e.g. 10 to 20 times) greater than the
sampling period of the source time series used to construct a correlation
network.

We have applied the above statistical methodology to compare the two mentioned
local clustering coefficient formulations (calculated for an evolving
correlation network of mean sea level pressure) in terms of their association
to trajectories of tropical cyclone centers over the northern part of the
Indian Ocean. The statistical analysis (i) confirms the association of
anomalously high values of the local clustering coefficient to tropical
cyclones, and (ii) shows that the newly suggested measure demonstrates a
stronger association than the conventional one.

We expect the newly suggested local clustering coefficient formulation to be a
powerful indicator for predicting extreme climate events, including tropical
cyclones. The prediction problem is likely to benefit from a simultaneous use
of several indicators, including network measures and non-network climate
variables. This can be achieved e.g. by using learnable multi-input
classifiers and machine learning techniques.

Our statistical methodology of association analysis can be applied to a
variety of problems in climate studies, such as analyzing connections between
various climate network measures (and their combinations) and particular
climate events. It can also find application beyond the context of climate
science, in establishing association between a spatially inhomogeneous random
field and an ensemble of spatiotemporal trajectories or binary-valued
spatiotemporal data.

## Acknowledgement

This research was supported by the Russian Ministry of Science and Education
under Agreement No. 075-15-2020-808.

## References

  * Langfelder and Horvath (2008) P. Langfelder and S. Horvath, BMC bioinformatics 9, 1 (2008).
  * Bruggeman and Westerhoff (2007) F. J. Bruggeman and H. V. Westerhoff, TRENDS in Microbiology 15, 45 (2007).
  * Van der Greef  _et al._ (2007) J. Van der Greef, S. Martin, P. Juhasz, A. Adourian, T. Plasterer, E. Verheij, and R. McBurney, Journal of proteome research 6, 1540 (2007).
  * Friedman and Alm (2012) J. Friedman and E. J. Alm, PLoS Comput Biol 8, e1002687 (2012).
  * Rosato  _et al._ (2018) A. Rosato, L. Tenori, M. Cascante, P. R. D. A. Carulla, V. A. M. Dos Santos, and E. Saccenti, Metabolomics 14, 1 (2018).
  * Fransson and Marrelec (2008) P. Fransson and G. Marrelec, Neuroimage 42, 1178 (2008).
  * Bassett and Sporns (2017) D. S. Bassett and O. Sporns, Nature neuroscience 20, 353 (2017).
  * Tsonis  _et al._ (2006) A. A. Tsonis, K. L. Swanson, and P. J. Roebber, Bulletin of the American Meteorological Society 87, 585 (2006).
  * Donges  _et al._ (2009) J. F. Donges, Y. Zou, N. Marwan, and J. Kurths, The European Physical Journal Special Topics 174, 157 (2009).
  * Gozolchiani  _et al._ (2011) A. Gozolchiani, S. Havlin, and K. Yamasaki, Physical review letters 107, 148501 (2011).
  * Ozturk  _et al._ (2018) U. Ozturk, N. Marwan, O. Korup, H. Saito, A. Agarwal, M. Grossman, M. Zaiki, and J. Kurths, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 075301 (2018).
  * Ozturk  _et al._ (2019) U. Ozturk, N. Malik, K. Cheung, N. Marwan, and J. Kurths, Climate dynamics 53, 521 (2019).
  * Traxl  _et al._ (2016) D. Traxl, N. Boers, A. Rheinwalt, B. Goswami, and J. Kurths, Geophysical Research Letters 43, 9939 (2016).
  * Gupta  _et al._ (2021) S. Gupta, N. Boers, F. Pappenberger, and J. Kurths, Climate Dynamics 57, 3355 (2021).
  * Bockmayr  _et al._ (2013) M. Bockmayr, F. Klauschen, B. Györffy, C. Denkert, and J. Budczies, BMC systems biology 7, 1 (2013).
  * Newman (2016) M. E. J. Newman, in  _The New Palgrave Dictionary of Economics_ (Palgrave Macmillan UK, London, 2016) pp. 1–8.
  * Masuda  _et al._ (2018) N. Masuda, M. Sakaki, T. Ezaki, and T. Watanabe, Frontiers in Neuroinformatics 12, 7 (2018).
  * Bonett and Wright (2000) D. G. Bonett and T. A. Wright, Psychometrika 65, 23 (2000).
  * Hersbach  _et al._ (2020) H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers, _et al._ , Quarterly Journal of the Royal Meteorological Society 146, 1999 (2020).
  * Papoulis and Pillai (2002) A. Papoulis and S. U. Pillai,  _Probability, random variables, and stochastic processes_ , 4th ed. (McGraw-Hill Higher Education, 2002).
  * Mohapatra  _et al._ (2012) M. Mohapatra, B. Bandyopadhyay, and A. Tyagi, Natural Hazards 63, 1285 (2012).
  * Kendall (1938) M. G. Kendall, Biometrika 30, 81 (1938).
  * Amerise and Tarsitano (2015) I. L. Amerise and A. Tarsitano, Journal of Applied Statistics 42, 2584 (2015).
  * Adler (1957) L. M. Adler, Journal of the American Statistical Association 52, 33 (1957).
  * Kruskal (1958) W. H. Kruskal, Journal of the American Statistical Association 53, 814 (1958).
  * Gibbons and Chakraborti (2003) J. Gibbons and S. Chakraborti,  _Nonparametric Statistical Inference: Fourth Edition, Revised and Expanded. Statistics: A Series of Textbooks and Monographs_ (New York, NY: Marcel Dekker, Inc., 2003).
  * Goodman and Kruskal (1954) L. A. Goodman and W. H. Kruskal, Journal of the American Statistical Association 49, 732 (1954).
  * Xu  _et al._ (2013) W. Xu, Y. Hou, Y. Hung, and Y. Zou, Signal Processing 93, 261 (2013).
  * Note (1) We developed an optimized algorithm to compute Kendall’s rank correlation coefficient over a sliding time window, which benefits from the reuse of computation results obtained at the previous positions of the window.
  * Watts and Strogatz (1998) D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998).
  * Gilbert (1959) E. N. Gilbert, The Annals of Mathematical Statistics 30, 1141 (1959).
  * Erdős and Rényi (1959) P. Erdős and A. Rényi, Publicationes Mathematicae 6, 290 (1959).
  * Newman (2003) M. E. Newman, SIAM review 45, 167 (2003).
  * Barrat  _et al._ (2004) A. Barrat, M. Barthelemy, R. Pastor-Satorras, and A. Vespignani, Proceedings of the national academy of sciences 101, 3747 (2004).
  * Onnela  _et al._ (2005) J.-P. Onnela, J. Saramäki, J. Kertész, and K. Kaski, Physical Review E 71, 065103 (2005).
  * Zhang and Horvath (2005) B. Zhang and S. Horvath, Statistical applications in genetics and molecular biology 4 (2005).
  * Holme  _et al._ (2007) P. Holme, S. M. Park, B. J. Kim, and C. R. Edling, Physica A: Statistical Mechanics and its Applications 373, 821 (2007).
  * Opsahl and Panzarasa (2009) T. Opsahl and P. Panzarasa, Social networks 31, 155 (2009).
  * Kalna and Higham (2006) G. Kalna and D. J. Higham, in  _Symposium on network analysis in natural sciences and engineering_ (2006) p. 45.
  * Grindrod (2002) P. Grindrod, Physical Review E 66, 066702 (2002).
  * Ahnert  _et al._ (2007) S. E. Ahnert, D. Garlaschelli, T. M. Fink, and G. Caldarelli, Physical Review E 76, 016101 (2007).
  * Saramäki  _et al._ (2007) J. Saramäki, M. Kivelä, J.-P. Onnela, K. Kaski, and J. Kertesz, Physical Review E 75, 027105 (2007).
  * Wang  _et al._ (2017) Y. Wang, E. Ghumare, R. Vandenberghe, and P. Dupont, Neural computation 29, 313 (2017).
  * Whittaker and Whittaker (1990) J. Whittaker and J. Whittaker,  _Graphical models in applied multivariate statistics_ , Vol. 19 (Wiley Chichester, 1990).
  * Brier  _et al._ (2015) M. R. Brier, A. Mitra, J. E. McCarthy, B. M. Ances, and A. Z. Snyder, NeuroImage 121, 29 (2015).
  * Conover and Iman (1981) W. J. Conover and R. L. Iman, The American Statistician 35, 124 (1981).
  * Zhuang  _et al._ (2020) H. Zhuang, X. Wang, M. Bendersky, and M. Najork, in  _Proceedings of the 43rd international ACM SIGIR conference on research and development in information retrieval_ (2020) pp. 1649–1652.
  * Embrechts and Hofert (2013) P. Embrechts and M. Hofert, Mathematical Methods of Operations Research 77, 423 (2013).
  * Birkhoff (1931) G. D. Birkhoff, Proceedings of the National Academy of Sciences 17, 656 (1931).
  * von Neumann (1932) J. von Neumann, Proceedings of the National Academy of Sciences 18, 70 (1932).
  * Birkhoff (1939) G. Birkhoff, Duke Mathematical Journal 5, 19 (1939).
  * Boyles and Gardner (1983) R. Boyles and W. Gardner, IEEE Transactions on information theory 29, 105 (1983).
  * Anh and Lunney (1991) V. Anh and K. Lunney, Bulletin of the Australian Mathematical Society 44, 49 (1991).
  * Bhagavan (1985) C. Bhagavan, in  _Time Series in the Time Domain_, Handbook of Statistics, Vol. 5 (Elsevier, 1985) pp. 311–320.
  * Sugi  _et al._ (2002) M. Sugi, A. Noda, and N. Sato, Journal of the Meteorological Society of Japan 80, 249 (2002).
  * Cornfeld  _et al._ (2012) I. P. Cornfeld, S. V. Fomin, and Y. G. Sinai,  _Ergodic theory_ (Springer, 2012).
  * Coudène (2016) Y. Coudène,  _Ergodic theory and dynamical systems_ (Springer, 2016).
  * Wiggins and Ottino (2004) S. Wiggins and J. Ottino, Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 362, 937 (2004).
  * Tél  _et al._ (2020) T. Tél, T. Bódai, G. Drótos, T. Haszpra, M. Herein, B. Kaszás, and M. Vincze, Journal of Statistical Physics 179, 1496 (2020).
  * Wiener (1939) N. Wiener, Duke Math. J. 5, 1 (1939).
  * Pitt (1942) H. R. Pitt, Mathematical Proceedings of the Cambridge Philosophical Society 38, 325 (1942).
  * Becker (1981) M. E. Becker, The Annals of Probability 9, 504 (1981).
  * Rosenblatt (1956) M. Rosenblatt, The annals of mathematical statistics 27, 832 (1956).
  * Parzen (1962) E. Parzen, The annals of mathematical statistics 33, 1065 (1962).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Journal of Machine Learning Research 12, 2825 (2011).