New Exploratory Tools for Extremal Dependence: Chi Networks and Annual Extremal Networks

TL;DR

This paper introduces novel network-based tools, the Chi network and annual extremal network, for analyzing extremal dependence in environmental data, with applications to hurricane-related precipitation extremes.

Contribution

It develops new estimators and bias correction methods for extremal dependence networks, and applies them to real hurricane data to reveal spatial and temporal extremal dependence patterns.

Findings

Long-distance extremal dependence exists in precipitation extremes.

Dependence strength varies with regional meteorological conditions.

The methods effectively capture temporal and spatial extremal dependence.

Abstract

Understanding dependence structure among extreme values plays an important role in risk assessment in environmental studies. In this work we propose the $\chi$ network and the annual extremal network for exploring the extremal dependence structure of environmental processes. A $\chi$ network is constructed by connecting pairs whose estimated upper tail dependence coefficient, $\hat \chi$, exceeds a prescribed threshold. We develop an initial $\chi$ network estimator and we use a spatial block bootstrap to assess both the bias and variance of our estimator. We then develop a method to correct the bias of the initial estimator by incorporating the spatial structure in $\chi$. In addition to the $\chi$ network, which assesses spatial extremal dependence over an extended period of time, we further introduce an annual extremal network to explore the year-to-year temporal variation of extremal connections. We illustrate the $\chi$ and the annual extremal networks by analyzing the hurricane season maximum precipitation at the US Gulf Coast and surrounding area. Analysis suggests there exists long distance extremal dependence for precipitation extremes in the study region and the strength of the extremal dependence may depend on some regional scale meteorological conditions, for example, sea surface temperature.

Tables

Table 1. Table 1: Summary of TPR/PPV for the empirical χ ^ ^ 𝜒 \hat{\chi} and the bias-corrected χ ~ ~ 𝜒 \tilde{\chi} network (i.e., χ 𝜒 \chi > 0.3) estimators.

TPR percentile	5th	25th	50th	75th	95th
$\hat{\chi }$	0.76	0.80	0.83	0.86	0.91
$\tilde{\chi }$	0.88	0.97	1	1	1
PPV percentile	5th	25th	50th	75th	95th
$\hat{\chi }$	0.26	0.31	0.34	0.37	0.41
$\tilde{\chi }$	0.58	0.76	0.88	0.99	1

Equations

$$\chi =\underset{u\to 1^{-}}{\lim }\chi (u)=\underset{u\to 1^{-}}{\lim }\mathbb{P}(F_2(X_2)>u|F_1(X_1)>u),$$

$$\nu =\frac{1}{2}\mathbb{E}(|G_1(Z_1)-G_2(Z_2)|),$$

$$\hat{\nu }=\frac{1}{2}\frac{1}{m}\sum _{t=1}^m(|{\hat{G}}_1(z_{1,t})-{\hat{G}}_2(z_{2,t})|).$$

$$\theta =\frac{1+2\nu }{1-2\nu }.$$

$$\hat{\chi }=2-\frac{1+2\hat{\nu }}{1-2\hat{\nu }}.$$

