A Hierarchical Max-infinitely Divisible Process for Extreme Areal Precipitation Over Watersheds

Understanding the spatial extent of extreme precipitation is necessary for determining flood risk and adequately designing infrastructure (e.g., stormwater pipes) to withstand such hazards. While environmental phenomena typically exhibit weakening spatial dependence at increasingly extreme levels, limiting max-stable process models for block maxima have a rigid dependence structure that does not capture this type of behavior. We propose a flexible Bayesian model from a broader family of max-infinitely divisible processes that allows for weakening spatial dependence at increasingly extreme levels, and due to a hierarchical representation of the likelihood in terms of random effects, our inference approach scales to large datasets. The proposed model is constructed using flexible random basis functions that are estimated from the data, allowing for straightforward inspection of the predominant spatial patterns of extremes. In addition, the described process possesses max-stability as a special case, making inference on the tail dependence class possible. We apply our model to extreme precipitation in eastern North America, and show that the proposed model adequately captures the extremal behavior of the data.



page 8

page 15

page 16


Penultimate modeling of spatial extremes: statistical inference for max-infinitely divisible processes

Extreme-value theory for stochastic processes has motivated the statisti...

Exploration and inference in spatial extremes using empirical basis functions

Statistical methods for inference on spatial extremes of large datasets ...

A multivariate spatial skew-t process for joint modeling of extreme precipitation indexes

To study trends in extreme precipitation across US over the years 1951-2...

Bayesian Modeling of Air Pollution Extremes Using Nested Multivariate Max-Stable Processes

Capturing the potentially strong dependence among the peak concentration...

Modeling short-ranged dependence in block extrema with application to polar temperature data

The block maxima approach is an important method in univariate extreme v...

Projecting Flood-Inducing Precipitation with a Bayesian Analogue Model

The hazard of pluvial flooding is largely influenced by the spatial and ...

A Model-Free Selection Criterion For The Mixing Coefficient Of Spatial Max-Mixture Models

One of the main concerns in extreme value theory is to quantify the depe...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The risk of precipitation-induced flooding (pluvial flooding) is strongly determined by the spatial extent of severe storms, and therefore, there is a need to adequately describe the spatial dependence properties of extreme precipitation. With this goal in mind, we propose a scalable model for spatial extremes that relaxes the rigid dependence structure of asymptotic max-stable models, characterizes the main modes of spatial variability using interpretable spatial factors, and allows for easy prediction at unobserved locations. The areal aspect of extreme precipitation plays a role in flood risk assessment. Precipitation falling over a single drainage basin flows into a common outlet, the aggregate effects of which can be devastating in large volumes. In 2006, heavy precipitation over the Susquehanna River basin in New York and Pennsylvania caused record high discharges along the Susquehanna River and flooding in the region, ultimately leading to federal-level disaster declarations and disaster-recovery assistance from the US Federal Emergency Management Agency (FEMA) in excess of $227 million (Suro et al., 2009).

The last decade has seen a considerable amount of research on the spatial dependence modeling of extremes, in part because of the hazard that extreme weather events pose to human life and property. For recent reviews, see Davison et al. (2012, 2013, 2019) and Davison and Huser (2015). The classical geostatistical Gaussian process models that are ideal for modeling the bulk of a distribution have weak tail-dependence and do not enforce the specific type of positive dependence structure inherent to extremes. Two classes of models, max-stable processes (de Haan and Ferreira, 2006) and generalized Pareto processes (Ferreira and de Haan, 2014; Thibaud and Opitz, 2015)

, have proven to be useful tools for the modeling of spatial extremes. Max-stable process models are infinite-dimensional generalizations of the limiting models for componentwise maxima. They are asymptotically justified models for pointwise maxima over an infinite collection of independent processes after suitable renormalization, a property which has made them prime candidates for the modeling of spatial extremes. In practice, maxima are taken over large, but finite blocks (e.g., months, years). An approximation error is incurred when applying limiting models to pointwise maxima over finite blocks, and the degree of this error will depend on the rate of convergence of the modeled process as the block size grows. Furthermore, the approximation error is more pronounced when the observed process exhibits weakening spatial dependence at increasingly high quantiles, as the spatial dependence of limiting max-stable processes is the same across all levels of the distribution, and hence would overestimate the level of dependence in the data. For more discussion, see, e.g.,

Wadsworth and Tawn (2012). Empirical evidence has shown that environmental processes often exhibit weakening spatial dependence at more extreme levels (Huser et al., 2017, 2018; Huser and Wadsworth, 2019). Misspecifying the dependence structure can have more significant consequences when estimating aggregate quantities such as the total precipitation volume discharged by the collection of annual maxima over a region.

In this paper, we aim to extend a class of max-stable models, in order to flexibly capture spatial dependence characteristics for sub-asymptotic block maxima data, while still retaining the positive dependence structure inherent to distributions for maxima. The general class of models that we consider, which nests the class of max-stable models, are known as max-infinitely divisible (max-id) processes (Resnick, 1987, Chapter 5)

. Suppose a random vector

has joint distribution

, then the distribution of maxima of independent and identically distributed (i.i.d.) replicates , taken componentwise, has distribution function . The max-id property applies to the converse statement. Suppose that is a random vector of componentwise maxima, composed from a collection of i.i.d. vectors. Then if has distribution function , there exists some root distribution such that , or equivalently such that . By continuous extension of the relation for , we say that a distribution is max-id if and only if is a valid distribution for all real . This is always the case for univariate distributions, but may not necessarily be so for multivariate distributions. An appealing property of max-id distributions is that they allow a change of temporal support. For example, from a fitted max-id model to annual maxima, conclusions might directly be drawn for monthly maxima (as is a valid distribution), or even daily maxima ( is also valid), provided temporal dependence and non-stationarity have properly been taken into account. Informally, max-id distributions are those which arise from taking componentwise maxima of i.i.d. random vectors and are therefore an appropriate class to constrain ourselves to if the goal is to model componentwise maxima. By slight abuse of language, we say that a spatial process is max-id if all its finite-dimensional distributions are max-id. Necessary and sufficient conditions for max-infinitely divisibility of a distribution function in were first given by Balkema and Resnick (1977). More recently, mixing conditions for stationary max-id processes were explored by Kabluchko and Schlather (2010), and minimality of their spectral representations were described in Kabluchko and Stoev (2016).

