Covariate-driven extremes with margEVT

Model

Exceedances of a threshold $u$ by a series $Y_t$ are modelled as a point process with integrated intensity \[ \Lambda_t\{(y, \infty)\} = \frac{1}{n_y} \left[1 + \xi_t \frac{y - \mu_t}{\sigma_t}\right]_+^{-1/\xi_t}, \qquad y > u, \] where $n_y$ is the number of observations per year, so that $\mu_t$, $\sigma_t$ and $\xi_t$ are the parameters of the generalized extreme value (GEV) distribution of the annual maximum. Each parameter is linear in covariates ($\log\sigma_t$ for the scale), and the covariate coefficients are estimated by penalized maximum likelihood. With penalty_weights = "adaptive" the penalty is the adaptive LASSO, with weights $1/|\tilde\theta_k|^\gamma$ from an unpenalized fit. The shape is kept at $\xi_t \ge -1/2$ (shape_lower), the region where maximum likelihood is regular.

Given one year of covariates $\mathbf{v}$, the annual non-exceedance probability is $G(z \mid \mathbf{v}) = \exp\{-\sum_t \Lambda_t(z)\}$. The marginal return level of period $T$ solves \[ G_\Pi(z_T) = \mathrm{E}_\Pi\{G(z_T \mid \mathbf{V})\} = 1 - 1/T , \] where $\Pi$ is a law for whole years of covariates. margEVT offers three choices: fixed scenarios (approach A, which gives conditional levels), a fitted vector autoregressive generator (approach B), and resampled observed years (approach C).

Simulated data

We simulate 40 years of weekly data with an annual cycle, an autocorrelated climate index x that drives the location, and an irrelevant covariate z.

set.seed(2024)
n_years <- 40; n_y <- 52
n <- n_years * n_y
week <- rep(seq_len(n_y), n_years)
df <- cbind(
  data.frame(year = rep(1981:2020, each = n_y),
             x = as.numeric(arima.sim(list(ar = 0.8), n)) * 0.6,
             z = rnorm(n)),
  seasonal_harmonics(week, period = n_y)
)
df$y <- 20 + 3 * df$x + 4 * df$cos1 + 3 * rexp(n)
u <- quantile(df$y, 0.9)

seasonal_harmonics() builds the cos1, sen1, … columns with the same convention used later, when annual covariate blocks are rebuilt.

Exceedances of an autocorrelated series come in clusters, which the Poisson model does not allow for. extremal_index() estimates the extremal index $\theta$, the reciprocal of the mean cluster size (1 means no clustering):

extremal_index(df$y, u)
#>       theta  n_clusters  run_length 
#>   0.6645427 121.0000000   2.0000000

With $\hat\theta$ well below 1, we keep only the largest value of each cluster. decluster() sets the other exceedances to the threshold, so the result can be fitted directly; the raw series is kept for validation.

df$y_raw <- df$y
df$y <- decluster(df$y_raw, u)
c(exceedances = sum(df$y_raw > u), clusters = sum(df$y > u))
#> exceedances    clusters 
#>         208         121

Estimation

Harmonics are declared in free_vars, so they are never penalized. The penalty level is selected by the Bayesian information criterion (BIC) along a regularization path. With refit = TRUE, the selected model is refitted without penalty (the relaxed LASSO): the penalty only selects covariates, and the reported coefficients are free of shrinkage.

fit <- fit_nhpp(df, threshold = u,
                loc_vars   = c("cos1", "sen1", "x", "z"),
                scale_vars = c("x", "z"),
                free_vars  = c("cos1", "sen1"),
                penalty = "lasso", lambda = "bic",
                penalty_weights = "adaptive", refit = TRUE,
                obs_per_year = n_y, verbose = FALSE)
fit
#> -- nhpp_fit --------------------------------------
#>   Threshold    : 29.11
#>   Penalty      : lasso
#>   Lambda       : mu=4.8454, sigma=4.8454, xi=4.8454
#>   Weights      : adaptive: mle
#>   Alpha        : 1
#>   Converged    : TRUE
#>   Refit        : unpenalized on the selected set
#>   nllh (raw)   : 215.8259
#>   nllh (pen)   : 220.9725
#>   obs/year     : 52.00
#>   Active params: 6 of 9
#>   Coefficients (non-zero):
#>    mu.(Intercept)           mu.cos1           mu.sen1              mu.x 
#>          31.54392           3.76378           0.42864           2.60607 
#> sigma.(Intercept)    xi.(Intercept) 
#>           1.18793          -0.14583

