Spatial Modeling of Heavy Precipitation by Coupling Weather Station Recordings and Ensemble Forecasts with Max-Stable Processes

03/12/2020 ∙ by Marco Oesting, et al. ∙ 0

Due to complex physical phenomena, the distribution of heavy rainfall events is difficult to model spatially. Physically based numerical models can often provide physically coherent spatial patterns, but may miss some important precipitation features like heavy rainfall intensities. Measurements at ground-based weather stations, however, supply adequate rainfall intensities, but most national weather recording networks are often spatially too sparse to capture rainfall patterns adequately. To bring the best out of these two sources of information, climatologists and hydrologists have been seeking models that can efficiently merge different types of rainfall data. One inherent difficulty is to capture the appropriate multivariate dependence structure among rainfall maxima. For this purpose, multivariate extreme value theory suggests the use of a max-stable process. Such a process can be represented by a max-linear combination of independent copies of a hidden stochastic process weighted by a Poisson point process. In practice, the choice of this hidden process is non-trivial, especially if anisotropy, non-stationarity and nugget effects are present in the spatial data at hand. By coupling forecast ensemble data from the French national weather service (Météo-France) with local observations, we construct and compare different types of data driven max-stable processes that are parsimonious in parameters, easy to simulate and capable of reproducing nugget effects and spatial non-stationarities. We also compare our new method with classical approaches from spatial extreme value theory such as Brown-Resnick processes.



There are no comments yet.


page 6

page 10

page 12

page 15

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

A recent review report on weather and climate extremes Chen et al. (2018) has emphasized the need of “theoretical frameworks and statistical methods” for modeling complex extreme events. This brings the particular question of how to statistically model the dependence among extremes. The field of multivariate extreme value theory (EVT), see the books by Resnick (2008); Embrechts, Klüppelberg and Mikosch (1997); De Haan and Ferreira (2006), for instance, offers a particular mathematical framework to address such a question.

In this work, we deal with modelling extreme precipitation. There exists a series of articles that have developed parametric EVT models to capture the multivariate extremal dependence structure of daily or subdaily precipitation (see, e.g. Cooley, Nychka and Naveau, 2007; Keef, Svensson and Tawn, 2009; Bechler, Bel and Vrac, 2015; Buhl and Klüppelberg, 2016; Saunders et al., 2017; De Fondeville and Davison, 2018). Keeping in mind that heavy rainfall events are difficult to model due to their large spatial and temporal variability (see, e.g. Sillmann et al., 2017)

, the choice of a parametric model that adequately characterizes rainfall heterogeneity is a non-trivial problem. For example, heavy rainfall events in the south of France are modulated by various effects due to the French orography, see the bottom panel in Figure

1, and different types of atmospheric patterns (see Carreau et al., 2012; Carreau, Naveau and Neppel, 2017; Blanchet, Molinié and Touati, 2018, for instance). These effects raise the question if relevant spatial information about rainfall patterns can be found to complement observations measured by weather stations. This inquiry moves the focus from finding and fitting complex parametric multivariate EVT structures to the problem of coupling local precipitation records with data at different spatial scales.

Nowadays, numerical weather models based on assimilation schemes are able to provide accurate ensemble forecasts of various atmospheric variables, in particular temperature fields (see, e.g. Taillardat et al., 2016)

. Due to their large variability, heavy rainfall intensities remain, even in terms of their probability distributions, a challenge for weather forecasts. In particular, forecasted rainfall extremes may strongly differ from precipitation records measured at weather stations. The later provide high quality data at the local spatial scale (the weather station), but high quality and well maintained observational networks have a spatial resolution which is much worse than current ensemble forecasts. To illustrate this discrepancy of spatial scales between station networks and weather numerical models, the top panel of Figure

1 superimposes the well monitored 110 Météo France weather station locations on top of the PEARP grid (see Taillardat et al., 2016, for instance) which is used for the numerical weather prediction system operated by Météo France.

Despite the challenges in modeling heavy precipitation intensities accurately, ensemble forecasts of rainfall data still provide relevant information in terms of spatial rainfall patterns, (see Taillardat et al., 2019, for instance). Compared to climate models, numerical weather models contain fine scale features and have complex parametrizations, and throughout their assimilation schemes, they spatially track storm patterns. This makes them good candidates as proxies of realistic rainfall patterns although their intensities can be misrepresented.