Unlike limiting max-stable process models, which have a rigid spatial dependence structure, sub-families of the broader class of max-id processes do not impose such constraints and can accommodate different spatial dependence characteristics across various levels of a distribution see, e.g. Padoan, 2013. It is this feature that can cause max-stable processes to fit poorly, as asymptotically independent processes may exhibit spatial dependence at finite levels. Extrapolation of max-stable fits to higher quantiles in this scenario can cause overestimation of the risk of concurrent extremes (Davison et al., 2013). Furthermore, the challenge of performing conditional simulation from max-stable models given observed values at many locations is a limiting factor for their use in practice (Dombry et al., 2013). The Bayesian model that we develop in the remainder of the paper permits a conditional, hierarchical representation in terms of random effects that facilitates fast conditional simulation, which is useful for prediction at unobserved locations, and for handling missing values.

2 Hierarchical Construction of Spatial Max-Id Models

2.1 Max-Stable Reich and Shaby (2012) Model

Our proposed approach is an extension of the Bayesian hierarchical model developed by Reich and Shaby (2012), which we review here. The Reich and Shaby (2012) model has found widespread interest for the modeling of spatial extremes because it possesses the max-stability property while being tractable in high-dimensions due to its conditional representation in terms of positive-stable variables see also Fougères et al., 2009 and Stephenson, 2009. Let and consider a set of independent

-stable random variables

, where generically the Laplace transform of has the form: . Then we construct the spatial process as the product of two independent processes,



is a white noise process (i.e., an everywhere-independent multiplicative nugget effect) with

-Fréchet marginals, , and is a spatially dependent process defined as an -norm (for ) of scaled, spatially-varying basis functions , :


The white noise process functions as a nugget effect, and accounts for measurement error occurring independently of the underlying process of interest. For small , the contribution of dominates that of the nugget effect, and vice-versa for large .

Reich and Shaby (2012) used fixed, deterministic spatial basis functions. In other words, they assumed a Dirac prior on the space of valid basis functions, based on the following construction: let be a collection of spatial knots over our spatial domain of interest , and , , be Gaussian densities centered at each knot , normalized such that for all . The Gaussian density basis functions may be replaced with normalized functions from a much broader class while still giving a valid construction for in (2). A more flexible prior for the kernels , , is discussed in Section 2.3.

The process has finite-dimensional distributions


see Tawn, 1990, which follows from the Laplace transform of an -stable variable. From (3) and the sum-to-one constraint, the marginal distributions are unit Fréchet, i.e., for all ,

Max-stability follows from (3) by checking that


The max-stability property of makes it suitable for modeling spatial extremes in scenarios of strong, non-vanishing upper tail dependence. In Section 2.2, we propose a generalized max-id model, which can better cope with weakening tail dependence.

Inference may be efficiently performed by taking advantage of the inherent hierarchical structure of the Reich and Shaby (2012) model, noticing that the data are independent conditional on the latent variables , and may be written in terms of the Fréchet distribution with scale parameter and shape parameter :


for all ; that is, , .

2.2 Sub-Asymptotic Modeling Based on a Max-Infinitely Divisible Process

Despite the appealing properties of the Reich and Shaby (2012) model, its deterministic basis functions and its max-stability make it fairly rigid in practice. Max-id processes are natural, flexible, sub-asymptotic models, that extend the class of max-stable processes while still possessing desirable properties reflecting the specific positive dependence structure of maxima. From (4), we can see that max-stable processes are always max-id. Therefore, the former form a smaller subclass within the latter.

The tail dependence class strongly determines how the probability of joint exceedances of a high threshold extrapolates to extreme quantiles. A random vector