summary(fit)
#> -- nhpp_fit summary --------------------------------------
#>   Threshold      : 29.11
#>   Penalty        : lasso
#>   Lambda         : mu=4.8454, sigma=4.8454, xi=4.8454
#>   Weights        : adaptive: mle
#>   Alpha          : 1
#>   Penalize shape : TRUE
#>   obs/year       : 52.00
#>   Converged      : TRUE
#>   KKT violation  : 6.36e-07 (per exceedance)
#>   nllh (raw)     : 215.8259
#>   nllh (pen)     : 220.9725
#>   Exceedances    : 121 | BIC: 460.426
#> 
#>   Active coefficients (6 of 9):
#>                   Estimate
#> mu.(Intercept)    31.54392
#> mu.cos1            3.76378
#> mu.sen1            0.42864
#> mu.x               2.60607
#> sigma.(Intercept)  1.18793
#> xi.(Intercept)    -0.14583
#> 
#>   Shrunk to zero : mu.z, sigma.x, sigma.z
active_covariates(fit)
#> [1] "x"

The irrelevant covariate z is dropped. The solution is checked against the Karush-Kuhn-Tucker conditions of the penalized problem (fit$kkt, the largest violation per exceedance). The regularization path and the BIC curve are stored in fit$path:

plot(fit, type = "path")

Diagnostics

By the time-rescaling theorem, the interarrival residuals are standard exponential under a correct model; the magnitude residuals check the fitted generalized Pareto tail. The Q-Q plots show exact pointwise 95% bands.

op <- par(mfrow = c(1, 2))
plot(fit, type = "qq", residuals = "interarrival")
plot(fit, type = "qq", residuals = "magnitude")

par(op)

The usual generics work: coef(), logLik(), AIC(), BIC(), nobs(), predict(), fitted(), residuals(), simulate() and update(). vcov() and confint() give Wald intervals, conditional on the selected model:

logLik(fit)
#> 'log Lik.' -215.8259 (df=6)
confint(fit)
#>                        2.5 %      97.5 %
#> mu.(Intercept)    30.6090664 32.47878242
#> mu.cos1            2.5049683  5.02258394
#> mu.sen1           -0.6077105  1.46499670
#> mu.x               1.7520137  3.46012005
#> sigma.(Intercept)  1.0275074  1.34835410
#> xi.(Intercept)    -0.2643596 -0.02730657
predict(fit, newdata = data.frame(cos1 = 1, sen1 = 0, x = c(-1, 0, 1),
                                  z = 0))
#> $mu
#> [1] 32.70163 35.30770 37.91377
#> 
#> $sigma
#> [1] 3.280286 3.280286 3.280286
#> 
#> $xi
#> [1] -0.1458331 -0.1458331 -0.1458331

Marginal return levels

Approach B needs a model for whole years of covariates. fit_var_generator() removes the seasonal cycle, fits a VAR(p) with the order chosen by the Schwarz criterion, and checks that it is stable; simulate_covariates() then draws stationary annual trajectories.

gen <- fit_var_generator(fit, df, vars = "x", lag_max = 4)
gen
#> VAR(1) covariate generator
#>   Covariates      : x
#>   Seasonal model  : 2 harmonics, period 52
#>   Spectral radius : 0.7902 (stable)
sims <- simulate_covariates(gen, n_mc = 300, seed = 1)

marginalize() computes the levels under the three approaches. Monte Carlo standard errors of the roots are reported for B and C.

rl <- marginalize(fit, df, TRs = c(2, 10, 50, 100),
                  approaches = c("A", "B", "C"),
                  scenarios = list(neutral = list(x = 0),
                                   positive = list(x = 1)),
                  mc_sample = sims, n_boot = 300, seed = 1)
rl_table(rl)
#>   approach            scenario      T2     T10     T50    T100
#> 1        A             neutral 33.6514 39.0272 42.7479 44.0932
#> 2        A            positive 36.2574 41.6333 45.3539 46.6993
#> 3        B       parametric_mc 34.3762 40.0934 44.0829 45.5499
#> 4        C empirical_bootstrap 34.4502 40.1142 44.0686 45.5218
plot(rl)