In this context, our main question is how to couple ensembles of rainfall forecasts within the construction of models suggested by multivariate EVT in order to simulate coherent spatial fields of extreme precipitation, while preserving the spatial structure observed in weather stations. As a direct by-product, any solution of our research problem is supposed to provide extreme precipitation fields at a very fine scale, the one of the ensemble forecast, see Figure 1 (top). To reach these objectives, Section 2 provides some theoretical background on the underlying statistical models used in our further analysis, i.e. max-stable models from EVT. Our two available precipitation datasets, observational records at weather stations in the south of France and forecast ensembles on a grid, are presented in Section 3. In Section 4, we present different ways to couple both types of data within a max-stable framework. More precisely, four data driven models are introduced, fitted to the region of interest, see Figure 1, and also compared to a classical stationary spatial max-stable model, the Brown–Resnick process. The paper closes with a discussion in Section 5.

2 Max-stable models

In this section, we will provide the theoretical background on the models we will use to model heavy precipitation events.

We start with one of the main results from univariate extreme value theory, the Fisher–Tippett Theorem (Fisher and Tippett, 1928)

. Assume that we have independently and identically distributed random variables

, e.g. different precipitation measurements at the same station or different forecasts for the same grid cell. If there exist normalizing sequences and , , such that the normalized maximum converges in distribution, i.e.

for some non-degenerate limit distribution , then is necessarily a Generalized Extreme Value (GEV) distribution

for a location , a scale and a shape parameter . Therefore, the distribution of maxima over certain time periods (so-called blocks) at a single station or grid cell are typically modeled by a GEV distribution. Note that the result above still holds true if the observations are not independent, but satisfy certain mixing conditions (see Leadbetter, Lindgren and Rootzén, 1983).

For modeling extreme events in space, we need an extension of the above limit result to stochastic processes: Let be independent copies of a stochastic process on some countable index set , e.g. a dense grid. Then, provided that there exist sequences of normalizing functions and , , , such that the process of normalized pointwise maxima

converges in distribution to a stochastic process with non-degenerate marginal distributions, this process is necessarily max-stable. From the univariate result, we immediately obtain that follows a GEV distribution , .

In this paper, we focus on the extremal spatial dependence structure. Therefore, by appropriate marginal transformations, we assume that possesses unit Fréchet margins, i.e.  for all . Such a process is called simple max-stable. By De Haan (1984), every simple max-stable process can be represented as


where are the points of a Poisson point process on with intensity and, independently from the Poisson points, , , are independent copies of a stochastic process such that for all .

As the intensity of the Poisson point process is fixed for given marginal distributions, the spatial dependence structure is fully determined by the multivariate distributions of the so-called spectral process . Many classes of max-stable models are given by specific choices of . For instance, a process of the form

for some centered Gaussian process , leads to a so-called Brown–Resnick process – one of the most popular max-stable models. If is a grid and has stationary increments, the corresponding Brown–Resnick process is stationary and its law depends on the variogram


only. Therefore, is also called Brown–Resnick process associated to variogram (Brown and Resnick, 1977; Kabluchko, Schlather and de Haan, 2009).

If the spectral process follows a discrete uniform distribution on some finite set of nonnegative functions

on with for all , representation (1) simplifies to the max-linear model (Wang and Stoev, 2011)


where are independent unit Fréchet random variables.

As an alternative to the max-linear model (3) which is given by the maximum over a finite number of basis functions , Reich and Shaby (2012) developed a max-stable model that can be written as a sum of , see also Oesting (2018) for a further generalization:


where is a noise process with -Fréchet marginal distributions and are i.i.d. -stable random variables whose distribution is given by the Laplace transform

for some . Compared to the max-linear model (3), the basis functions in the Reich and Shaby (2012) allow a multiplicative random nugget effect in model (4). This nugget effect and the simple additive form based on a finite sum make this model attractive as a spatial model in environmental applications.

A popular dependence measure for a max-stable process with marginal distribution is the extremal coefficient defined via


Note that, here, does not depend on and, if a.s., then , while corresponds to the independence case. For a simple max-stable process with general representation (1), the extremal coefficient can also be expressed as

This allows us to make the link with max-stable models studied in this work,


denotes the standard normal distribution function.

In classical geostatistics, spatial dependence is often summarized via variograms which correspond to -distances, see Equation (2

). As both variance and expectation are non-finite in case of unit Fréchet margins, other distances have to be used in such cases. For example,