with marginal distributions and is said to be asymptotically independent if as , and asymptotically dependent otherwise (Coles et al., 1999). We say that a spatial process { is asymptotically independent if and are asymptotically independent for all , . Max-stable processes are always asymptotically dependent (except in the case of complete independence) and, therefore, they lack flexibility to adequately capture the tail behavior of asymptotically independent data. In this section, we propose an asymptotically independent max-id model that possesses the max-stable Reich and Shaby (2012) model on the boundary of its parameter space. Dependence properties are further detailed in Section 2.5.

To extend the Reich and Shaby (2012) model to a more flexible max-id formulation, we can change the distribution of the underlying random basis coefficients . The heavy-tailedness of the distribution yields asymptotic dependence and, by construction, max-stability. To achieve asymptotic independence while staying within the class of max-id processes, we can consider a lighter-tailed, exponentially tilted, positive-stable distribution,


which was first introduced by Hougaard (1986), and has Laplace transform


Denote the density by . The density may be expressed in terms of the positive-stable density as


for , , and (Hougaard, 1986). An efficient algorithm for simulating from is given by Devroye (2009). A simple rejection sampler for the case when is not large is given in the Supplementary Material. When and , we recover the positive-stable distribution . The parameter controls the tail decay, with smaller values of corresponding to heavier-tailed distributions. Moreover, the density becomes increasingly concentrated around one as . When

, the gamma distribution with shape

and rate is obtained as .

Upon reparameterization in terms of , and , we see from (8) that is a scale parameter, which does not affect the dependence structure of our new model. Therefore, in the remainder of this paper, we set (i.e., ) and use throughout without any loss in flexibility.

When and , is an exponentially tilted form of , where the parameter has the effect of exponentially tapering the tail of at rate . Other extensions of the positive-stable distribution may also be interesting avenues for future research (e.g., polynomial tilting (Devroye, 2009)). However, our choice of (6) preserves the simplicity of the model while introducing a single parameter, the exponential tilting parameter , that is directly connected to the dependence properties of the resulting process, while allowing for inference that is computationally tractable.

Proposition 2.1.

Let be defined as in (1) with , , . Then, is max-id.


From (7), the finite-dimensional distributions for based on (6) are



the finite-dimensional distributions, denoted , from this new process satisfy
for all , and thus the process is max-id. This also confirms that the process is max-stable if and only if . ∎

Marginal distributions are no longer unit Fréchet when ; they may be expressed as


Bayesian and likelihood-based inference may be performed similarly as before, so this process enjoys the same computational benefits as the Reich and Shaby (2012) model, while having the traditional max-stable Reich and Shaby (2012) process as a special case on the boundary of the parameter space (i.e., when ). Note that unlike the Reich and Shaby (2012) model, here the marginal distributions depend on the dependence parameters and , however, this is not a problem for inference as we adopt a copula-based approach, in which we separate the treatment of the marginal distributions and the dependence structure. Marginal modeling is described in greater detail in Section 2.4. Finally, a spectral representation for the proposed max-id model is described in the Supplementary Material, which makes a link with the max-id models of Huser et al. (2018).

2.3 Prior Specification for the Spatial Kernels Based on Flexible Log-Gaussian Process Factors

The basis functions used in Reich and Shaby (2012), constructed from Gaussian densities, are radial functions, decaying symmetrically from their knot centers. While it is possible to approximate a wide range of extremal functions by considering a large collection of Gaussian density basis functions as in (2), the resulting process is overly smooth and artificially non-stationary for fixed . In this section, we propose an alternative prior for the basis functions, which allows for a parsimonious, yet flexible, stationary representation that can give insights into the predominant modes of spatial variability among of the underlying process.

More precisely, we extend the Reich and Shaby (2012)

model by replacing the Dirac prior on the Gaussian density basis functions with flexible log-Gaussian process priors, which more closely approximate the features of natural phenomena than radial basis functions. This choice of basis functions is analogous to the construction of the Brown-Resnick process

(Brown and Resnick, 1977; Kabluchko et al., 2009), which itself can be represented as the pointwise maximum over an infinite collection of scaled log-Gaussian processes. Let , , be i.i.d. mean-zero stationary Gaussian processes, each with exponential covariance function,

, whose variance and range are

and , respectively. We take the th basis to be the constant function equal to the mean of the Gaussian process, i.e., for all . Fixing the th term ensures invertibility of the , which is necessary for making posterior draws of (see Supplementary Material). Other prior choices for the basis functions that may also be worth exploring include using a more general Matérn class of covariance functions or Gaussian processes with stationary increments and an unbounded variogram (i.e., fractional Brownian motions), akin to the Brown-Resnick process. Application of a fractional Brownian motion prior in this context would require a choice of origin for each basis function, which would increase the computational cost if one wanted to marginalize over that unknown origin, and so we do not pursue it here. To satisfy the sum-to-one constraint for each spatial location , we set


The variance parameter controls the long-range spatial dependence of the max-id process , with smaller values corresponding to stronger long-range dependence (see Davison et al. (2012) for a similar discussion of geometric Gaussian processes). When is large, the difference in relative magnitudes of the unnormalized log-Gaussian processes at any given location is likely to be larger than when is small. Normalizing the basis functions when the difference in magnitudes is great gives way to more volatile fluctuations between dominating basis functions, and hence less long-range dependence. The Gaussian process range parameter governs the short-range dependence, now with larger values corresponding to stronger short-range dependence. Because the proposed basis functions provide greater flexibility in adapting to the data than the fixed Gaussian density basis, fewer basis functions are needed. In the data application presented in Section 3, we choose the number of basis functions using an out-of-sample log-score criterion. Increasing the number of basis functions allows for greater flexibility in capturing spatially dependent subregions that tend have extreme events together at the cost of greater computational burden.

When the deterministic basis functions used by Reich and Shaby (2012) are replaced with random ones, the max-stability (when ) and max-infinite divisibility properties should be interpreted conditionally on the basis functions. Both the conditional and unconditional dependence properties are described in Section 2.5.

2.4 Marginal Modeling and Realizations

For marginal distribution modeling, we use the Generalized Extreme-Value (GEV) distribution, which is the asymptotic distribution for univariate block maxima. The distribution function has the following form:

where , for some location , scale , and shape parameters, with support when , and when . Since monotone increasing transformations of the marginal distributions do not change the max-id or max-stable dependence structure, we allow for general GEV marginal distributions that are possibly different for each spatial location. In other words, we set , where is the marginal distribution of , which in the case of the Reich and Shaby (2012) model is , and in the case is given in (10), and is the quantile function for a GEV distribution with location , scale , and shape . We treat as our response. In subsequent sections, Gaussian process priors are assumed for the GEV parameters , , and

, and Markov chain Monte Carlo (MCMC) methods are used to draw posterior samples for this model. The details of the MCMC sampler are given in the Supplementary Material.

To visualize some of the features of our model, we present some sample paths in Figure 1. Realizations of on the unit square constructed using the Gaussian density (

evenly spaced basis functions, with standard deviation

) and log-Gaussian process (variance and range ) basis functions are shown in Figure 1. For illustration, the realizations have standard Gumbel margins everywhere in space, i.e., and for all . The figure illustrates the role of in controlling the relative contribution of the nugget process, and the impact of on the asymptotic dependence structure. Weaker tail dependence is present in the max-id models () than their max-stable counterparts (). Moreover, the general shapes of the Gaussian density basis model realizations appear less resemblant of natural processes than do those from the log-Gaussian process basis model.

While we have only developed the model for a single realization of the process so far, the model can easily be generalized to accommodate multiple replicates in time, which we will use in Section 3. In particular, treating time replicates of the process to be independent, we denote the maxima process observed at spatial location and time by , . We assume the marginal GEV parameters and basis functions do not vary in time, but allow the relative contribution of each basis function to be different for different time replicates of the process by taking the random basis coefficients to be , , and .

Figure 1: Realizations of the max-stable () and max-id () processes with Gaussian density (top) and log-Gaussian process (bottom) basis functions, plotted on Gumbel margins.

2.5 Dependence Properties

In this section, we explore the dependence properties of the proposed max-id model. The parameter plays a crucial role in determining the asymptotic dependence class. Reich and Shaby (2012) show that is asymptotically dependent and max-stable for , . However, when , this is no longer the case.

Proposition 2.2.

The process defined in Sections 2.22.3 using the log-Gaussian process basis prior in (11) is an asymptotically independent process when and asymptotically dependent when and .

For a proof, see Appendix A. Figure 2 displays two common dependence measures, and , (Coles et al., 1999) to illustrate the role of and in controlling the dependence properties of the tail process. Although notationally we have omitted the dependence of on and , will also depend on the locations in the (non-stationary) Gaussian density basis case. Nevertheless, while is non-stationary for the Reich and Shaby (2012) model, it is approximately stationary for a dense set of spatial knots. An attractive feature of the proposed model is that as , and transition smoothly from weak dependence to strong dependence for all .

Figure 2: Dependence measures and for the max-stable () and max-id () models for , using Gaussian density () and log-Gaussian process (, ) basis functions for and . The knots of the Gaussian density basis functions are evenly spaced between and . The figures in the bottom row correspond to after marginalizing over the log-Gaussian process basis functions based on Monte Carlo draws.

The extremal coefficient , studied by Schlather and Tawn (2003), is a measure of spatial dependence along the diagonal of the finite-dimensional distributions of max-stable processes. It takes on values from when the components are perfectly dependent to when they are independent, and therefore can be interpreted as the effective number of independent variables. The finite-dimensional distributions of a max-stable process with unit-Fréchet margins at level can be written in the form


where determines the spatial dependence and does not depend on the level . The rigidity of the dependence structure across all quantiles limits the applicability of max-stable models to processes that exhibit varying spatial dependence types at different quantiles. From (9), we can see that the max-id extension of the Reich and Shaby (2012) model does not possess this property for .

Figure 3 contrasts the spatial dependence features of the proposed models. We examine how the conditional probability of jointly exceeding a fixed quantile decays with increasing distance. Each panel shows the spatial decay of as a function of increasing spatial lag for several quantiles. We see qualitatively different behavior in the spatial decay of dependence at different quantiles between the max-stable and max-id models. In the max-stable cases, the conditional exceedance probability at short spatial lags is very similar at all levels of the distribution. The max-id models allow for more flexibility, as can be seen by the attenuated curves for higher quantiles and wider array of spatial decay types. From Figure 3, it can be seen that for , the parameter plays a role in how precipitous the decay in spatial dependence is with increasing distance, with smaller corresponding to steeper decay. Also, just as in Reich and Shaby (2012), determines the contribution of the nugget effect, which is greater when is large and lesser when is small.

Figure 3: Dependence measure between and for as a function of for max-stable (left column) and max-id (right column) models on , with Gaussian density basis functions with (top row) and log-Gaussian process basis functions with and (bottom row) basis functions for varying and . Gaussian density basis functions are evenly spaced between 0 and 1. Estimates of in the log-Gaussian process basis model are based on 50,000 Monte Carlo replicates. Horizontal dash-dot gray lines representing the values of for independent and are plotted for reference.

To confirm that our MCMC algorithm produces reliable results, and to evaluate the algorithm’s ability to infer the parameters under different regimes, we conduct a simulation study for both the Gaussian density basis and the log-Gaussian process basis models. The simulation study design and results are described in detail in the Supplementary Material. In all scenarios considered, credible intervals achieve nearly nominal levels, confirming the reliability of our MCMC algorithm.

3 Application to Extreme Precipitation

3.1 Data and Motivation

In this section, we apply our model to extreme precipitation over the northeastern United States and Canada. Our aim is twofold: (a) to understand the spatial dependence of extreme precipitation while accounting for measurement uncertainty, and (b) to predict precipitation-induced flood-risk. The data for this application were obtained from https://hdsc.nws.noaa.gov/hdsc/pfds/pfds_series.html, which is maintained by the National Oceanic and Atmospheric Administration (NOAA). Observations consist of annual maximum daily precipitation accumulations (in inches) observed between 1960 and 2015 at gauge stations (see Figure 4). The observation at gauge location , , and year , is denoted by .

Figure 4: Precipitation gauge locations () across the northeastern United States and Canada are plotted as black dots and Gaussian density basis knot locations () are plotted as red crosses.

3.2 Model Fitting and Validation

The precipitation data are analyzed by applying the four max-id models described in Section 2, namely (M1) Gaussian density basis, ; (M2) Gaussian density basis, ; (M3) log-Gaussian process basis, ; and (M4) log-Gaussian process basis, , where realizations of the process for each year are treated as i.i.d. replicates. Although further temporal dependence and trends could be modeled in both the GEV marginal parameters and basis scaling factors , Kwiatkowski-Phillips-Schmidt-Shin tests (Kwiatkowski et al., 1992) for temporal non-stationarity among the annual maxima were performed separately for each station, and 85% of stations yielded no evidence for temporal non-stationarity at confidence level . The proposed model would be more complex and computationally demanding to fit if one were to account for temporal non-stationarity. Therefore, for the sake of simplicity, and since overall the data do not appear to be highly non-stationary over time, we will ignore this aspect in our analysis. Accounting for temporal non-stationarity would be an interesting avenue of future research to further develop this model.

In particular, both the dependence model and GEV marginal distributions are assumed to be constant over time. We assume independent Gaussian process priors, each with constant mean and stationary exponential covariance function , , on the location and log-scale marginal parameters of the GEV distribution, with half-normal priors for and . Due to the difficulty in estimating the shape parameter (Cooley et al., 2007; Opitz et al., 2018), we use a spatially constant prior, . The dependence parameter priors are as follows: For and , we take and . For the Gaussian density basis models, we use knot locations on an evenly spaced grid (see Figure 4). A half normal prior is put on the Gaussian density bandwidth parameter . In the case of the log-Gaussian process basis models, we consider and basis functions, putting priors and

on the exponential covariance parameters. Handling missing values is straightforward using the proposed approach. For each iteration of the MCMC algorithm, missing values are sampled from the posterior predictive distribution; this is detailed in the Supplementary Material. We run each MCMC chain under two different parameter initializations for 40,000 iterations using a burn-in of 10,000 with data from

stations, reserving stations for model evaluation. In all four cases, the posterior densities were similar across the two initializations.

It is currently not possible to fit existing max-stable, inverted-max-stable (Wadsworth and Tawn, 2012), and other max-id models see, e.g. Huser et al., 2018, Padoan, 2013 using a full likelihood or Bayesian approach when the number of spatial locations is large; see Castruccio et al. (2016), Dombry et al. (2017) and Huser et al. (2019). Under these constraints, a natural alternative for comparison is the model for block maxima proposed by Sang and Gelfand (2010), which also belongs to the asymptotic independence class. Specifically, let be a mean-zero Gaussian process with exponential correlation function and unit variance. The annual maxima are then modeled as , where the location and each follow mean zero Gaussian processes with exponential covariance functions, with the same priors as above, and

denotes the standard normal distribution function. We refer to this as the the GEV-Gaussian process copula model.

To compare models, we calculate out-of-sample log-scores (Gneiting and Raftery, 2007), for annual maxima at the 100 holdout stations, which is simply the log-likelihood of the holdout data for each model based on conditional predictive simulations of the latent model parameters at the unobserved sites. Since the log-scores are calculated on holdout data, they implicitly account for model complexity. We also emphasize that because the predictions are based on the joint likelihood, the log-scores reflect not only the marginal fits, but also how well the model captures the dependence characteristics of the observed data. The best log-score (higher scores are better) of the two initializations for each model is reported in Table 1. The max-id models () outperform their max-stable counterparts (). The log-score for the GEV-Gaussian process copula model is worse than the other models considered. The estimated marginal surfaces are similar across all of the models considered, indicating that the misspecification is due to differences in the dependence model for the annual maxima.

The max-id, log-Gaussian process basis model with and basis functions has the highest log-score (shown in bold), suggesting it should be preferred among the considered models for this data application, and as such we focus on this model for the remainder of our analysis. For this model, the posterior mean (95% credible interval) estimates of the dependence parameters are for , for , and for the spatial basis functions for and miles for , suggesting the presence of some residual dependence beyond that explained by spatially-varying marginal parameters. Also, while we have specified vague priors on the model parameters, the posterior distributions are highly concentrated around their corresponding posterior means. Although the proposed inference scheme does not allow for jumps between and , the posterior samples of are still somewhat informative about the asymptotic dependence class. In particular, since the dependence properties of our model are smooth in at zero, the fact that the 95% credible interval for is relatively symmetric and distant from 0 gives support for asymptotic independence among precipitation extremes.

To validate the decision of having the same dependence parameters and over the entire region, log-Gaussian process basis models with were also separately fitted to four subregions, two inland and two coastal. The 95% credible intervals for and overlap with those fitted to the entire region, suggesting homogeneous spatial dependence of the process over the study region.

. Higher log-scores correspond to better fit. The max-id, log-Gaussian process basis model has the highest log-score (shown in bold). Gaussian Density Basis log-Gaussian Process Basis GEV-Gaussian Process Copula L 60 10 15 20 -5292.5 -5410.7 -5406.4 -5415.2 -6097.048 -5218.3 -5194.6 -5172.6 -5207.9

Table 1: Log-scores estimated from annual maxima observed at the holdout stations are used to compare the four models presented in Section 2, and the GEV-Gaussian process copula model

Further, to examine the model fit, we compare empirical and model-based estimates of as a function of spatial lag and threshold for the holdout stations (Figure 5).

Figure 5: The left panel shows as a function of for fixed spatial lags miles calculated for the 100-holdout stations. Empirical estimates are shown as a solid black line, and max-id, log-Gaussian process basis model 95% credible intervals are shown as gray ribbons. The decay of towards zero as suggests that daily precipitation are asymptotically independent. To understand the spatial dependence of extreme precipitation at increasingly extreme levels, empirical (solid lines) and model 95% credible intervals (ribbons) of for the holdout stations are plotted for several quantiles (right panel). Horizontal dash-dot gray lines representing the values of under an everywhere-independent model are plotted for reference. The plot shows good overall agreement between the model fits and empirical estimates, except at very short distances.

The left panel shows as a function of for at fixed lags miles, and the right panel shows the spatial decay of as a function of spatial lag for several fixed marginal quantiles . Empirical estimates are represented by solid lines and credible intervals for each model by shaded ribbons. From the left panel, we can see that the max-id model captures the asymptotic independence behavior of the precipitation data quite well, except at very short distances where our model slightly underestimates the strong dependence present in the data. The max-stable model also underestimates the relatively strong dependence at shorter distances, but with comparable coverage to the max-id model at other distances (see Supplementary Material). The slight discrepancy at shorter distances may be due to the phenomenon described by Robins et al. (2000) wherein intervals from posterior summaries like that are calculated from MCMC draws are too narrow. From the right panel, we deduce that the annual maximum precipitation data exhibit quite strong spatial dependence up to about miles, with weaker spatial dependence at higher quantiles. Moreover, decays towards its independence level as a function of distance faster at the and quantiles than at the and quantiles.

In order to assess the spatial prediction skill of our model, we display in Figure 6 quantile-quantile (QQ)-plots for group-wise summaries of the annual maxima taken over the holdout stations (see Davison et al. (2012) for a similar analysis). The results show adequate correspondence between the model-based and empirical quantiles of the group-wise means, whereas the observed group-wise minima (maxima, respectively) appear to be slightly underestimated (overestimated, respectively) by the model. Corresponding QQ-plots when (not shown) give similar patterns with minima (maxima, respectively) lying slightly further above (below, respectively) the 95% credible intervals.

Figure 6: QQ-plots of the observed and predicted group-wise minima (left), mean (center), and maxima (right) taken over the annual maxima from all 100 holdout stations. The dashed lines represent 95% credible intervals. The plots reflect reasonable correspondence between the empirical and modeled multivariate distributions. To account for the fact that the marginal GEV distributions vary across stations, observations are first transformed to unit Gumbel scale using the probability integral transform for the GEV marginal distributions at each station from the fitted model.

Maps of the marginal posterior predictive means and standard deviations of the 0.99 quantile of annual maxima (i.e., 100-year return level) for the max-id, log-Gaussian process basis model are shown in Figure 7. The posterior mean surfaces are consistent with marginal quantile surfaces for the region as reported in NOAA Atlas 14 (Perica et al., 2013). The posterior standard deviation surface shows the greatest variability in Maine, Long Island, and along the boundary of the observation region where there are relatively few gauge locations. For illustration, observed maxima in 2012, a single posterior predictive draw, and the posterior predictive mean for that year are plotted in Figure 8. The posterior predictive plots appear to capture the general spatial characteristics of the maxima observed in 2012 well.

Figure 7: Pointwise posterior predictive mean (left) and standard deviation (right) of the 100-year return level of daily precipitation.
Figure 8: Observed precipitation accumulations (left), a single posterior predictive draw (middle), and posterior predictive means (right) for the year 2012. Missing values are shown in gray.

3.3 Principal Modes of Spatial Variability Among Precipitation Extremes

Spatial principal component analysis (PCA)

(Demsar et al., 2013; Jolliffe, 2002) and Empirical Orthogonal Functions (Hannachi et al., 2007) have proven to be useful methods for exploring the main large scale features of spatial processes. However, aside from recent work by Morris (2016) and Cooley and Thibaud (2018), little has been done to this end for spatial extremes. The model we have proposed allows for an exploratory visualization that is very similar to a spatial PCA method that Demsar et al. (2013) refers to as Atmospheric Science PCA in their review of Spatial PCA methods, where the data consist of time replicates of a univariate spatial process observed at several locations.

An attractive feature of the log-Gaussian process basis model is that it provides a low-dimensional representation of the predominant modes of spatial variability among extremes. Analogously to factor analysis, the primary spatial trends among extreme precipitation can be described by a subset of the spatial basis functions that contribute the most to the overall process. To achieve this, motivated by PCA factorization, which finds the directions of maximum variance in the data, we rank the spatial basis functions , by the posterior year-to-year variation of their corresponding basis coefficients (i.e., higher posterior variance corresponds to lower rank). Arguably, both the means and variances of the coefficients play a role in the relative contribution of the corresponding basis function to the overall process. However, from inspection, the basis coefficients with the highest posterior variance also have the highest posterior means. Examining the variance of the basis coefficients for each , against their ranks give a rough indication of the number of basis functions with sizable contributions to the overall process. Also, while label switching is possible, from inspection of the MCMC samples of the basis functions, this does not appear to be a major concern for this application. If label switching is present, application of the pivotal reordering algorithm proposed by Marin et al. (2005); Marin and Robert (2007) can be used to permute the labels of the basis functions and scaling factors before ranking the basis functions. Posterior means of the first six spatial basis functions are shown in Figure 9. Most of the top ranked factor means in the basis function case were also identified as top ranked functions in the and case (see Supplementary Material).

Figure 9: First six spatial basis functions ordered by the variance of their corresponding random basis coefficients from largest to smallest (left to right, top to bottom) for the basis function model. The year-to-year variation among the coefficients of these first six basis functions accounts for 97% of the total year-to-year variation among all of the basis coefficients. The shapes of the latent factors have reasonable interpretations in terms of geographic coastal and mountain features.

Unlike the pointwise marginal surfaces, which do not provide any information about the joint dependence of extremes, these basis functions capture spatial regions of simultaneous (in this case, merely the same year) extreme precipitation. The proportion of the total variation among the accounted for by variation in the coefficients of each of the first six basis functions is and respectively. This does not imply that the top ranked factor is the dominating kernel of the time. Rather, if the variance of the scaling coefficients for the th factor is high, then the year-to-year differences in the spatial modes of extremes should be well described by the peaks and troughs of the th factor. For example, if has a peak around some location then the conditional GEV distribution (given the factors and scaling coefficients) will be stochastically larger at in years when is large and smaller when is small. Therefore, the low ranked factors describe regions where precipitation tends to be extreme together or more moderate together. The latent factors in Figure 9 have reasonable physical interpretations that are reflective of natural geographic features. In particular, they resemble observed patterns in extreme precipitation events occurring along the coast and mountain range borders.

3.4 Drainage Basin Flood Risk Analysis

To understand flood-risk, it is necessary to account for the Earth’s topographic features that dictate the flow of rainwater. Drainage basins, also commonly referred to as catchment areas or watersheds, delineate a spatial region into hydrological units based on the natural flow of water. Precipitation falling over a single drainage basin collects into a common outlet (e.g., river, lake, or bay).

The United States is divided into increasingly granular hydrological units, each identified by a unique hydrological unit code (HUC). There are six levels to the hierarchy: regions, sub-regions, basins, sub-basins, watersheds, and sub-watersheds, each level of the hierarchy possessing a two digit code. To reconcile the analysis of componentwise maxima—which for a given year, may occur on different days at different stations—with the aggregate quantities needed for a flood analysis, we consider small drainage areas. Within a sufficiently small area, it is reasonable to assume that the pointwise maximum daily precipitation event for a single year occurs on the same day at every point, making it possible to calculate the volume of precipitation due to a single maximum precipitation event for the entire area. With this goal in mind, and in light of the long-range dependence of extreme precipitation in this region (recall Figure 5), we consider the basin (HUC6) and sub-basin (HUC8) resolutions. A majority of stations in the sub-basin had their annual maximum event on the same day of the year in over 60% of the sub-basin-year pairs during the observation period. However, it is also common for the annual maxima among stations in a single sub-basin to occur on subsequent days. In the analysis that follows, our forecasts will treat annual maximum events as occurring on the same day over a given sub-basin, and are therefore conservative.

The basins we consider are motivated by historical floods. In 1938, the Great New England Hurricane, one of the most powerful in recorded history, made landfall in southern New England, dropping over six inches of rainfall in some areas, and causing at that time what was the greatest amount property damage that had ever occurred due to a single storm (Brooks, 1939). For this reason, we consider Merrimack River basin covering northern Massachusetts and central New Hampshire. In addition, we also consider the Susquehanna River basins in southern New York and northeastern Pennsylvania, some of the most flood-prone areas in the country. The left panel of Figure 10 illustrates the watershed boundaries of the chosen basins.

Figure 10: Merrimack River (red) and Susquehanna River (blue) basins and estimated survivor curves of the total daily precipitation volume. The middle column shows area averaged, pointwise estimates, which are calculated by dividing the total predicted basin precipitation volume by the basin area.

To estimate the total volume of precipitation over a basin, we make posterior predictive draws over each sub-basin on a grid of 1 cells by applying the prediction at the cell center uniformly over the entire cell. Specifically, denoting the prediction cell centers by , each with area , we repeat the following procedure 1,000 times for each MCMC iteration to estimate the quantiles of total precipitation volume for each basin:

1:With equal probability, randomly select one of the observation years, .
2:Draw for , from the posterior predictive distribution.
3:Estimate the total volume of precipitation over the basin by .

Estimated survivor curves of the annual maximum total daily precipitation volume are shown in Figure 10. The right panels show estimates the total precipitation volume survivor curves, and the middle panels show the total volume estimates normalized by the total area of each sub-basin to give a measure of the flood risk relative to the size of the basin. The information from this analysis gives a sense of the flood risk due to annual precipitation events occurring on a single day. An infrastructure planner can use such estimates of the -year return level of maximum precipitation volume (the level exceeded on average once every years or, equivalently, the th quantile of annual maximum precipitation volume) directly to incorporate flood risk into their design requirements. We see that the Merrimack River basin, while having a stochastically smaller estimated distribution for total precipitation volume has a stochastically larger distribution once it is normalized by basin size. We perform the same analysis for all of the sub-basins of the Merrimack River and Susquehanna River basins (Figure 11), where the assumption of simultaneous annual maxima across the entire area is even more convincing. The sub-basin survivor curves give a more granular assessment of the flood risk of the local tributaries of the main rivers. The estimated 10, 30, and 50-year return-levels for the area averaged precipitation (in/mi) and total precipitation volume (ft) are reported in Table 2.

Figure 11: Merrimack River (top-left) and Susquehanna River (bottom-left) sub-basins and estimated survivor curves of the total daily precipitation volume. The middle column shows area averaged, pointwise estimates, which are calculated by dividing the total predicted sub-basin precipitation volume by the sub-basin area. The first two digits of the eight digit Hydrological Unit Code (HUC) identify the region (top: New England (01) and bottom: Mid-Atlantic (02)) and have been omitted for readability.
Basin/Sub-basin Area Avg. Precipitation (in/mi) Precipitation Volume (billion ft)
HUC 6     10-year 30-year 50-year     10-year 30-year 50-year
010700     3.52 4.12 4.35     40.85 47.82 50.53
020501     2.52 2.81 2.92     66.35 73.98 76.74
HUC 8     10-year 30-year 50-year     10-year 30-year 50-year
01070001     3.17 3.70 3.92     7.62 8.89 9.42
01070002     3.54 4.26 4.46     4.03 4.86 5.09
01070003     3.39 4.06 4.31     6.08 7.27 7.72
01070004     3.44 4.06 4.29     4.20 4.96 5.24
01070005     3.60 4.28 4.52     3.34 3.96 4.19
01070006     3.77 4.32 4.56     15.58 17.87 18.87

    2.49 2.71 2.79     13.35 14.51 14.95
02050102     2.41 2.67 2.75     9.10 10.10 10.41
02050103     2.38 2.68 2.79     5.80 6.54 6.80
02050104     2.43 2.74 2.88     7.88 8.87 9.33
02050105     2.39 2.69 2.82     6.71 7.55 7.93
02050106     2.52 2.87 2.96     11.65 13.25 13.69
02050107     2.91 3.24 3.35     11.86 13.17 13.63

Table 2: Area averaged precipitation (in/mi) and total precipitation volume (billion ft) return-level estimates for 10, 30, and 50-year return-periods are reported for the Merrimack and Susquehanna basins and sub-basins.

4 Discussion

In this paper, we extend the max-stable model for spatial extremes developed by Reich and Shaby (2012) in several ways. First, by using flexible log-Gaussian process basis functions, our model provides a more realistic low-dimensional factor representation that can be used to visualize the main modes of spatial variability among extremes. Second, our approach relaxes the rigid spatial dependence structure imposed by max-stable models, while possessing the positive dependence inherent to distributions for maxima. Inference on the tail dependence class is also possible, as our model can capture asymptotic independence when , while having an asymptotically dependent, max-stable model on the boundary of the parameter space (when ).

We apply our model to extreme precipitation over the northeastern United States and Canada. Because it accounts for the spatial dependence among maxima and we are able to efficiently make conditional draws from our fitted model, it is possible to estimate total precipitation volume survivor curves for annual maxima over hydrologically defined sub-basins. The precipitation predictions from our model could be incorporated into a hydrological model for the flow path dynamics that incorporates factors like drainage basin topography, land use, and land cover to describe how precipitation falling over a common catchment translates into drainage and potential flooding. The precipitation analysis does not account for the cumulative effect of heavy precipitation over several days, which can overload an urban stormwater drainage system that is already operating at capacity. Further temporal modeling of the marginal distributions and space-time dependence characteristics would facilitate such an analysis; see, e.g., Huser and Davison (2014) for space-time modeling of precipitation extremes using max-stable processes.

For future work, adding a point mass at in the prior and proposal distributions would make it possible to account for model uncertainty and simultaneously perform model selection directly within the MCMC. Moreover, an interesting extension of the model described here might be to assume heavier-tailed random scaling factors, such that the resulting process is asymptotically dependent, while possibly losing the max-id property. Finally, while our focus in this paper has been on flexible sub-asymptotic modeling of maxima, another avenue for research is to investigate relaxing the rigid dependence structure of limiting generalized Pareto process models for peaks-over-threshold data (see, e.g., Castro Camilo and Huser, 2018; Huser and Wadsworth, 2019).

Appendix A Model Tail Dependence Properties

Since the marginal distributions of are the same when constructed using the log-Gaussian process basis, and are asymptotically independent if as . The marginal distribution of the process at location conditional on the basis functions is , and the joint distribution at two locations and is

For brevity, we will drop the indices , and write, e.g., . By ’s rule, we obtain

when . Finally, by application of the Dominated Convergence Theorem, since , we obtain for all For more detail, see the Supplementary Material.

In the case of and , Reich and Shaby (2012) showed that is max-stable with extremal coefficient . Using the relation for max-stable processes with unit Fréchet margins that , and by the Dominated Convergence Theorem, we have when for all . So, when and , is asymptotically dependent, both conditionally on and unconditionally.


  • Balkema and Resnick (1977) Balkema, A. A. and S. I. Resnick (1977). Max-infinite divisibility. Journal of Applied Probability 14(2), 309–319.
  • Brooks (1939) Brooks, C. F. (1939). Hurricanes into new england: meteorology of the storm of september 21, 1938. Geographical Review 29(1), 119–127.
  • Brown and Resnick (1977) Brown, B. M. and S. I. Resnick (1977). Extreme values of independent stochastic processes. Journal of Applied Probability 14(4), 732–739.
  • Castro Camilo and Huser (2018) Castro Camilo, D. and R. Huser (2018). Local likelihood estimation of complex tail dependence structures, applied to U.S. precipitation extremes. arXiv preprint 1710.00875.
  • Castruccio et al. (2016) Castruccio, S., R. Huser, and M. G. Genton (2016). High-order composite likelihood inference for max-stable distributions and processes. Journal of Computational and Graphical Statistics 25(4), 1212–1229.
  • Coles et al. (1999) Coles, S., J. Heffernan, and J. Tawn (1999). Dependence measures for extreme value analyses. Extremes 2(4), 339.
  • Cooley et al. (2007) Cooley, D., D. Nychka, and P. Naveau (2007). Bayesian spatial modeling of extreme precipitation return levels. Journal of the American Statistical Association 102(479), 824–840.
  • Cooley and Thibaud (2018) Cooley, D. and E. Thibaud (2018). Decompositions of dependence for high-dimensional extremes. arXiv preprint 1612.07190.
  • Davison and Huser (2015) Davison, A. C. and R. Huser (2015). Statistics of extremes. Annual Review of Statistics and Its Application 2, 203–235.
  • Davison et al. (2013) Davison, A. C., R. Huser, and E. Thibaud (2013). Geostatistics of dependent and asymptotically independent extremes. Mathematical Geosciences 45(5), 511–529.
  • Davison et al. (2019) Davison, A. C., R. Huser, and E. Thibaud (2019). Spatial extremes. In A. E. Gelfand, M. Fuentes, and R. L. Smith (Eds.), Handbook of Environmental and Ecological Statistics. CRC Press. (forthcoming).
  • Davison et al. (2012) Davison, A. C., S. A. Padoan, and M. Ribatet (2012). Statistical modeling of spatial extremes. Statistical Science. A Review Journal of the Institute of Mathematical Statistics 27(2), 161–186.
  • de Haan and Ferreira (2006) de Haan, L. and A. Ferreira (2006). Extreme value theory. Springer Series in Operations Research and Financial Engineering. Springer, New York. An introduction.
  • Demsar et al. (2013) Demsar, U., P. Harris, C. Brunsdon, A. S. Fotheringham, and S. McLoone (2013). Principal component analysis on spatial data: an overview. Annals of the Association of American Geographers 103(1), 106–128.
  • Devroye (2009) Devroye, L. (2009). Random variate generation for exponentially and polynomially tilted stable distributions. ACM Transactions on Modeling and Computer Simulation (TOMACS) 19(4), 18.
  • Dombry et al. (2017) Dombry, C., S. Engelke, and M. Oesting (2017). Bayesian inference for multivariate extreme value distributions. Electronic Journal of Statistics 11(2), 4813–4844.
  • Dombry et al. (2013) Dombry, C., F. Éyi-Minko, and M. Ribatet (2013). Conditional simulation of max-stable processes. Biometrika 100(1), 111–124.
  • Ferreira and de Haan (2014) Ferreira, A. and L. de Haan (2014). The generalized Pareto process; with a view towards application and simulation. Bernoulli. Official Journal of the Bernoulli Society for Mathematical Statistics and Probability 20(4), 1717–1737.
  • Fougères et al. (2009) Fougères, A.-L., J. P. Nolan, and H. Rootzén (2009). Models for dependent extremes using stable mixtures. Scandinavian Journal of Statistics. Theory and Applications 36(1), 42–59.
  • Gneiting and Raftery (2007) Gneiting, T. and A. E. Raftery (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association 102(477), 359–378.
  • Hannachi et al. (2007) Hannachi, A., I. Jolliffe, and D. Stephenson (2007). Empirical orthogonal functions and related techniques in atmospheric science: A review. International journal of climatology 27(9), 1119–1152.
  • Hougaard (1986) Hougaard, P. (1986). Survival models for heterogeneous populations derived from stable distributions. Biometrika 73(2), 387–396.
  • Huser and Davison (2014) Huser, R. and A. C. Davison (2014). Space-time modelling of extreme events. Journal of the Royal Statistical Society. Series B. Statistical Methodology 76(2), 439–461.
  • Huser et al. (2019) Huser, R., C. Dombry, M. Ribatet, and M. G. Genton (2019). Full likelihood inference for max-stable data. Stat 8, e218.
  • Huser et al. (2017) Huser, R., T. Opitz, and E. Thibaud (2017). Bridging asymptotic independence and dependence in spatial extremes using Gaussian scale mixtures. Spatial Statistics 21(part A), 166–186.
  • Huser et al. (2018) Huser, R., T. Opitz, and E. Thibaud (2018). Max-infinitely divisible models and inference for spatial extremes. arXiv preprint 1801.02946.
  • Huser and Wadsworth (2019) Huser, R. and J. L. Wadsworth (2019). Modeling spatial processes with unknown extremal dependence class. Journal of the American Statistical Association. To appear.
  • Jolliffe (2002) Jolliffe, I. T. (2002). Principal component analysis (Second ed.). Springer Series in Statistics. Springer-Verlag, New York.
  • Kabluchko and Schlather (2010) Kabluchko, Z. and M. Schlather (2010). Ergodic properties of max-infinitely divisible processes. Stochastic Processes and their Applications 120(3), 281–295.
  • Kabluchko et al. (2009) Kabluchko, Z., M. Schlather, and L. de Haan (2009). Stationary max-stable fields associated to negative definite functions. The Annals of Probability 37(5), 2042–2065.
  • Kabluchko and Stoev (2016) Kabluchko, Z. and S. Stoev (2016). Stochastic integral representations and classification of sum- and max-infinitely divisible processes. Bernoulli. Official Journal of the Bernoulli Society for Mathematical Statistics and Probability 22(1), 107–142.
  • Kwiatkowski et al. (1992) Kwiatkowski, D., P. C. Phillips, P. Schmidt, and Y. Shin (1992).

    Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are we that economic time series have a unit root?

    Journal of econometrics 54(1-3), 159–178.
  • Marin et al. (2005) Marin, J.-M., K. Mengersen, and C. P. Robert (2005). Bayesian modelling and inference on mixtures of distributions. In Bayesian thinking: modeling and computation, Volume 25 of Handbook of Statistics, pp. 459–507. Elsevier/North-Holland, Amsterdam.
  • Marin and Robert (2007) Marin, J.-M. and C. P. Robert (2007).

    Bayesian core: a practical approach to computational Bayesian statistics

    Springer Texts in Statistics. Springer, New York.
  • Morris (2016) Morris, S. A. (2016). Spatial Methods for Modeling Extreme and Rare Events. Thesis (Ph.D.) North Carolina State University.
  • Opitz et al. (2018) Opitz, T., R. Huser, H. Bakka, and H. Rue (2018). INLA goes extreme: Bayesian tail regression for the estimation of high spatio-temporal quantiles. Extremes 21(3), 441–462.
  • Padoan (2013) Padoan, S. A. (2013). Extreme dependence models based on event magnitude.

    Journal of Multivariate Analysis

     122, 1–19.
  • Perica et al. (2013) Perica, S., D. Martin, S. Pavlovic, I. Roy, M. S. Laurent, C. Trypaluk, D. Unruh, M. Yekta, and G. Bonnin (2013). Noaa atlas 14 volume 9 version 2, precipitation-frequency atlas of the united states, southeastern states. NOAA, National Weather Service 9, 18.
  • Reich and Shaby (2012) Reich, B. J. and B. A. Shaby (2012). A hierarchical max-stable spatial model for extreme precipitation. The Annals of Applied Statistics 6(4), 1430–1451.
  • Resnick (1987) Resnick, S. I. (1987). Extreme values, regular variation, and point processes, Volume 4 of Applied Probability. A Series of the Applied Probability Trust. Springer-Verlag, New York.
  • Robins et al. (2000) Robins, J. M., A. van der Vaart, and V. Ventura (2000). Asymptotic distribution of values in composite null models. J. Amer. Statist. Assoc. 95(452), 1143–1167, 1171–1172. With comments and a rejoinder by the authors.
  • Sang and Gelfand (2010) Sang, H. and A. E. Gelfand (2010). Continuous spatial process models for spatial extreme values. Journal of Agricultural, Biological, and Environmental Statistics 15(1), 49–65.
  • Schlather and Tawn (2003) Schlather, M. and J. A. Tawn (2003). A dependence measure for multivariate and spatial extreme values: properties and inference. Biometrika 90(1), 139–156.
  • Stephenson (2009) Stephenson, A. G. (2009).

    High-dimensional parametric modelling of multivariate extreme events.

    Australian & New Zealand Journal of Statistics 51(1), 77–88.
  • Suro et al. (2009) Suro, T. P., G. D. Firda, and C. O. Szabo (2009). Flood of june 26-29, 2006, mohawk, delaware, and susquehanna river basins, new york. US Geological Survey Open-File Report 1063, 354.
  • Tawn (1990) Tawn, J. A. (1990). Modelling multivariate extreme value distributions. Biometrika 77(2), 245–253.
  • Thibaud and Opitz (2015) Thibaud, E. and T. Opitz (2015). Efficient inference and simulation for elliptical Pareto processes. Biometrika 102(4), 855–870.
  • Wadsworth and Tawn (2012) Wadsworth, J. L. and J. A. Tawn (2012). Dependence modelling for spatial extremes. Biometrika 99(2), 253–272.