$$\text{TPR}=\frac{\#\{(i,j):{\hat{\chi }}_{ij}>0.3\text{ and }{\chi }_{ij}>0.3\}}{\#\{(i,j):{\chi }_{ij}>0.3\}};$$

$$\text{PPV}=\frac{\#\{(i,j):{\hat{\chi }}_{ij}>0.3\text{ and }{\chi }_{ij}>0.3\}}{\#\{(i,j):{\hat{\chi }}_{ij}>0.3\}}.$$

$${\tilde{\chi }}_{ij}={\lambda }_{ij}{\hat{\chi }}_{ij}+(1-{\lambda }_{ij})\hat{\chi }(h_{ij}),i,j=1,\cdots ,d.$$

$${\hat{\chi }}_{ij}\sim \text{N}({\chi }_{ij},{\sigma }_{ij}^2);$$

$${\chi }_{ij}=\hat{\chi }(h_{ij})+{\omega }_{ij};$$

$${\omega }_{ij}\sim \text{N}(0,{\tau }^2(h_{ij})).$$

$$\text{log ratio}={\alpha }_0+{\alpha }_1\text{SST}+\epsilon .$$

$$\mathbb{E}[\log (\mu )|\text{SST}]={\beta }_0+{\beta }_1\text{SST}.$$

Full text

1 Introduction

The 2017 hurricane season was the costliest tropical cyclone season on record (Halverson, 2018). The season had 17 named storms, 10 hurricanes, and 6 major hurricanes. Great damage was done to US territories by three major hurricanes, Harvey, Irma, and Maria. These three storms were distinct both in location and timing: Harvey impacted Texas and Louisiana in mid-to-late August, Irma impacted Florida in early September, and Maria devastated Puerto Rico in late September. A standard treatment of the data from these three storms would likely treat them as independent events, due to their temporal and spatial distance. However, it is largely acknowledged that these storms are related, arising from conditions in 2017 conducive to tropical cyclone formation and intensification (Emanuel, 2011; Trenberth et al., 2018). While it is important to assess local impacts of each storm, understanding the connection of hurricane-induced extreme precipitation is also important for insurance companies or government agencies (e.g., Federal Emergency Management Agency (FEMA)) who need to prepare for the possibility of responding to multiple extreme events in a year. Thus we seek to use an approach that studies such phenomena in context of each other, specifically using a spatial network type approach.

The concept of networks has been used in climate science for about 15 years. The connection was first established by Tsonis and Roebber (Tsonis and Roebber, 2004; Tsonis et al., 2006), who introduced the concept of climate networks to study the dynamics of Earth’s climate system. A climate network consists of nodes (e.g., grid cells from a climate model output/gridded reanalysis data product) and edges. Two nodes are connected by an edge if the strength of a (pair-wise) dependence measure between the corresponding pairs of time-series is “strong” enough. For example, Tsonis and Roebber (2004) connected nodes with correlation greater than 0.5. This approach has been applied to many climate applications (e.g. Tsonis et al., 2006, 2008; Donges et al., 2009; Steinhaeuser et al., 2011) to gain insights into the dynamics of the climate system over many spatial scales, e.g., teleconnection paths.

Other important network types in climate science are event synchronization networks (Quiroga et al., 2002; Malik et al., 2010, 2012; Boers et al., 2014), where connections are defined based on whether extreme events at one point commonly coincide (within the span of a selected time interval) with extreme events at another point; phase synchronization networks (Yamasaki et al., 2009), where the signals at each point are viewed as oscillations, and connections are defined based on the apparent coupling between those oscillations; and causal networks (Runge, 2014; Ebert-Uphoff and Deng, 2014; Zerenner et al., 2014; Kretschmer et al., 2016), where connections are defined to represent potential cause-effect relationships. Causal networks can be established by many different methods. For example, they can be based on probabilistic graphical models (Pearl, 2000; Spirtes et al., 2000) (for their use in climate science, see for example (Ebert-Uphoff and Deng, 2014)); on combining the concept of Granger causality with vector autoregression (VAR) (Lütkepohl, 2007) or LASSO (Tibshirani, 1996) methods, as known as Graphical Lasso or GLasso, (Yuan and Lin, 2007; Friedman et al., 2008) (for their use in climate science see for example (Bahadori and Liu, 2011)); or on Gaussian models (Zerenner et al., 2014). Note that the majority of these network types are closely related to the concept of correlation. Namely, correlation networks use correlation itself (Tsonis and Roebber, 2004); many probabilistic graphical models are derived using partial correlation (Ebert-Uphoff and Deng, 2014);VAR and GLasso approaches approximate regression coefficients, which can be expressed in terms of partial correlation; and Gaussian models approximate the precision matrix (i.e., inverse covariance matrix), whose elements can be interpreted as partial correlation coefficients. Only two network types mentioned above are not based on correlation, namely phase synchronization networks, which with their assumption of oscillators do not appear to be a good match to model extremes, and event synchronization networks, which are in fact most closely related to the concepts discussed here, but are not based on well understood statistical measures.

Correlation measures the linear relationship between a bivariate random vector around the center of its distribution and hence is poorly suited for describing dependence in the tail, i.e., the dependence between extreme events. In contrast, several dependence measures designed specifically for extremes have been proposed in the literature and those are likely to be more appropriate than correlation to construct networks that focus on analyzing relationships between extremes. Three examples are the upper tail dependence, aka tail dependence coefficient ($\chi$) (Coles et al., 1999), the extremal coefficient ($\theta$) (Smith, 1990; Schlather and Tawn, 2003), and the F-madogram ($\nu$) (Cooley et al., 2006; Naveau et al., 2009). See Chapter 2 of Yan and Dey (2015) for a review. The upper tail dependence is defined as

[TABLE]

where $X_{1}$ and $X_{2}$ are two random variables and $F_{1},F_{2}$ their corresponding cumulative distribution function (CDF). Thus the $\chi$ measure is the (limiting) conditional probability of one random variable being extreme, given that the other random variable is extreme. Note that the $\chi$ measure is symmetric in terms of which variable is chosen to be conditioned upon, due to the standardization of the marginal CDFs.

In Section 2, we introduce a new network type, the $\chi$ network, where we use the aforementioned upper tail dependence ($\chi$) to explore the spatial dependence structure in extremes. We briefly proposed this idea at a workshop (Ebert-Uphoff et al. (2018)), but it is fully developed for the first time here. Specifically, we construct a $\chi$ network by connecting the pairs of nodes whose $\hat{\chi}\text{'s}$ exceed some specified value (see Sec. 2.1 for details). However, we find that a straight forward implementation of this procedure does not perform well as it introduces many spurious network connections (see Sec. 2.2). This problem, which can be described as a network bias, is not unique to $\chi$ networks. In fact, we found that the same bias can occur in traditional correlation networks as well, although, to the best of our knowledge, this issue has never been discussed for correlation networks. In Sec. 2.3 we describe a method to alleviate this issue by incorporating the spatial structure on $\chi$ observed from the empirical estimates to reconstruct the $\chi$ network. Although we did not test this hypothesis yet, we conjecture that this bias correction method could be used to correct the bias in correlation networks as well.

The $\chi$ network assesses extremal dependence integrated over time, which is useful for assessing risk. However, as such it cannot be used to explore the “dynamics” (i.e., year-to-year) structure of extreme events in order to gain insights about the potential mechanisms that contribute to extreme behavior, such as the 2017 hurricane season. In Section 2.4 we develop another network type, the annual extremal network, that aims to explore the potential drivers of extreme behavior.

In Section 3, we create a $\chi$ network to hurricane-season maximum daily precipitation data collected in the US Gulf Coast to explore the spatial pattern of the extremal network connections. We also apply the annual extremal network to this precipitation data and examine whether sea surface temperature affects characteristics of the annual networks.

It should be noted that the purpose of the new network types proposed here is not to replace existing network types, such as correlation networks or event synchronization networks, but to supplement them. Namely, for any given application one might construct both a correlation network, which is more sensitive to identifying dependencies near the mean of the distribution, and a $\chi$ network, which is more sensitive to identifying dependencies between extremes. The two different networks are thus expected to provide two complementary viewpoints of the spatial dependencies. However, comparing these types of networks is not trivial, so is subject of future research, and here we simply focus on presenting the new type of network.

2 $\chi$ and Annual Extremal Networks

In this section we first describe how to construct a $\chi$ network, but also in tandem describe an illustrative simulation study. Our reason for interweaving these two discussions is that the simulation study exposes a bias issue with the $\chi$ network that arises due to the act of thresholding at a level of $\chi$ and creating edges for pairs of nodes whose estimated $\hat{\chi}$ exceeds this threshold. We will begin by introducing the needed background in order to describe the procedure for estimating $\chi$ for each pair of nodes/locations in Section 2.1. We then present a simulation study to illustrate the issue of network bias (Sec. 2.2), which we then correct by conducting a spatial regularization on $\chi$ (Sec. 2.3). We conclude the section by describing the annual extremal networks (Sec. 2.4).

2.1 An Empirical Estimator of $\chi$

We first discuss how we estimate $\chi$ from observations of a random sample of a bivariate random vector. To develop an estimator for $\chi$ one would need to first estimate the marginal distributions. Second, since $\chi$ is the limiting (conditional) probability, one would need to choose a threshold $u\in(0,1)$ to estimate $\chi$. There is a bias-variance trade-off to be made in the choice of $u$. Using a $u$ that is not sufficiently high results in increased bias, while using a $u$ that is too high results in increased variance due to paucity of data. In practice, often $\hat{\chi}(u)$ is estimated for many values of $u$ and plotted (see, for example, the chiplot function of the evd package (Stephenson, 2002) for R (R Core Team, 2018)).

Our particular application, with an additional assumption of max-stability, allows us to avoid the selection of $u$. Our motivation here is to better understand relationships between hurricane events in a given year, which leads us to investigate annual extremes. To make this precise, let $\bm{z}_{t}=(z_{1,t},z_{2,t})^{\text{T}}=\left(\max_{i=(t-1)*b+1}^{b*t}(y_{1,i}),\max_{i=(t-1)*b+1}^{b*t}(y_{2,i})\right)^{\text{T}},t=1,\cdots,m=n/b$, where $\bm{y_{i}}=(y_{1,i},y_{2,i})^{\text{T}},i=1,\cdots,n$ be the original observations with sample size $n$, $b$ be the block size, and $m=n/b$ be the sample size for the block maximum. For example, $\{\bm{y}_{i}\}_{i=1}^{n}$ could denote the daily rainfall amounts at two locations, and $i$ the time index (e.g., starting from January 1st for some year). If we use $b=365$ (or $366$ for leap years), then $\bm{z}_{t}$ would be the annual maximum precipitation amounts at these two locations in year $t$, where their dependence is of interest for risk assessment.

Below, we use estimators of extremal dependence that have been proposed for the max-stable setting to obtain our estimator for $\chi$. Max-stability is an important notion that underlies much of extreme value analysis and here we give a brief introduction, The reader is referred to (Coles et al., 2001; Beirlant et al., 2004; de Haan and Ferreira, 2006; Segers, 2012) for more details. A multivariate distribution $\bm{Z}\in\mathbb{R}^{d}$ is max–stable if the distribution is close under (component-wise) maximization. That is, each marginal distribution, by Extremal Types Theorem (Fisher and Tippett, 1928; Gnedenko, 1943), is belongs to the family of generalized extreme value (GEV) (Jenkinson, 1955; Coles et al., 2001). Furthermore, the dependence structure, which can be described by its copula $C_{*}$ (Joe, 1997; Nelsen, 2007) is satisfied the property that $C_{*}(u_{1}^{1/N},\cdots,u_{d}^{1/N})^{N}=C_{*}(u_{1},\cdots,u_{d})$ for all $N>0$. The copula $C_{*}$ is called max-stable copula which form the basis for much of multivariate and spatial extreme value analysis (Davison et al., 2012; Ribatet and Sedki, 2013). The implication of max-stability is that any pairwise $\chi$ for a max-stable random vector is a constant with respect to $u$ and hence there is no need to choose a threshold. Moreover, there exists deterministic relationships between $\chi$ and the other aforementioned summary statistics for extremal dependence, namely, $\theta$ and $\nu$ under the max-stable condition. We will make use of these relationships to develop our estimator for $\chi$ in the next paragraph.

We first estimate the marginal distributions of $Z_{1}$ and $Z_{2}$, denoted as $G_{1},G_{2}$, non-parametrically by their empirical distribution functions (EDFs) $\hat{G}_{1},\hat{G}_{2}$. We then estimate the F-madogram, the first order absolute difference in CDF

[TABLE]

using the following empirical estimator

[TABLE]

(Cooley et al., 2006) has shown that there is a one to one correspondence between $\nu$ and $\theta$

[TABLE]

Furthermore, by the max-stability assumption, we have $\chi=2-\theta$. Therefore we use the plug-in estimate $\hat{\theta}=(1+2\hat{\nu})/(1-2\hat{\nu})$ to obtain the resulting estimator of $\chi$:

[TABLE]

To create a $\chi$ network, we first apply this estimator to all the pairwise components $\{(i,j)\}$ of a $d$ dimensional random vector $\bm{Z}$. We then connect those pairs with $\{\hat{\chi}_{ij}\}_{i=1,\cdots,d}^{j=1,\cdots,d},i\neq j$ greater than a chosen threshold $\chi_{min}$ to construct the corresponding $\chi$ network.

2.2 Simulation Study

The purpose of this simulation study is to examine the performance of the empirical estimator introduced in the previous section for estimating the $\chi$ network, where the network consists of pairs with $\chi$ greater than some pre-specified threshold. To do that we will need to generate data from a random vector for which we know the true network for a given threshold. The estimated network can then be compared with the true network for evaluation. This simulated network thus serves as a benchmark for the proposed new network estimation procedure.

To make this simulation relevant in the context of environmental applications, we first simulate 100 spatial locations uniformly in a $[0,1]\times[0,1]$ domain and these are the nodes of the network. We then generate realizations from a stationary and isotropic Brown-Resnick process (Brown and Resnick, 1977; Kabluchko et al., 2009) at these locations. A Brown–Resnick process is a flexible class of stationary max-stable processes where every finite collection of those random variables has a multivariate max-stable distribution. The dependence structure of a Brown–Resnick process is induced by a Gaussian process via the spectral representation of max-stable process (de Haan, 1984). The interested reader is referred to (Kabluchko et al., 2009; Huser and Davison, 2013; Thibaud et al., 2016) for more details. The implied $\chi$ value of a Brown–Resnick process for each pair $(i,j)$ only depends on their spatial distance, denoted by $h_{ij}$, and the (Gaussian) process parameters $\rho$, the range, and $\kappa$, the smoothness: $\chi(h)=2-2\Phi(\sqrt{\frac{h^{\kappa}}{2\rho}})$, where $\Phi$ is the CDF of the standard normal distribution. Setting $\rho=0.05$, $\kappa=1$, and $\chi_{min}=0.3$, two nodes will be connected if their distance is less than 0.107, a relatively small distance compare with the spatial domain in our simulation study.

We simulate 100 realizations (the Monte Carlo sample size), each with sample size $m=50$ and dimension $d=100$ (i.e., a Brown-Resnick process observed at the aforementioned 100 spatial locations). For each realization in the simulation study, we estimate the $\chi$ network with $\chi_{min}=0.3$ using the empirical estimator described in Sec. 2.1. Comparing the estimated networks and the true network, it is apparent that the number of edges is overestimated (see Fig. 1 (a), (b) and (c)). We further evaluate the estimation performance in terms of the true positive rate (TPR) and positive predictive value (PPV), defined as

[TABLE]

[TABLE]

While TPR is reasonably high ($\sim$ 80 %), the estimation procedure identifies many more false positives than true positives as the PPV is rather low ($\sim 35\%$, see Fig. 1 (d)).

Fig. 2 provides insight into how this biased estimate of the number of connections arises. Plotted is the true $\chi$ as a function of spatial distance along with one simulation’s $\chi$ estimates. Although the estimates themselves appear unbiased as they “center” around the true values, it is the act of thresholding which introduces the bias. Clearly there are more false positives than false negatives due to the shape of the true $\chi$ function and the increasing variance of the estimates as $\chi$ decreases. Note that this problem is not unique to the $\chi$ networks proposed here. In fact, we confirmed in another simulation that the same effect can also occur for the well-established correlation networks, although, to the best of our knowledge, this issue has never been mentioned in the literature on climate networks.

2.3 Network Bias Correction

Here we develop a method to correct the bias of the network estimation as demonstrated in the previous subsection. We exploit the spatial structure revealed from $\{\hat{\chi}_{ij}\}$ (see Fig. 2) to spatially regularize the network estimate. In what follows we will describe a new estimator for $\{\chi_{ij}\}$, denoted by $\{\tilde{\chi}_{ij}\}$, which is a weighted average of empirical estimator $(\{\hat{\chi}_{ij}\})$ and an estimator that explicitly incorporates spatial information $(\hat{\chi}(h_{ij}))$. The weight for a given pair $(i,j)$, denoted as $\lambda_{ij}$, will depend on the estimation precision of $\hat{\chi}_{ij}$ and $\hat{\chi}(h_{ij})$. The resulting estimator, the weighted average of $\hat{\chi}_{ij}$ and $\hat{\chi}(h_{ij})$, has an empirical Bayes interpretation (see Loader and Switzer (1992) in the context of spatial covariance function estimation), can then be used to construct the $\chi$ network. Specifically, the bias-corrected estimator has the following form:

[TABLE]

In the following we first describe how to estimate $\{\chi(h_{ij})\}$ and then discuss how to come up with the weights $\{\lambda_{ij}\}$.

2.3.1 The estimator for $\chi(h)$

We smooth the $\{\hat{\chi}_{ij}\}$, the empirical estimates for all the pairs, by spatial distance to obtain $\hat{\chi}(h)$. Specifically, we divide $\{\hat{\chi}_{ij}\}$ into $K$ distance bins and we compute the mean value of $\{\hat{\chi}_{ij}\}$ within each bin, denoted by $\hat{\chi}^{k}$. We then fit a cubic smoothing spline regression (Wahba, 1990; Gu, 2013) to $\{h^{k},\hat{\chi}^{k}\}_{k=1}^{K}$, where $h^{k}$ denotes the average distance within the $k_{\text{th}}$ bin, using the sreg function in the fields R package (Nychka et al., 2015) to get the $\hat{\chi}(h)$ as a function of spatial distance only. Note that this approach effectively introduced a spatially stationary and isotropic structure on $\chi$, which is less tenable especially when $h$ is “large”. Therefore, we assume for any $h>0$, $\hat{\chi}(h)$ is normally distributed with variance $\tau^{2}(h)$ increases with spatial distance in order to weaken this somewhat strong and rigid assumption in order to learn potential nonstationary pattern on $\chi$.

2.3.2 The weights $\{\lambda_{ij}\}$

The weight $\lambda_{ij}$ for a given pair $(i,j)$ is determined by the estimation variability of $\hat{\chi}_{ij}$ ($\sigma^{2}_{ij}$), and the amount of departure from $\hat{\chi}(h)$, quantified by a variance function: $\tau^{2}(h)$. We perform a “spatial” block bootstrap (see e.g., Wang et al. (2016), Appendix C, or Huang et al. (2016), Appendix D) by re-sampling the vectors $\{\bm{z}_{t}\}_{t=1}^{m}$ to obtain the bootstrap variances $\{\sigma^{2*}_{ij}\}$ as the estimates for the corresponding $\{\sigma^{2}_{ij}\}$. We then use $\tau^{2}_{ij}/(\tau^{2}_{ij}+\sigma^{2*}_{ij})$ as the weight $\lambda_{ij}$.

The estimator (8) can be formulated as the posterior mean from the hierarchical model

[TABLE]

$\tilde{\chi}_{ij}$ can be interpreted as an (empirical) Bayes (Robbins, 1955; Efron and Morris, 1972; Casella, 1985) estimator, where the empirical aspect arises from the fact that we estimate $\chi(h)$ from $\{\hat{\chi}_{ij}\}$.

2.3.3 An illustration

Following the bias correction procedure described in previous subsections, we first bin the $\{\hat{\chi}_{ij}\}$ into $K=100$ bins to obtain $\{\hat{\chi}^{k}\}_{k=1}^{100}$. We then fit a cubic spline regression to $\{h^{k},\hat{\chi}^{k}\}_{k=1}^{100}$ to get the estimate of $\chi(h)$ (see Fig. 3, (a)), we then assume $\tau^{2}(h)=0.095/\left(1+\exp\left(-6*\left(h-0.72\right)\right)\right)$, which is an increasing (logistic) function with $h$. As described in previous subsection, the bootstrap standard deviations are used as the estimates of $\{\sigma_{ij}\}$ and we compare those estimates with the “oracle” standard errors (i.e., standard deviations (SDs) among those 100 Monte Carlo simulations from the true model). As would be expected the bootstrap SDs are (slightly) higher than their corresponding oracle SDs but they show similar behavior, at least qualitatively. We can then plug in the $\sigma^{2*}_{ij}$, $\hat{\chi}(h_{ij})$, and $\tau^{2}(h_{ij})$ to obtain the resulting estimate $\tilde{\chi}_{ij}$ for each pair to construct the bias-corrected network. This estimator performs well as $\{\tilde{\chi}_{ij}\}$ have been “shrunk” toward $\hat{\chi}(h_{ij})$, which is much closer to the truth than $\{\hat{\chi}_{ij}\}$ (Fig. 3, (c)).

The comparison of TPR/PPV for the empirical and the bias-corrected $\chi$ network estimators shows that our method not only identifies most of the true connections, i.e., TPR $\approx$ 1, but also avoids making many false discoveries, i.e., PPV for the bias-corrected estimator is substantially higher than that of the empirical estimator (see Tab. 1).

This bias correction method, with some modifications, is expected to apply also to correlation networks to compensate for their network bias, but that topic is beyond the scope of this work.

2.4 Annual Extremal Network

The estimation procedure of the $\chi$ network and its bias-corrected version assumes that the dependence structure of the block maximum, $\{\bm{Z}_{t}\}$, does not change with “time” (i.e., $\bm{z}_{t}$, $t=1,\cdots,m$ is a random sample from an underlying time invariant random vector $\bm{Z}$), so that one can use these “replicates” to infer the spatial extremal dependence. Moreover, we assume that $\bm{Z}$ is max-stable to avoid the need to choose a high threshold $u$ to estimate $\chi$. Here we introduce another type of network, the annual extremal network, to examine the year-to-year variation of the “observed” extremal networks and we seek to explain this variability. The main idea of the annual extremal network is to connect a pair of nodes in a given “year” $t$ if extreme events occur at both locations. For example, an extreme event can be defined as the value exceeding a high percentile, say the 95th percentile. In this case, two locations are connected if both locations exceed their 20-year return level in a given “year”. Specifically, we treat the EDFs of $\{\bm{z}_{t}\}$, denoted by $\{\bm{u}_{t}\}$, where $\bm{u}_{t}=\{u_{1,t},u_{2,t},\cdots,u_{d,t}\}^{T}$, as the “data” and for each $t=1,\cdots,m$, we connect pairs $(i,j)$ if both $u_{i,t}$ and $u_{j,t}$ are greater than $0.95$.

The definition of the annual extremal network defined here has similaraties to the definition of event synchronization networks (Malik et al., 2010, 2012; Boers et al., 2014), but the latter operate on a completely different time scale. Event synchronization networks in climate science typically consider daily data, and define extreme events at each location to be those exceeding a certain threshold, e.g. the 95th percentile of the overall distribution. Two locations are connected if an extreme event at one location is frequently followed by an extreme event at another location within a matter of days, which can be used to track the propagation of individual phenomena from one location to another, e.g., the propagation of storm systems. In contrast, the annual extremal network, due to its annual scale, investigates whether the overall conditions in a given year lead to several events in the same year, although those events may not have any direct connection to each other. Thus, due to the different time scale (days vs. years) these two network types focus on very different phenomena.

Summary statistics (e.g., the number of the connections, the average distance of connections at each location) of this annual network are time series (on an annual scale) and can be related to some physically meaningful covariates to explore its year-to-year variation. In Section 3.3 we use the number of “long distance” connected pairs to relate this statistic to sea surface temperature. Note that the annual extremal network described here is constructed directly from the “data”, namely the EDFs for all the locations, whereas the $\chi$ network is constructed from the $\chi$ statistic that was estimated from the EDFs. Therefore, by ignoring the estimation uncertainty of EDFs themselves, we do not need to develop a bias correction procedure as we did for the $\chi$ network.

3 Gulf Coast Extreme Precipitation Application

3.1 Data Set

We use a subset of the Global Historical Climatology Network (GHCN) data set (Menne et al., 2012), where the spatial domain includes Texas, Louisiana, Mississippi, Alabama, Florida, and Georgia and the time period is from Jan. 1949 to Oct. 2017. There are 339 GHCN stations being included in this study where the selecting criterion is that those stations should have at least 90 percentage of non-missing daily precipitation values. We further restrict our attention to the annual maximum of the daily precipitation values from June to October to investigate the extremal dependence structure for the hurricane season.

3.2 $\chi$ Network for Gulf Coast Extreme Precipitation

We first compute the $\chi$ network based on the empirical estimator as described in Sec. 2.1 and we assess the network estimation uncertainty using a spatial block bootstrap. The results obtained from the block bootstrap (Fig. 4, Left panel) clearly demonstrate the issue of overestimation. We then conduct a bias correction using the method described in Sec. 2.3. The right panel in Fig. 4 displays the $\{\hat{\chi}_{ij}\}$, $\chi(h)$, and $\{\tilde{\chi}_{ij}\}$ as a function of the spatial distance.

Fig. 5 shows that many long distance links in the $\chi$ network based on the empirical estimator disappear and may simply be due to estimation variability for using the empirical estimator. However, the bias-corrected $\chi$ network does identify some long distance edges, for example between the Houston and Tampa regions. It is important to note that although $\chi$ takes on values between 0 and 1, its values should not be interpreted in a similar manner to correlation values. A $\chi$ of 0.3 between two distant locations would imply strong tail dependence; that is, when one location experiences a very extreme event, the other would have a 30 percent chance of also experiencing an extreme event. Individual links are hard to interpret, and we do not believe there is some phenomenon that specifically links Houston to Tampa. Nevertheless, these long-distance network links at a minimum suggest that calculating annual risk between two distant regions like Houston and Tampa should not be done with an independence assumption as this would underestimate joint risk.

3.3 Annual Extremal Precipitation Network for Gulf Coast

In this subsection we construct annual extremal networks in the study region to examine the year-to-year structure. Specifically, for each hurricane season we identify the connections where their EDFs of the hurricane season maximum exceed 0.95 (i.e., both places experience at least a 20-year event in that year, see Fig. 6 for examples). To further sutdy the long range extremal dependence, we compute the number of pairs where their spatial distance are greater than 1,000km. We then study how the inter-annual variability of the numbers of these “long distance” extremal pairs might be explained by some meteorological covariate, for example, annual mean sea surface temperature (SST) averaged in a relevant spatial region. Here we use the Hadley Centre Global monthly average Sea Surface Temperature version 1.1 (HadISST 1.1, Rayner et al. (2003)) and we compute the spatial average of the annual mean SST over a “rectangular” region with longitude from $95^{\circ}$ west to $83^{\circ}$ west and latitude from $23^{\circ}$ north to $29^{\circ}$ north where the annual window is defined as from July of the previous year to June (i.e., the 2017 mean SST is computed by averaging the monthly SST from July 2016 to June 2017).

We explore the potential cause of the annual variation in the number of extremal connections as shown in Fig. 6. In Fig. 7 (left panel), we plot two time series together, the red one for the proportion of the long distance (i.e., > 1,000km) extremal connections on a log scale (log ratio hereafter), where this log ratio is computed as the logarithm of the number of the connections divided by the total number of the long distance pairs in the study region in a given year. The blue line is the SST time series. There exists “co-movement” between these two time series motivated us to perform regression analyses where we use the log ratio as the response and SST as the predictor (see Fig. 7, Right panel). Specifically, we perform two regression analyses to study this relationship. First, a linear regression (lm):

[TABLE]

Second, a generalized linear model (McCullagh and Nelder, 1989) assuming the number of the long distance extremal connections follows a Poisson distribution with its log intensity, denoted by $\log(\mu)$, depends linearly on SST (glm):

[TABLE]

Both the lm (although the homoscedasticity assumption is likely not appropriate here) and glm model fits suggest that the SST can explain some proportion of the year-to-year variation observed in the number of long distance extremal pairs, and the estimated slopes ($\hat{\alpha}_{1}$ and $\hat{\beta}_{1}$) are positive with p-values $0.008$ and $<2^{-16}$. The results here suggest that the annual mean Gulf Coast SST may partly contribute to the abnormally long distance extremal dependence observed in the 2017 hurricane season.

3.4 Scientific Implications

More intense hurricanes (i.e., named storms) under climate change lead to great US economic losses (e.g., Emanuel, 2011; Klotzbach et al., 2018). Studies already found that number of extreme precipitation events associated with hurricane has increased in recent decades (e.g., Kunkel et al., 2010). But it is still largely unknown how the spatial dependence of extreme precipitation may respond. Sparse summer convective precipitation events have typical radius around 10-100km, which lead to short-distance spatial dependence. However, major Atlantic hurricanes can grow to 800km measured with the outer radius of vanishing wind (Chavas et al., 2016). Such large low-pressure system deliver huge amount of water when it moves from the warm ocean to land. In this case, faraway places like some regions in Texas and Florida could all encounter extreme precipitation during a hurricane season. Here our biase-corrected $\chi$ network (see Fig. 5, Bottom panel) captures a range of precipitation spatial dependence, not only the common short-distance connectivity among nearby stations but also the long-distance connectivity across states associated with hurricane events. With this overall understanding for the spatial extremal dependence, we further examine the temporal variation of the annual extremal precipitation networks and we investigate the potential relationship with the US Gulf Coast SST. The analyses in Sec. 3.3 reveal there exists a robust relationship between the number of long distance extremal precipitation connection in southern states of US and pre-hurricane season (previous year July to current year June) SST in Gulf of Mexico. It agrees with the understanding that increased ocean heat content leads to more tropical cyclone (Trenberth et al., 2018), which in turn induces more long-distance preciptation extremes. In summary, by combining the extremal dependence estimates with the network analysis, we are able to quantitatively measure how the hurricane-induced extreme precipitation changes under climate change. The network perspective enables us to reveal the structural change among a suite of stations, which provide additional information to standard analysis based on individual locations.

4 Discussion

This work presents two new tools, the $\chi$ and annual extremal networks, for exploring dependence structure of extreme values for space-time data sets. In contrast to the widely used correlation network in the climate literature, the $\chi$ network uses $\chi$ statistics to summarize extremal dependence, and is hence suitable for studying extremal dependence. We demonstrate through a simulation study that a $\chi$ network constructed by using an empirical estimator could overestimate the degree of the $\chi$ network. We found that this issue can also arise if one uses the sample correlation to construct a correlation network. We propose a bias-correction method to incorporate the spatial dependence structure to alleviate the network estimation bias for the $\chi$ network. We expect that this bias correction procedure could also be used to overcome the bias in traditional correlation networks. We then apply the $\chi$ network to the Gulf Coast hurricane season rainfall extremes to study the spatial dependence structure of the annual maximum precipitation observed over the Gulf Coast.

We also apply the annual extremal network to further explore the potential non-stationary long-range extremal dependence. Our results suggest that the sea surface temperature may contribute to the abnormally long distance extremal dependence observed in the 2017 hurricane season.

Due to the exploratory nature of the $\chi$ network, there are some limitations in terms of providing the whole picture of the extremal dependence structure of an environmental process of interest. First, the $\chi$ statistic is a measure of pairwise tail dependence and hence not able to describe the higher order extremal dependence. Second, $\chi$ is an appropriate extremal dependence metric in the case of asymptotic dependence, it is less clear about the usefulness of this measure if the environmental process of interest exhibit weakening spatial dependence as events become more extreme (Huser et al., 2017; Wadsworth et al., 2017; Huser and Wadsworth, 2018; Bopp et al., 2018).

In this work, we make several simplifying assumptions for developing our bias-corrected method in order to obtain a close form expression for our estimator. One could try to relax these assumptions at the price that the computations needed for a sampling based approach, which can be substantial. Given that the number of pairs will be large for most environmental applications, we prefer a simple method such as the one we proposed here, that allows us to quickly explore the extremal dependence structure, which could provide some guidance on how to proceed the following confirmatory analysis.