Cooley, Naveau and Poncet (2006) studied a marginal free distance, called a -madogram,

To understand strong local dependences from a geostatistician perspective, one can notice that the extremal coefficient and the -madogram of a max-stable process are related via


This implies that the spatial structure in any max-stable process can be related to the spatial structure of the input . More precisely, the madogram of the spectral process in (1) can be linked to the F-madogram by


This one-to-one link between the generative input and the output leads, in the absence of any nugget effect, to

Keeping in mind that , this limiting result implies that distance between the distributions at two nearby locations in the input process becomes twice smaller in the output process . So, creating strong extremal dependences implies the need of strong dependence in the generating process. On the contrary, a nugget effect, i.e. imperfect dependence even at infinitesimal distances, may appear in the output process only if it is present in the input process.

Still, Equality (6) tells us that the extremal coefficient will be ideally modeled if and only if is well captured. The later condition can only be satisfied if a very similar dependence in is built in, see relation (7). This reasoning leads to our main modeling idea. Instead of building complex parametric models for with inference schemes that typically result in high-dimensional optimization problems, see, for instance, Padoan, Ribatet and Sisson (2010); Dombry, Engelke and Oesting (2017); Huser et al. (2019), for likelihood-based inference methods or Oesting, Schlather and Friederichs (2017); Einmahl, Kiriliouk and Segers (2018) for (weighted) least square fits of certain summary statistics, we can “just” plug forecast ensemble members as max-stable constructions like the ones defined by (1), (3) or (4). Equalities (6) and (7) indicate that, if a subset of ensemble forecast members is well-chosen, then extremal coefficients measured from the weather stations should be well reproduced. By construction, such a model is based on a very small number of parameters only, and, thus, is easy to fit and simulate.

Before closing this section, we can note that extremal coefficients and madograms only provide specific information about pairwise dependences and do not capture multivariate features. Still, the same ideas could be extended to multivariate versions and complete dependence (see Marcon et al., 2017; Naveau et al., 2009).

3 Description of Rainfall Data

The mainland French territory witnesses complex spatial rainfall weather patterns due to its changing orography, the influence of the Atlantic Ocean, the Mediterranean Sea and different small and large climatological factors such as NOA. Another important aspect when modeling rainfall distributions is the quality of the data at hand. Here, our goal is to build our statistical analysis from the reference network of Météo France that is composed of 110 well kept weather stations, see the white and black dots in Figure 1 (top). These high quality stations recorded daily rainfall amounts over the time period 1980–2017. The elevation map in the bottom panel of Figure 1 displays a strong orography over the south of France, see e.g. the Cèvennes region north of Montpellier, the Pyrenees in the southwest and the Alps in the east. These geographical features suggest that either anisotropy, non-stationarity or both can be expected to be present in the spatial component of heavy rainfall, even after removing spatial trends at the marginal levels. Autumnal moisture brought from the Mediterranean sea can lead to severe convective storms with a specific spatial structure, different from northern weather patterns stopped by the Pyrenees.

As seen in Section 2, max-stable statistical models can capture strong dependence among maxima, but they are not appropriate to model weak dependencies for extremes (cf. Wadsworth and Tawn, 2012, for instance). In contrast to temperature extremes like heat waves, heavy rainfall are not likely to be dependent over very large regions. Precipitation extremes recorded at two stations more than 500 kilometers apart are likely to be independent. For this reason, we reduce our area of interest from the whole mainland territory to the southern part of France, see the black dots and the box in the top panel and the corresponding altitude map in the bottom panel of Figure 1. We have chosen this particular region because most severe heavy rainfall events occur there.

Figure 1: Top panel: our 39 locations of interest from the 110 reference Météo France weather stations network; the finer and homogeneous grid in grey corresponds to PEARP points of ensemble forecasts. Bottom: altitude map of our region of interest to highlight the strong orography from sea level to 4000 meters; the elevation data are retrieved from the Terrain Tiles on Amazon Web Services via the R package elevatr (Hollister and Shah, 2017).

3.1 Daily rainfall recorded by weather stations

We consider fall (SON) daily precipitation at stations in the south of France, see black dots in Figure 1, over 38 years, from 1980 to 2017. Each fall season has been divided into five blocks of length 18. Maxima were computed over each of these blocks. The top panel of Figure 2

