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 3.0 of the code has the same features as version 2.1, as described below, but makes three changes (see below).

Changes:

• A new parameter that captures the spatially averaged value of the GEV shape parameter has been introduced. This enables BAYEX to model a wider range of shape values.
• A new parameterization has been used to enforce the support of the GEV distribution. This leads to more robust estimates as well as predictions.
• This version of the code runs on CmdStanPy, instead of PyStan.

The code was originally developed and described in the paper:
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).

This version is being released to support the following paper:
Morim, J., Wahl, T., Rasmussen, DJ, Calafat, F. M., Vitousek, S., Dangendorf, S., Kopp, R. E., Oppenheimer, M. (2024). Observations reveal changing coastal storms along US coastlines. Submitted to Nature Climate Change.

General features:

• Users can 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 model temporal variability in the GEV location parameter using either an integrated random walk or a linear trend.
• Users can choose to fit the model using either full Gaussian processes (GPs) or Nearest Neighbor Gaussian Process (NNGP) for the distribution parameters. 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 interation in NNGPs is linear in the number of data sites and NNGPs do not require storing or inverting large matrices.
• Users can choose to build the residual spatial process at each site using either all spatial knots or a smaller subset of knots (a form of Nearest Neighbor approximation). The latter reduces the computational complexity of the model.