README

margEVT is an R package developed to conduct non-stationary extreme value analysis under covariate-driven regimes. The package implements the statistical framework developed in Villa (2026) under the supervision of Prof. Dr. Flavio Ziegelmann, coupling a covariate-driven Non-Homogeneous Poisson Process (NHPP) with an Elastic-Net penalized maximum likelihood estimation framework and stochastically marginalized return-level estimation.

1. Mathematical Framework

1.1 Non-Homogeneous Poisson Process (NHPP) Representation

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}\]

To capture non-stationarity, the parameters are modeled as linear combinations of covariate vectors \(\mathbf{x}_t\):

The logarithmic parameterization of \(\sigma(t) = \exp(\phi(t))\) guarantees positivity.

1.2 Elastic-Net Regularization and Optimization

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 an Elastic-Net penalty:

\[Q(\boldsymbol{\theta}; \lambda, \alpha) = -\ell(\boldsymbol{\theta}) + \lambda \sum_{k \in \mathcal{P}} \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). - \(\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 or L-BFGS-B algorithm supplied with exact analytical gradients:

\[\nabla_{\boldsymbol{\theta}} Q(\boldsymbol{\theta}; \lambda, \alpha) = -\nabla_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}) + \nabla_{\boldsymbol{\theta}} P_{\lambda,\alpha}(\boldsymbol{\theta})\]

The optimal regularization path parameter \(\lambda^{\ast}\) is selected by minimizing the Bayesian Information Criterion (BIC):

\[\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().

2. Return Level Inference

Under non-stationarity, traditional definitions of a \(T\)-year return level are conceptually ill-defined. This package implements three distinct frameworks:

2.1 Approach A: Static Conditional Return Level

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.

2.2 Approach B: Unconditional Parametric Stochastic Marginalization

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)\). The fitted generator verifies that the companion-root spectral radius is below one before simulation. Independent synthetic annual trajectories are generated with separate burn-in periods 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:

\[\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:

2.3 Approach C: Empirical Annual-Block Marginalization

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.

3. Installation

install.packages("margEVT")

# install.packages("devtools")
devtools::install_github("rodrigosqrt3/margEVT")

4. Quick-Start Example

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)

margEVT: Regularized Point Processes and Stochastic Marginalization for Return Level Inference