displays the estimated GEV shape parameters obtained by a probability-weighted moment method

(Diebolt et al., 2008) separately for each station, assuming independence among blocks. Kolmogorov-Smirnov tests per stations indicate a good fit of the estimated GEV distributions with an average -value of . The value range and the spatial pattern of estimated is roughly similar to the ones observed in previous studies (see, e.g. Carreau, Naveau and Neppel, 2017; Carreau et al., 2012; Blanchet, Molinié and Touati, 2018).

Figure 2: Top: Map of measurement stations and estimated GEV shape parameters; Bottom: Pairwise empirical extremal coefficients, see Equation (6), plotted against the distance between the stations.

With the analysis of the marginal distributions ensuring the compatibility of the station data with a max-stable model, we will henceforth focus on our main objective, i.e. modeling the spatial structure among heavy rainfall. To this end, we first investigate the extremal dependence structure between rainfall data at different stations as described in Section 2. More precisely, for each pair of stations, we estimate the pairwise extremal coefficient via a rank-based empirical version of the weighted -madogram, see Marcon et al. (2017). In the bottom panel of Figure 2, the dependence captured by the estimated pairwise extremal coefficient decreases as the distance between stations increases. A nugget effect seems to be present because extremal coefficients for very small distances do not appear to be close to one, but rather around 1.2.

The main task is to produce max-stable processes that can reproduce such spatial features at finer scales. To this end, we will make use of the second type of data, namely precipitation forecasts.

3.2 Gridded daily precipitation produced from weather forecast center

The national French weather service, Météo France, produces, on a daily basis, an ensemble of 35 members with forecasted daily precipitation at 17596 cells of size covering the mainland of France, see the grid displayed in Figure 1. To merge these forecasted data with our observational network described in Section 3.1, we extract a subset of grid cells over a rectangular region that contains our stations (see the box in the bottom panel of Figure 1) and over a time period that overlaps, more precisely Fall seasons from 2012 to 2017, i.e. 546 days. Let denote the forecast of the th ensemble member for grid cell and day . In order avoid potential scaling problems, for each ensemble member and each grid cell , all the forecasts are transformed to a unit Fréchet scale via a rank transform

where denotes the number of days for which the forecast of ensemble member is available at grid cell . This results in transformed forecast maps. As we will treat all the ensemble members in the same way, henceforth, for simplicity, we will denote these maps , with , of transformed forecasted daily precipitation. These forecasts are used to construct data driven max-stable models for extreme precipitation.

4 Statistical models

In this section, five different max-stable models are studied and compared. For simplicity and analogously to Section 2, all models are defined in a standardized way with unit Fréchet margins. While one of the five max-stable models is fully parametric and based on a Brown–Resnick process with a nugget effect, the other four will be data driven by the maps from Section 3.2. To assess all five models, we will focus on their extremal dependence structure. In particular, the extremal coefficients associated to each model will be compared to their empirically estimated counterparts in Figure 3. In this comparison, each station is spatially identified to its closest grid cell. The resulting root mean squared error (RMSE) will be used to evaluate the quality of the model fit. Estimation uncertainty will be assessed by a parametric bootstrap. More precisely, 190 block maxima are simulated from each of the fitted models 500 times. For each of these 500 simulations, the pairwise extremal coefficients are estimated. The intervals between empirical - and

-quantiles of these samples are displayed in gray, indicating the region in which the empirical estimates would be expected to be if the fitted model was correct.

The five models are listed below.

  • Model A: As a starting point, we assume that all maps contain relevant information about rainfall maxima. Hence, the standard max-linear model (3) is implemented by defining

    as normalized basis functions based on all the forecast vectors


According to the top left panel of Figure 3, Model A does not capture accurately the extremal dependence structure observed in biweekly maxima recorded at weather stations. This shortcoming can be explained by the incorrect assumption that all grid points of all daily fields are linked to extreme rainfall. Days with little or no rain, however, should be not used to build the basis functions. To account for this fact, we exploit the theory of generalized Pareto processes (Ferreira and De Haan, 2014). According to the theory, extremal dependence in the forecasts should fully described by the max-spectral functions for those such that for some high threshold . Thus, we will use these spectral functions to build new basis functions. Assuming, without loss of generality, that the vectors are sorted w.r.t. their maximum, i.e. , this results in a number of forecasts to be taken into account. For simplicity, for each of the following models, Model B, C and D, we will choose a fixed threshold as the empirical -quantile of the vector . This choice leads to maps used for the construction of the models.

  • Model B: As an improvement of Model A, we consider a max-linear model with max-spectral functions built from those forecast exceeding the threshold . The resulting normalized basis functions are given by