A level computed at a fixed covariate value answers a different question from the marginal one: it holds the climate index at that value for the whole year, while the marginal level averages over the years the index can produce. When the covariates also act on the scale, the conditional level at the median state tends to fall below the marginal one at long return periods, and longer records do not remove the difference, because it comes from the definition of the level, not from estimation error.

Risk over a planning horizon

A return period describes one year. For a structure designed to last $L$ years, horizon_risk() gives the probability that a level is exceeded at least once in the horizon, the expected number of years with an exceedance, and the expected waiting time. The trajectories can be any annual covariate blocks: here, 10 years resampled from the record, 200 times.

years <- unique(df$year)
traj <- replicate(200, lapply(sample(years, 10, replace = TRUE),
                              function(yr) df[df$year == yr, "x", drop = FALSE]),
                  simplify = FALSE)
hr <- horizon_risk(fit, z = c(40, 45, 50), trajectories = traj)
hr
#> Exceedance risk over a horizon of 10 years (200 trajectories)
#>   z prob_exceed    se_mc expected_years  wait
#>  40    0.665453 0.003424        1.03698 6.405
#>  45    0.121580 0.001503        0.12878 9.437
#>  50    0.007744 0.000177        0.00777 9.965
plot(hr)

design_life_level() inverts this: it gives the level exceeded with a given probability during the horizon (Rootzén and Katz, 2013).

design_life_level(fit, traj, prob = c(0.5, 0.1))
#>   prob horizon     DLL      se_mc
#> 1  0.5      10 41.1830 0.02512763
#> 2  0.1      10 45.4223 0.02757573

Non-stationary scenarios, such as covariates from climate projections or a trend that keeps growing, enter in the same way, as trajectories whose years differ.

Uncertainty

bootstrap_nhpp() refits the selected specification (same penalty level and weights) on parametric bootstrap samples, and returns percentile intervals for the return levels and for every coefficient in one pass. Each replicate uses its own L’Ecuyer-CMRG random-number stream, so the results do not depend on the number of cores.

b <- bootstrap_nhpp(fit, df, TRs = c(10, 50), R = 20, approach = "C",
                    marginalize_args = list(n_boot = 100, seed = 1),
                    verbose = FALSE)
b
#> Conditional parametric bootstrap | R = 20 | 95% percentile intervals | occurrence: fixed
#> 
#> Return levels (approach C):
#>  TR RL_est CI_low CI_high n_ok
#>  10 40.168 38.348  41.854   20
#>  50 44.134 41.952  47.242   20
#> 
#> Coefficients:
#>          parameter estimate   CI_low  CI_high p_zero n_ok
#>     mu.(Intercept) 31.54392 31.11676 31.88356      0   20
#>            mu.cos1  3.76378  3.07237  4.21117      0   20
#>            mu.sen1  0.42864  0.20407  0.47338      0   20
#>               mu.x  2.60607  2.19865  3.09588      0   20
#>               mu.z  0.00000  0.00000  0.00000      1   20
#>  sigma.(Intercept)  1.18793  1.03276  1.35426      0   20
#>            sigma.x  0.00000  0.00000  0.00000      1   20
#>            sigma.z  0.00000  0.00000  0.00000      1   20
#>     xi.(Intercept) -0.14583 -0.22210 -0.04628      0   20

The intervals are conditional on the selected model and do not include the uncertainty of the selection. In applications, use several hundred replicates and cores > 1.

Walk-forward validation

backtest() refits the model on expanding training windows, computes return levels from the training years only, and compares the observed annual maxima of the following years (of the raw series, varname) with the nominal exceedance rates. Declustering is redone inside each window, on the training years only, and so are the adaptive weights: no validation data enter the training fit. reselect_lambda = TRUE also re-selects the penalty, and trend_col handles a linear time trend.

bt <- backtest(fit, df, varname = "y_raw", TRs = c(2, 5),
               window_years = 10, min_train_years = 20, n_boot = 100,
               decluster = function(y) decluster(y, u), verbose = FALSE)
bt
#> Walk-forward backtest: 20 validation years
#> 
#> Exceedance rates (binomial comparison with 1/T):
#>  approach TR expected observed n_exceedances n_total p_value not_rejected
#>         C  2      0.5     0.45             9      20  0.8238         TRUE
#>         C  5      0.2     0.10             2      20  0.4019         TRUE
#> 
#> Tail-conditional PIT, Kolmogorov-Smirnov (exploratory):
#>  approach  n D_statistic p_value not_rejected
#>         C 20      0.2428  0.1598         TRUE