Calculate GPD Return Levels from Per-Station Fits — calculate_gpd_return

Computes return levels from a GPD (Generalized Pareto Distribution) fit produced by mev::fit.gpd(). Uses the standard GPD return level formula with delta method confidence intervals.

The formula is: $$z_T = u + \frac{\sigma}{\xi}\left[(\lambda T)^{\xi} - 1\right]$$ where $u$ is the threshold, $\sigma$ is the scale, $\xi$ is the shape, $\lambda$ is the exceedance rate, and $T$ is the return period in years.

When shape is approximately zero, the exponential fallback is used: $$z_T = u + \sigma \log(\lambda T)$$

Usage

calculate_gpd_return_levels(
  gpd_fit,
  return_periods = c(1, 5, 10),
  n_obs_per_year = 8760,
  n_total = NULL,
  exceedance_rate = NULL,
  conf_level = 0.95
)

Arguments

gpd_fit: A list from the per-station GPD fitting targets, with elements: u (threshold), scale, shape, n_exceed, and optionally se_scale, se_shape. If it contains an error field, NA rows are returned.
return_periods: Numeric vector of return periods in years (default: c(1, 5, 10))
n_obs_per_year: Number of observations per year for exceedance rate calculation (default: 8760 for hourly data)
n_total: Total number of observations. If NULL, estimated from n_exceed and threshold percentile.
exceedance_rate: Pre-computed exceedance rate (lambda). If NULL, estimated as n_exceed / n_total.
conf_level: Confidence level for intervals (default: 0.95)

Value

Data frame with columns: return_period, return_level, lower, upper, threshold_value, method. Returns NA return levels if the fit has an error or missing parameters.

Examples

fit <- list(u = 5.0, scale = 1.2, shape = 0.1, n_exceed = 500)
calculate_gpd_return_levels(fit, return_periods = c(10, 50, 100))
#>   return_period return_level lower upper threshold_value method
#> 1            10     20.75420    NA    NA               5    GPD
#> 2            50     25.60060    NA    NA               5    GPD
#> 3           100     27.94046    NA    NA               5    GPD