The top right panel in Figure 3 indicates that the extremal coefficients for pairs of strongly dependent stations are still systematically smaller than their empirical counterparts, i.e. dependence in the model is stronger than dependence in the data – a phenomenon that is often present when comparing forecasts to observations and that can be explained by the presence of some nugget effect in the observed data. Models C and D provide two ways of incorporating nugget effects into max-stable models.

  • Model C: This model is based on Reich and Shaby (2012) construction defined by (4). The spectral functions are identical to the ones used in Model B.

  • Model D: To handle a possible nugget effect, we combine Model B with a noise process throughout a max-linear operator

    where is a mixture parameter, is a noise process with unit Fréchet marginal distributions and, independently of , is the max-linear process in Model B.

The additional parameters, in Model C and in Model D, respectively, are chosen such that the RMSE is minimized. By a first visual inspection, see Figure 3, models C and D appear to capture the observed extremal coefficients well.

The quality of these plots has to be interpreted in regard to the number of parameters inferred. This number is zero, one and one for model B, C and D, respectively. This highlights that our approach is very parsimonious in terms of parameter number and inference complexity. The threshold was set to the -quantile of for all models, and consequently, the RMSE could be even lower if the threshold choice was optimized for each of the three models separately. But, as our goal is to propose a straightforward estimation scheme, we refrain from optimizing the choice of individual thresholds for models B, C and D.

As previously mentioned, our last model is different and based on a classical and fully parametric max-stable model, the Brown–Resnick process.

  • Model E: We use a Brown–Resnick process associated to the variogram

    for some , , , . Here, the five parameters are chosen such that the RMSE of the estimated pairwise extremal coefficients is minimized.

By construction, this model cannot capture non-stationarity and increasing the number of parameters to do so will be non-trivial. Still, Model E offers some flexibility in terms of anisotropy and nugget effects.

Figure 3: For each of the fitted models, Model A–E, the theoretical pairwise extremal coefficients are plotted against the empirical estimates (black); the gray intervals indicate the uncertainty of the estimates based on simulations from the fitted models.

In the following, we will compare the classical Model E to the two data driven models C and D which provide the best fit in terms of pairwise extremal coefficients as indicated by the root mean squared errors displayed in Figure 3. It can be seen that the theoretical extremal coefficients according to models C and D are almost identical apart from the coefficients for the strongly dependent pairs where dependence seems to be slightly stronger in Model D than in Model C. The coefficients according to Model E are also largely similar to the ones obtained by the data driven models. The main difference is that the fitted model E exhibits slightly larger extremal coefficients than the empirical ones for strongly dependent pairs, while the theoretical coefficients are slightly smaller than their empirical counterparts for weakly dependent pairs of stations. These inaccuracies are also reflected by the root mean squared error for model E. The relative gain from model E to model D is 16% in terms of RMSE. Analogously to the pairwise dependence structures, one may also compare the fits for summary statistics of higher order. The corresponding results for the triplewise extremal coefficients can be found in the appendix.

While the summary statistics considered so far provide some information about extremal dependence along the diagonal of bivariate and trivariate distributions, respectively, some further insight in the three models may be gained by regarding artificial precipitation fields obtained by simulations. Such realizations from the three models are displayed in Figure 4. Although such graphical comparisons are only qualitative, models C and D appear to be able provide a wide range of spatial features with specific regional behavior, having a slightly stronger nugget effect than Model E. They can reproduce very localized events, but also generate simulated fields that appear to have well defined spatial structures along geographical features, see the last two rows.

The stationary model E with five parameters appear to simulate spatial structures similar to the simpler models C and D with one parameter only. This fact leads to less expensive inference schemes for the data driven models. Furthermore, due to their simple structure as given in Equations (3) and (4), respectively, the data driven models are easy to simulate, while simulation of Brown–Resnick or other parametric spatial max-stable models on dense grids is typically computationally intensive (cf. Dombry, Engelke and Oesting, 2016; Oesting, Schlather and Zhou, 2018, for instance).

