BAYEX: Spatiotemporal Bayesian hierarchical modeling of extremes with max-stable processes

BAYEX offers spatiotemporal Bayesian hierarchical modeling of extremes using max-stable and latent processes. As a key feature, BAYEX makes estimates of both the GEV parameters (location, scale and shape) and the annual maxima at any arbitrary location, either gauged or ungauged, while providing realistic uncertainty estimates. Inference in BAYEX is performed using Hamiltonian Monte Carlo as implemented by the Stan probabilistic programming language.

This version 2.0 of the code adds new features to the previous release (version 1.0), as described below. Please cite the following paper when using this code:

Calafat, F. M., and M. Marcos (2020), Probabilistic reanalysis of storm surge extremes in Europe. Proc. Natl. Acad. Sci. U. S. A. 117 (4), 1877-1883.

The Bayesian hierarchical model implemented by this code is based on the approach developed in the paper below, but with several modifications to how the spatiotemporal evolution of the GEV parameters is modelled:
Reich, B. J., and B. A. Shaby, A hierarchical max-stable spatial model for extreme precipitation. Ann. Appl. Stat. 6, 1430–1451 (2012).

New features:

• Users can now choose to estimate the GEV shape parameter from the data or to set its value to zero effectively assuming a Gumbel distribution
• Allow the GEV shape parameter to vary in space.
• Users can now model temporal variability in the GEV location parameter using either an integrated random walk or a linear trend.
• Add Nearest Neighbor Gaussian Process (NNGP) for the GEV parameters. Users can now choose to fit the model using either full Gaussian processes (GPs) or NNGPs. Full GPs require ~O(n^3) flops and ~O(n^2) storage (n: number of data sites), and so they can rapidly become computationally infeasible. The total flop count per iteration in NNGPs is linear in the number of data sites and NNGPs do not require storing or inverting large matrices. This means that the users are now able to fit the model to larger datasets. The implementation of the NNGP is based on the formulation given by Datta et al. (2016) (Datta, A., S. Banerjee, A. O. Finley, and A. E. Gelfand 2016. Hierarchical Nearest-Neighbor Gaussian Process Models for Large Geostatistical Datasets. J. Am. Stat. Assoc. 111, 800–812).
• Add Nearest Neighbor approximation for the residual spatial process. Users can now also choose to build the residual spatial process at each site using either all spatial knots or a smaller subset of knots. The latter further reduces the computational complexity of the model.