margEVT is an R package for non-stationary extreme value
analysis under covariate-driven regimes. It couples a covariate-driven
Non-Homogeneous Poisson Process (NHPP) with (adaptive) elastic-net
penalized maximum likelihood and stochastically marginalized
return-level estimation, with conditional bootstrap intervals and
walk-forward validation.
Exceedances of an extreme threshold \(u\) over an observational domain are modeled as a Non-Homogeneous Poisson Process on the product space \(\mathcal{S} = [0, 1] \times (u, \infty)\). The time-varying intensity function \(\lambda(t, y)\) is parameterized using the location \(\mu(t)\), scale \(\sigma(t) > 0\), and shape \(\xi(t)\) parameters of the Generalized Extreme Value (GEV) distribution:
\[\lambda(t, y) = \frac{1}{\sigma(t)} \left[ 1 + \xi(t) \left( \frac{y - \mu(t)}{\sigma(t)} \right) \right]_{+}^{-1/\xi(t)-1}\]
where \([a]_+ = \max(a, 0)\).
To capture non-stationarity, the parameters are modeled as linear combinations of covariate vectors \(\mathbf{x}_t\):
\[\mu(t) = \mu_0 + \sum_{k=1}^{K} \mu_k x_{k,t}\]
\[\log(\sigma(t)) = \phi_0 + \sum_{j=1}^{J} \phi_j x_{j,t}\]
\[\xi(t) = \xi_0 + \sum_{l=1}^{L} \xi_l x_{l,t}\]
The logarithmic parameterization of \(\sigma(t) = \exp(\phi(t))\) guarantees
positivity. The shape is restricted to \(\xi(t) \ge -1/2\) by default
(shape_lower), the region where the maximum likelihood
estimator is regular (Smith, 1985); fits at the bound are handled by an
active-set solver.
In high-dimensional covariate settings, parameter estimation is conducted via penalized maximum likelihood. The objective function minimizes the negative log-likelihood \(\ell(\boldsymbol{\theta})\) subject to a weighted elastic-net penalty:
\[Q(\boldsymbol{\theta}; \lambda, \alpha) = -\ell(\boldsymbol{\theta}) + \lambda \sum_{k \in \mathcal{P}} w_k \left[ \alpha \sqrt{\theta_k^2 + \varepsilon} + \left( \frac{1 - \alpha}{2} \right) \theta_k^2 \right]\]
where: - \(\lambda \ge 0\) controls
the overall regularization strength. - \(\alpha \in [0, 1]\) controls the balance
between the \(\ell_1\) (LASSO) and
\(\ell_2\) (Ridge) components (setting
\(\alpha = 1\) recovers the pure LASSO
penalty). - \(w_k\) are coordinate
weights: \(w_k = 1\) by default, or
\(w_k = 1/|\tilde\theta_k|^\gamma\)
from an unpenalized fit with penalty_weights = "adaptive"
(the adaptive LASSO). - \(\varepsilon >
0\) is a small smoothing parameter ensuring differentiability at
the origin. - \(\mathcal{P}\) is the
index set of penalized parameters. The baseline intercepts (\(\mu_0, \phi_0, \xi_0\)) and deterministic
seasonal Fourier harmonics are kept unpenalized by construction (\(\mathcal{P} \cap \text{unpenalized} =
\emptyset\)).
The optimization is solved via the quasi-Newton BFGS algorithm (with L-BFGS-B as a fallback) supplied with exact analytical gradients, and convergence is certified by the Karush-Kuhn-Tucker conditions of the penalized problem:
\[\nabla_{\boldsymbol{\theta}} Q(\boldsymbol{\theta}; \lambda, \alpha) = -\nabla_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}) + \nabla_{\boldsymbol{\theta}} P_{\lambda,\alpha}(\boldsymbol{\theta})\]
The regularization path parameter \(\lambda^{\ast}\) is selected by minimizing
the Bayesian Information Criterion (BIC). Automatic selection can either
use a common scalar penalty for all blocks
(lambda_scaling = "common", the default) or score-magnitude
calibration across the location, log-scale, and shape blocks
(lambda_scaling = "gradient"):
\[\text{BIC}(\lambda) = -2\ell(\hat{\boldsymbol{\theta}}_\lambda) + k_\lambda \log(m)\]
where \(m\) is the exceedance count
supplied to the fitted point-process model and \(k_\lambda\) is the number of active
parameters under the model’s operational tolerance. The factor \(\log(m)\) is an event-based model-selection
convention, not a universal effective-sample-size identity.
Declustering, when required by the application, is performed before
calling fit_nhpp().
Under non-stationarity, traditional definitions of a \(T\)-year return level are conceptually ill-defined. This package implements three distinct frameworks:
The user-specified stochastic covariates are fixed at a scenario \(\mathbf{x}^{\ast}\), representing a hypothetical frozen climate state while any deterministic predictors retained in the fitted design are evaluated over the annual block. If every predictor is constant, the conditional return level reduces to the analytical expression:
\[z_T(\mathbf{x}^{\ast}) = \mu(\mathbf{x}^{\ast}) + \frac{\sigma(\mathbf{x}^{\ast})}{\xi(\mathbf{x}^{\ast})} \left[ \left(-\log\left(1 - \frac{1}{T}\right)\right)^{-\xi(\mathbf{x}^{\ast})} - 1 \right]\]
The package evaluates the general conditional annual path numerically, so this approach also accommodates deterministic seasonal terms.
To capture long-run risk over the natural variability of the climate system, the non-stationary intensity is integrated over the stationary joint distribution \(\Pi\) of the covariate trajectories \(\mathbf{v}\):
\[G_{\Pi}(z) = \mathbb{E}_{\{\mathbf{v} \sim \Pi\}} \left\lbrack G(z \mid \mathbf{v}) \right\rbrack = \mathbb{E}_{\{\mathbf{v} \sim \Pi\}} \left\lbrack \exp \left\lbrace -\frac{1}{n_{y}} \sum_{j=1}^{n_{y}} \left\lbrack 1 + \xi(t_{j}) \left\lparen \frac{z - \mu(t_{j} \mid \mathbf{v})}{\sigma(t_{j} \mid \mathbf{v})} \right\rparen \right\rbrack_{{+}}^{-1/\xi(t_{j})} \right\rbrace \right\rbrack\]
The joint distribution \(\Pi\) is modeled via a stationary Vector Autoregressive process, \(\text{VAR}(p)\), fitted by least squares with Schwarz order selection (a single covariate is a \(\text{VAR}(p)\) with \(K = 1\)). The fitted generator verifies that the companion-root spectral radius is below one before simulation. Independent synthetic annual trajectories are initialized from the fitted stationary Gaussian state distribution at the fitted temporal resolution, Fourier seasonality is re-injected, and \(G_{\Pi}(z)\) is estimated via Monte Carlo integration over \(n_{mc}\) simulated years. A zero-state burn-in initialization remains available as a sensitivity option:
\[\hat{G}_{\Pi, B}(z) = \frac{1}{n_{mc}} \sum_{r=1}^{n_{mc}} G(z \mid \mathbf{v}^{(r)})\]
The marginalized plug-in return level is recovered numerically from a target crossing satisfying:
\[\hat{G}_{\Pi, B}(z_T^{\Pi}) - \left( 1 - \frac{1}{T} \right) = 0\]
The trajectory law \(\Pi\) is replaced by the empirical distribution of eligible historical annual blocks. If \(K\) complete blocks are available, the exact empirical plug-in distribution is
\[G_{\hat\Pi_{\mathrm{emp}}}(z) = \frac{1}{K} \sum_{k=1}^{K} G(z \mid \mathbf{V}_k).\]
The implementation approximates this finite empirical average by
sampling n_boot complete years with replacement. This
preserves the observed within-year covariate paths and is used as a
sensitivity benchmark for the parametric trajectory generator.
bootstrap_rl(), bootstrap_coef() and
bootstrap_nhpp() give parametric bootstrap intervals for
return levels and coefficients, conditional on the selected penalty
(with occurrence = "simulated", the number and timing of
exceedances are also redrawn). backtest() refits the model
on expanding training windows and compares the observed annual maxima of
the following years with the nominal exceedance rates of approaches B
and C. residuals() and plot(fit, type = "qq")
provide time-rescaling diagnostics.
You can install the stable release version of margEVT
from CRAN:
install.packages("margEVT")Alternatively, you can install the development version from GitHub:
# install.packages("devtools")
devtools::install_github("rodrigosqrt3/margEVT")This reproducible example simulates a generic non-stationary extreme value process, fits the regularized NHPP model using the LASSO penalty (\(\alpha = 1\)) with BIC-based \(\lambda\) selection, and computes return levels under all three major marginalization frameworks (including the stochastically marginalized Approach B).
library(margEVT)
# =============================================================================
# 1. Simulate Non-Stationary Extreme Value Data
# =============================================================================
set.seed(3)
n_years <- 10L
n_per_year <- 100L
n <- n_years * n_per_year
years <- rep(2011:(2011 + n_years - 1L), each = n_per_year)
# Simulate 4 generic standardized covariates
x1 <- rnorm(n, mean = 0, sd = 1)
x2 <- rnorm(n, mean = 0, sd = 1)
x3 <- rnorm(n, mean = 0, sd = 1)
x4 <- rnorm(n, mean = 0, sd = 1)
# Define non-stationary GEV parameters
# Location varies with x1 and x3; Scale (log) varies with x2 and x4
mu <- 15 + 2.0 * x1 + 1.5 * x3
sigma <- exp(1.2 + 0.3 * x2 + 0.1 * x4)
xi <- -0.10 # Bounded upper tail (Weibull domain)
# Generate non-stationary GEV observations via the inverse transform method
u_rand <- runif(n)
y <- mu + sigma * (((-log(u_rand))^(-xi) - 1) / xi)
sim_data <- data.frame(
y = y,
x1 = x1,
x2 = x2,
x3 = x3,
x4 = x4,
year = years
)
# =============================================================================
# 2. Fit the Penalized Non-Stationary NHPP Model
# =============================================================================
fit <- fit_nhpp(
df = sim_data,
threshold = 19,
loc_vars = c("x1", "x2", "x3"),
scale_vars = c("x1", "x2", "x4"),
shape_vars = NULL, # Stationary shape parameter
penalty = "lasso",
alpha = 1.0,
lambda = "bic",
obs_per_year = 100, # 100 observations per annual block
calc_hessian = TRUE,
verbose = TRUE
)
# Inspect estimated parameter coefficients and Hessian conditioning
summary(fit)
# Extract raw coefficients
coef(fit)
# =============================================================================
# 4. Fit the Parametric VAR(p) Stochastic Generator
# =============================================================================
# Automatically deseasonalizes active covariates, selects the optimal lag p
# via BIC, and fits the multivariate VAR model to the anomalies.
generator <- fit_var_generator(fit, sim_data)
# Simulate 50 years of stochastic trajectories from the fitted generator,
# automatically re-injecting the seasonal harmonics.
mc_sample <- simulate_covariates(generator, n_mc = 50L, n_obs = 100L, seed = 123L)
# =============================================================================
# 5. Compute Unconditional Return Levels (Core Contribution)
# =============================================================================
# Rather than conditioning on a "frozen" covariate state (Approach A),
# we integrate out covariate uncertainty under both the parametric generator
# (Approach B) and empirical annual-block marginalization (Approach C).
rl_long <- marginalize(
fit = fit,
data = sim_data,
TRs = c(2, 5, 10, 50), # Return periods (years)
n_obs = 100L, # 100 observations per year
approaches = c("A", "B", "C"),# Evaluate all three approaches
scenarios = list("baseline" = list(x1 = 0, x2 = 0, x3 = 0, x4 = 0)),
mc_sample = mc_sample, # Supplied simulated sample from our VAR(p)
n_boot = 100L, # 100 bootstrap years for Approach C
seed = 42L,
year_col = "year"
)
# Convert the long-format return levels into a wide-format risk table
rl_wide <- rl_table(rl_long)
print(rl_wide)To cite margEVT in publications, please use:
Villa, R. F. (2026). Regularized Poisson Point Process Models and Stochastic Marginalization for Return-Level Inference under Covariate-Driven Extremes (Master’s dissertation, Instituto de Matemática e Estatística, Universidade Federal do Rio Grande do Sul, Porto Alegre, Brazil. Advisor: Flavio Ziegelmann).
Alternatively, run citation("margEVT") in R once the
package is installed.