Figure 4: Simulations from Models C, D and E.

5 Discussion and conclusion

Although climatologists, weather forecasters and statisticians have been collaborating extensively the last decades, the field of data assimilation being a successful example of such a joint research effort, the extreme value theory community has been slower at integrating new data sources within their multivariate extremal models. For example, there are many high quality methodological articles on heavy rainfall analysis, but most are based on complex parametric models applied to one unique data source. This leads to non-trivial inference problems and such approaches can be difficult to transfer to researchers outside of this particular domain.

In this work, our goal was to show that the framework of max-stable processes, one pillar of spatial extreme value analysis, can be easily coupled with other data sources. More specifically, simple max-stable processes that integrate ensemble forecast rainfall data as spectral profiles were able to reproduce the main spatial features of heavy rainfall over a complex climatological region. We also show that our approach compares favorably with a more complicated parametrized model. In addition, it is simple to handle non-stationarity, and both inference and simulation are straightforward and quick. One obvious limitation of our method is that, as expected from the theory of max-stable processes, the efficiency of our approach directly depends on the quality of the input data. If weather services were unable to reproduce adequately important spatial features of storms and fronts in their forecast ensembles, our strategy will naturally be inefficient.

To conclude, the production of numerical models outputs, their capabilities, their associated resolution and their sizes appear to have increased rapidly these last few years, and this trend is likely to continue. In the context of extreme value analysis, it is always a delicate question to know if such numerical models can simulate adequately extreme events, or produce even unobserved ones. From a geophysical point of view, most of these numerical models are physically consistent, and consequently should contain some meaningful information about spatial and temporal structures. Besides the case study presented in this paper, it would be interesting to determine if other extremal models, e.g. asymptotic independence, could benefit of such additional information to improve the estimation of very high quantiles, and if so, how to combine them to produce spatio-temporal extremal fields in compliance with EVT in a physically coherent way.


We would like to thank Météo France for providing the data used in this work. Part of P. Naveau’s work was supported by the European DAMOCLES-COST-ACTION on compound events, and also benefited from French national programs, in particular FRAISE-LEFE/INSU, MELODY-ANR, and ANR-11-IDEX-0004 -17-EURE-0006.


  • Bechler, Bel and Vrac (2015) [author] Bechler, AurélienA., Bel, LilianeL. Vrac, MathieuM. (2015). Conditional simulations of the extremal process: application to fields of extreme precipitation. Spat. Stat. 12 109–127.
  • Blanchet, Molinié and Touati (2018) [author] Blanchet, JulietteJ., Molinié, GillesG. Touati, JulienJ. (2018). Spatial analysis of trend in extreme daily rainfall in southern France. Clim. Dyn. 51 799–812.
  • Brown and Resnick (1977) [author] Brown, Bruce M.B. M. Resnick, Sidney I.S. I. (1977). Extreme Values of Independent Stochastic Processes. J. Appl. Probab. 14 732–739.
  • Buhl and Klüppelberg (2016) [author] Buhl, SvenS. Klüppelberg, ClaudiaC. (2016). Anisotropic Brown–Resnick space-time processes: estimation and model assessment. Extremes 19 627–660.
  • Carreau, Naveau and Neppel (2017) [author] Carreau, JulieJ., Naveau, PP. Neppel, LL. (2017). Partitioning into hazard subregions for regional peaks-over-threshold modeling of heavy precipitation. Water Resour. Res. 53 4407–4426.
  • Carreau et al. (2012) [author] Carreau, JulieJ., Neppel, LucL., Arnaud, PatrickP. Cantet, PhilippeP. (2012). Extreme rainfall analysis at ungauged sites in the South of France: comparison of three approaches. J. Soc. Fr. Statistique.
  • Chen et al. (2018) [author] Chen, YangY., Moufouma-Okia, WilfranW., Masson-Delmotte, ValerieV., Zhai, PanmaoP. Pirani, AnnaA. (2018). Recent Progress and Emerging Topics on Weather and Climate Extremes Since the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Annu. Rev. Environ. Resour. 43 35–59.
  • Cooley, Naveau and Poncet (2006) [author] Cooley, DanielD., Naveau, PhilippeP. Poncet, PaulP. (2006). Variograms for spatial max-stable random fields. In Dependence in Probability and Statistics 373–390. Springer, New York.
  • Cooley, Nychka and Naveau (2007) [author] Cooley, DanielD., Nychka, DouglasD. Naveau, PhilippeP. (2007). Bayesian spatial modeling of extreme precipitation return levels. J. Am. Stat. Assoc. 102 824–840.
  • De Fondeville and Davison (2018) [author] De Fondeville, RaphaëlR. Davison, Anthony CA. C. (2018). High-dimensional peaks-over-threshold inference. Biometrika 105 575–592.
  • De Haan (1984) [author] De Haan, LaurensL. (1984). A spectral representation for max-stable processes. Ann. Probab. 12 1194–1204.
  • De Haan and Ferreira (2006) [author] De Haan, LaurensL. Ferreira, AnaA. (2006). Extreme Value Theory: An Introduction. Springer, New York.
  • Diebolt et al. (2008) [author] Diebolt, JeanJ., Guillou, ArmelleA., Naveau, PhilippeP. Ribereau, PierreP. (2008). Improving probability-weighted moment methods for the generalized extreme value distribution. Revstat Stat. J. 6 33–50.
  • Dombry, Engelke and Oesting (2016) [author] Dombry, C.C., Engelke, S.S. Oesting, M.M. (2016). Exact Simulation of Max-Stable Processes. Biometrika 103 303–317.
  • Dombry, Engelke and Oesting (2017)

    [author] Dombry, ClémentC., Engelke, SebastianS. Oesting, MarcoM. (2017). Bayesian inference for multivariate extreme value distributions. Electron. J. Stat. 11 4813–4844.

  • Einmahl, Kiriliouk and Segers (2018) [author] Einmahl, John HJJ. H., Kiriliouk, AnnaA. Segers, JohanJ. (2018). A continuous updating weighted least squares estimator of tail dependence in high dimensions. Extremes 21 205–233.
  • Embrechts, Klüppelberg and Mikosch (1997) [author] Embrechts, PaulP., Klüppelberg, ClaudiaC. Mikosch, ThomasT. (1997). Modelling Extremal Events for Insurance and Finance. Springer-Verlag, Berlin.
  • Ferreira and De Haan (2014) [author] Ferreira, AnaA. De Haan, LaurensL. (2014). The generalized Pareto process; with a view towards application and simulation. Bernoulli 20 1717–1737.
  • Fisher and Tippett (1928) [author] Fisher, Ronald AylmerR. A. Tippett, Leonard Henry CalebL. H. C. (1928). Limiting forms of the frequency distribution of the largest or smallest member of a sample. Math. Proc. Camb. Philos. Soc. 24 180–190.
  • Hollister and Shah (2017) [author] Hollister, JeffreyJ. Shah, TarakT. (2017). elevatr: Access Elevation Data from Various APIs R package version 0.1.3.
  • Huser et al. (2019) [author] Huser, RaphaëlR., Dombry, ClémentC., Ribatet, MathieuM. Genton, Marc GM. G. (2019). Full likelihood inference for max-stable data. Stat 8 e218.
  • Kabluchko, Schlather and de Haan (2009) [author] Kabluchko, ZakharZ., Schlather, MartinM. de Haan, LaurensL. (2009). Stationary Max-Stable Fields Associated to Negative Definite Functions. Ann. Probab. 37 2042–2065.
  • Keef, Svensson and Tawn (2009) [author] Keef, CarolineC., Svensson, CeciliaC. Tawn, Jonathan AJ. A. (2009). Spatial dependence in extreme river flows and precipitation for Great Britain. J. Hydrol. 378 240–252.
  • Leadbetter, Lindgren and Rootzén (1983) [author] Leadbetter, M. RossM. R., Lindgren, GeorgG. Rootzén, HolgerH. (1983). Extremes and Related Properties of Random Sequences and Processes. Springer-Verlag, New York.
  • Marcon et al. (2017) [author] Marcon, GiuliaG., Padoan, Simone AS. A., Naveau, PhilippeP., Muliere, PietroP. Segers, JohanJ. (2017). Multivariate nonparametric estimation of the Pickands dependence function using Bernstein polynomials. J. Stat. Plan. Inference 183 1–17.
  • Naveau et al. (2009) [author] Naveau, PhilippeP., Guillou, ArmelleA., Cooley, DanielD. Diebolt, JeanJ. (2009). Modelling pairwise dependence of maxima in space. Biometrika 96 1–17.
  • Oesting (2018) [author] Oesting, MarcoM. (2018). Equivalent representations of max-stable processes via -norms. J. Appl. Probab. 55 54–68.
  • Oesting, Schlather and Friederichs (2017) [author] Oesting, MarcoM., Schlather, MartinM. Friederichs, PetraP. (2017). Statistical post-processing of forecasts for extremes using bivariate Brown-Resnick processes with an application to wind gusts. Extremes 20 309–332.
  • Oesting, Schlather and Zhou (2018) [author] Oesting, MarcoM., Schlather, MartinM. Zhou, ChenC. (2018). Exact and Fast Simulation of Max-Stable Processes on a Compact Set Using the Normalized Spectral Representation. Bernoulli 24 1497–1530.
  • Padoan, Ribatet and Sisson (2010) [author] Padoan, Simone AS. A., Ribatet, MathieuM. Sisson, Scott AS. A. (2010). Likelihood-Based Inference for Max-Stable Processes. J. Am. Stat. Assoc. 105 263–277.
  • Reich and Shaby (2012) [author] Reich, Brian JB. J. Shaby, Benjamin AB. A. (2012). A hierarchical max-stable spatial model for extreme precipitation. Ann. Appl. Stat. 6 1430–1451.
  • Resnick (2008) [author] Resnick, Sidney IS. I. (2008). Extreme Values, Regular Variation and Point Processes. Springer, New York.
  • Saunders et al. (2017) [author] Saunders, KateK., Stephenson, Alec GA. G., Taylor, Peter GP. G. Karoly, DavidD. (2017). The spatial distribution of rainfall extremes and the influence of El Niño Southern Oscillation. Weather Clim. Extrem. 18 17–28.
  • Sillmann et al. (2017) [author] Sillmann, JanaJ., Thorarinsdottir, ThordisT., Keenlyside, NoelN., Schaller, NathalieN., Alexander, Lisa VL. V., Hegerl, GabrieleG., Seneviratne, Sonia IS. I., Vautard, RobertR., Zhang, XuebinX. Zwiers, Francis WF. W. (2017). Understanding, modeling and predicting weather and climate extremes: Challenges and opportunities. Weather Clim. Extrem. 18 65–74.
  • Taillardat et al. (2016) [author] Taillardat, MaximeM., Mestre, OlivierO., Zamo, MichaëlM. Naveau, PhilippeP. (2016). Calibrated ensemble forecasts using quantile regression forests and ensemble model output statistics. Mon. Weather Rev. 144 2375–2393.
  • Taillardat et al. (2019) [author] Taillardat, MaximeM., Fougères, Anne-LaureA.-L., Naveau, PhilippeP. Mestre, OlivierO. (2019). Forest-based and semi-parametric methods for the postprocessing of rainfall ensemble forecasting. Weather and Forecasting 617-633.
  • Wadsworth and Tawn (2012) [author] Wadsworth, Jennifer LJ. L. Tawn, Jonathan AJ. A. (2012). Dependence modelling for spatial extremes. Biometrika 99 253–272.
  • Wang and Stoev (2011) [author] Wang, YizaoY. Stoev, Stilian A.S. A. (2011). Conditional sampling for spectrally discrete max-stable random fields. Adv. Appl. Probab. 43 461–483.

Trivariate dependence structure

While pairwise extremal coefficients are used for model fitting, we can also consider extremal coefficients of higher order. Analogously to Equation (5), for a max-stable process with marginal distribution and , the extremal coefficient can be defined via

By construction, – a quantity that is often interpreted as the number of independent random variables among . Likewise the pairwise coefficients, the higher order coefficients can be estimated via the empirical multivariate weighted -madogram, see Marcon et al., 2017).

In Figure 5, results for the estimated triplewise extremal coefficients are shown and compared to the theoretical ones for the fitted models. Similarly to the results for the pairwise coefficients displayed in Figure 3, the root mean square error for Model E is slightly worse than the errors for Model C and Model D, respectively. Furthermore, it can be seen that the triplewise extremal coefficients of Models C and D are nearly identical, while the coefficients of Model E show stronger deviations.

Figure 5: Triplewise extremal coefficients for three fitted models, Model C, Model D and Model E. Left: analogously to Figure 3. Right: Comparison of the three models by plotting the theoretical triplewise extremal coefficients against one another.