Package {margEVT}

Type:	Package
Title:	Regularized Point Processes and Stochastic Marginalization for Extremes
Version:	0.4.0
Description:	Implements a non-stationary extreme value analysis framework by coupling a covariate-driven non-homogeneous Poisson process (NHPP) with (adaptive) elastic-net regularization and analytical gradients. Provides methods for estimating conditional return levels and unconditional (marginalized) return levels via parametric stochastic integration over stable vector autoregressive VAR(p) covariate trajectories, or non-parametric annual-block resampling. Supports block-specific penalty controls, BIC selection along a regularization path, Karush-Kuhn-Tucker convergence checks, time-rescaling residual diagnostics, conditional parametric bootstrap inference, and walk-forward assessment.
License:	GPL (≥ 3)
Encoding:	UTF-8
URL:	https://github.com/rodrigosqrt3/margEVT
BugReports:	https://github.com/rodrigosqrt3/margEVT/issues
RoxygenNote:	7.3.3
Imports:	graphics, grDevices, parallel, stats, utils
Suggests:	knitr, rmarkdown, testthat (≥ 3.0.0), vars
VignetteBuilder:	knitr
Config/testthat/edition:	3
NeedsCompilation:	no
Packaged:	2026-10-13 17:17:13 UTC; rodri
Author:	Rodrigo Fonseca Villa [aut, cre]
Maintainer:	Rodrigo Fonseca Villa <rodrigo03.villa@gmail.com>
Repository:	CRAN
Date/Publication:	2026-10-13 23:00:15 UTC

Extract names of non-seasonal covariates from a fitted model

Description

Returns the covariate names that appear in the model with a non-negligible coefficient - i.e. those that actually affect the parameters and need to be supplied during marginalization. Seasonal harmonic columns (cos<k>, sin<k>, sen<k>) are always excluded.

Usage

active_covariates(fit, free_cols = NULL, tol = NULL)

Arguments

Value

Character vector of active covariate names, or NULL if the model has no active non-seasonal covariates.

Examples

set.seed(1)
df <- data.frame(x = rnorm(2000), z = rnorm(2000))
df$y <- 10 + 2 * df$x + rexp(2000, 1 / 3)
fit <- fit_nhpp(df, threshold = 18, loc_vars = c("x", "z"),
                penalty = "lasso", lambda = 5, verbose = FALSE)
active_covariates(fit)

Walk-forward backtesting for NHPP extreme value models

Description

Evaluates predictive performance of a fitted nhpp_fit model via walk-forward cross-validation. For each window, the model is refit on training data at the original penalty and lambda, return levels are computed, and exceedance rates are compared to nominal rates via binomial tests. A Kolmogorov-Smirnov comparison on tail-conditional PIT values is returned as an exploratory diagnostic; it is not a formal test of global predictive calibration. Validation years with fewer than min_obs_year observed values are excluded from binomial denominators.

Usage

backtest(
  fit,
  data,
  varname,
  TRs = c(2, 5, 10),
  n_obs = NULL,
  period = NULL,
  window_years = 5L,
  min_train_years = 10L,
  n_boot = 200L,
  interactions = list(),
  year_col = "year",
  min_obs_year = NULL,
  seed = 2024L,
  trend_col = NULL,
  fixed = list(),
  reselect_lambda = FALSE,
  cores = 1L,
  verbose = TRUE
)

Arguments

Value

A list with three elements:

See Also

marginalize, fit_nhpp

Examples

set.seed(1)
n_years <- 30; n_obs <- 50
df <- data.frame(x = rnorm(n_years * n_obs),
                 year = rep(1991:2020, each = n_obs))
df$y <- 5 + df$x + rexp(nrow(df), 1 / 2)
fit <- fit_nhpp(df, threshold = 10, loc_vars = "x", penalty = "none",
                obs_per_year = n_obs, verbose = FALSE)

bt <- backtest(fit, df, varname = "y", TRs = c(2, 5), n_obs = n_obs,
               window_years = 5, min_train_years = 15, n_boot = 50,
               verbose = FALSE)
bt$binom_tests

Compute BIC for a fitted nhpp_fit model

Description

Computes -2\ell + k \log m, where \ell is the unpenalized log-likelihood at the operational estimator, k the number of active parameters, and m the number of exceedances. Equivalent to BIC(fit) for fits produced by margEVT 0.4.0 or later.

Usage

bic_nhpp(fit, tol = NULL)

Arguments

Value

Numeric scalar. BIC value.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)), x = rnorm(330))
fit <- fit_nhpp(df, threshold = 4, loc_vars = "x", penalty = "none",
                verbose = FALSE)
bic_nhpp(fit)

Parametric bootstrap confidence intervals for model coefficients

Description

Simulates R bootstrap samples, refits at the original penalty, lambda, and penalty weights, and returns empirical percentile intervals for each operational coefficient. See bootstrap_rl for the resampling scheme and the random-number streams.

Usage

bootstrap_coef(
  fit,
  data,
  R = 200L,
  level = 0.95,
  seed = 42L,
  cores = 1L,
  verbose = TRUE
)

Arguments

Value

A data frame with columns parameter, estimate, CI_low, CI_high, n_ok, with one row for every coefficient in the fitted model. Intervals require at least five successful replicates. The replicate matrix is attached as attribute "replicates".

See Also

bootstrap_rl

Examples

set.seed(1)
df <- data.frame(x = rnorm(1500))
df$y <- 10 + df$x + rexp(1500, 1 / 3)
fit <- fit_nhpp(df, threshold = 18, loc_vars = "x", penalty = "none",
                verbose = FALSE)

bootstrap_coef(fit, df, R = 50, verbose = FALSE)

Parametric bootstrap confidence intervals for return levels

Description

Simulates R bootstrap samples from the fitted model, refits at the original model design, penalty, lambda, penalty weights, and activity tolerance (starting from the original estimate), recomputes return levels via marginalize, and returns empirical percentile intervals. Exceedance count and times are held fixed, so the result is a nominal conditional summary rather than unconditional or post-selection inference.

Usage

bootstrap_rl(
  fit,
  data,
  TRs = c(2, 5, 10, 20, 50, 100),
  R = 200L,
  approach = "C",
  marginalize_args = list(),
  level = 0.95,
  seed = 42L,
  cores = 1L,
  verbose = TRUE
)

Arguments

Details

Each replicate uses its own L'Ecuyer-CMRG random-number stream derived from seed, so the output does not depend on cores. The caller's random-number state is left untouched when seed is supplied.

Value

A data frame with columns TR, RL_est (the return level from the original fitted model under the requested approach), CI_low, CI_high, n_ok (number of successful replicates). The R x length(TRs) matrix of replicate return levels is attached as attribute "replicates".

See Also

bootstrap_coef, marginalize

Examples

set.seed(1)
df <- data.frame(x = rnorm(1500))
df$y <- 10 + df$x + rexp(1500, 1 / 3)
fit <- fit_nhpp(df, threshold = 18, loc_vars = "x", penalty = "none",
                verbose = FALSE)

bootstrap_rl(fit, df, TRs = c(10, 50), R = 50, approach = "A",
             verbose = FALSE)

Build a covariate data frame for marginalization

Description

Constructs a data frame with n_obs rows containing all columns required by a fitted nhpp_fit model. The default seasonal columns (cos1, sen1, cos2, sen2) and any other harmonic column required by the model (cos<k>, sin<k>, sen<k>) are computed automatically with period (see seasonal_harmonics). Any other required column is filled from cov_vals, or set to zero with a warning if absent.

Usage

build_cov_annual(
  fit,
  cov_vals = list(),
  n_obs = NULL,
  period = NULL,
  interactions = list()
)

Arguments

Value

A data frame with n_obs rows and all columns required by fit.

Examples

set.seed(1)
n <- 3650
df <- cbind(data.frame(x = rnorm(n)), seasonal_harmonics(rep(1:365, 10)))
df$y <- 10 + df$x + 3 * df$cos1 + rexp(n, 1 / 3)
fit <- fit_nhpp(df, threshold = 20, loc_vars = c("x", "cos1", "sen1"),
                free_vars = c("cos1", "sen1"), penalty = "none",
                verbose = FALSE)
head(build_cov_annual(fit, list(x = 0.5), n_obs = 365))

Build design matrices for the NHPP point process model

Description

Build design matrices for the NHPP point process model

Usage

build_design_matrices(
  df,
  loc_vars = NULL,
  scale_vars = NULL,
  shape_vars = NULL,
  free_vars = NULL
)

Arguments

Value

A list with elements:

Examples

df <- cbind(data.frame(y = rexp(10), x = rnorm(10)),
            seasonal_harmonics(1:10))
dm <- build_design_matrices(df, loc_vars = c("x", "cos1", "sen1"),
                            scale_vars = "x", free_vars = c("cos1", "sen1"))
colnames(dm$X_mu)
dm$idx_pen_mu

Extract coefficients of a fitted nhpp_fit model

Description

Extract coefficients of a fitted nhpp_fit model

Usage

## S3 method for class 'nhpp_fit'
coef(object, operational = FALSE, ...)

Arguments

Value

Named numeric vector of coefficients.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)), x = rnorm(330))
fit <- fit_nhpp(df, threshold = 4, loc_vars = "x", penalty = "none",
                verbose = FALSE)
coef(fit)

Fit a non-homogeneous point process model for extremes

Description

Fits a peaks-over-threshold point process model via penalized maximum likelihood with quasi-Newton optimization and exact analytical gradients. The location, log-scale, and shape parameters are linear in the supplied covariates; seasonal harmonics and other free_vars are never penalized.

Usage

fit_nhpp(
  df,
  threshold,
  loc_vars = NULL,
  scale_vars = NULL,
  shape_vars = NULL,
  free_vars = NULL,
  penalty = c("none", "lasso", "ridge", "elnet"),
  alpha = 0.5,
  lambda = "bic",
  penalize_shape = TRUE,
  obs_per_year = 365.25,
  active_tol = 0.01,
  lambda_scaling = c("common", "gradient"),
  penalty_weights = "none",
  adaptive_gamma = 1,
  lambda_grid = NULL,
  init = NULL,
  grad_tol = 5e-04,
  maxit = 10000L,
  calc_hessian = FALSE,
  verbose = TRUE
)

Arguments

Details

Convergence is assessed with the Karush-Kuhn-Tucker (KKT) conditions of the smoothed elastic-net objective: unpenalized and active coordinates must have a near-zero gradient, while penalized coordinates numerically at zero must satisfy the subgradient condition. The violation is scaled by the number of exceedances and compared with grad_tol.

When lambda = "bic", the regularization path is evaluated on a coarse grid followed by a refined grid around the provisional optimum (or on a user-supplied lambda_grid). The returned model is the path solution at the selected penalty, so the reported fit is exactly the one scored by the BIC.

Value

An object of class nhpp_fit. Besides the smooth optimizer solution (par) and the operational activity-thresholded solution (par_oper), it stores the penalty weights, the scaled KKT violation (kkt), and, for lambda = "bic", the evaluated regularization path (path).

See Also

marginalize, bootstrap_coef, plot.nhpp_fit

Examples

set.seed(42)
n <- 2000
df <- data.frame(x1 = rnorm(n), x2 = rnorm(n))
df$y <- 10 + 2 * df$x1 + rexp(n, 1 / 3)

# Unpenalized stationary-shape model
fit0 <- fit_nhpp(df, threshold = 18, loc_vars = c("x1", "x2"),
                 penalty = "none", verbose = FALSE)
coef(fit0)


# Adaptive LASSO with BIC-selected penalty
fit <- fit_nhpp(df, threshold = 18, loc_vars = c("x1", "x2"),
                penalty = "lasso", lambda = "bic",
                penalty_weights = "adaptive", verbose = FALSE)
summary(fit)

Fit a stationary VAR generator of covariate trajectories

Description

Deseasonalizes each base variable by OLS on n_harmonics Fourier harmonic pairs, standardizes the residuals, and fits a stable vector autoregression VAR(p) without intercept by least squares, with the order selected by the Schwarz criterion on a common estimation sample (the convention of vars::VARselect). A univariate base variable gives a stable AR(p) model.

Usage

fit_var_generator(
  fit,
  data,
  vars = NULL,
  period = NULL,
  lag_max = NULL,
  n_harmonics = 2L,
  derive = NULL,
  derive_lead = 0L
)

Arguments

Details

By default the base variables are the active covariates of fit. When model covariates are deterministic transformations of a smaller set of base variables (lags, rolling means, differences, interactions, standardizations), supply the base variables in vars and a derive function: simulate_covariates then simulates the base variables and recomputes the model covariates on each trajectory, preserving their deterministic relations.

Seasonal harmonics present in data (cos<k>, sen<k> or sin<k>) are reused; otherwise they are computed from the row index, which assumes contiguous rows starting at the first position of a year.

Value

An object of class nhpp_var_generator, including the coefficient matrix, innovation covariance, root moduli and spectral radius, or NULL (with a warning) if there are no variables to model. An unstable selected VAR causes an error.

Examples

set.seed(1)
n <- 3650
t <- rep(1:365, 10)
x <- as.numeric(stats::arima.sim(list(ar = 0.7), n)) + 2 * cos(2 * pi * t / 365.25)
df <- cbind(data.frame(x = x), seasonal_harmonics(t))
df$y <- 10 + df$x + rexp(n, 1 / 3)
fit <- fit_nhpp(df, threshold = 20, loc_vars = "x", penalty = "none",
                verbose = FALSE)
gen <- fit_var_generator(fit, df)
gen$p_opt
gen$spectral_radius

# Model covariates derived from a base variable: the lagged value of x
df$x_lag1 <- c(NA, head(df$x, -1))
gen_d <- fit_var_generator(
  fit, df[-1, ], vars = "x", derive_lead = 1,
  derive = function(b) data.frame(x = b$x, x_lag1 = c(NA, head(b$x, -1))))
head(simulate_covariates(gen_d, n_mc = 2, seed = 1)[[1]])

Check if an object is an nhpp_fit

Description

Check if an object is an nhpp_fit

Usage

is_nhpp_fit(x)

Arguments

Value

Logical.

Examples

is_nhpp_fit(list())

Compute unconditional return levels by marginalizing over covariates

Description

Given a fitted nhpp_fit model, computes return levels by integrating out covariate uncertainty under one or more approaches:

Usage

marginalize(
  fit,
  data,
  TRs = c(2, 5, 10, 20, 50, 100),
  n_obs = NULL,
  period = NULL,
  approaches = c("A", "C"),
  scenarios = NULL,
  mc_sample = NULL,
  interactions = list(),
  fixed = list(),
  year_col = "year",
  n_boot = 500L,
  seed = 2024L,
  tol_root = 1e-04,
  z_hi = NULL
)

Arguments

Details

A - conditional

Covariates fixed at user-specified scenarios (e.g. quantiles of observed distributions). Fast; useful for sensitivity analysis.

B - parametric

Covariates drawn from a user-supplied Monte Carlo sample of complete annual trajectories (e.g. from simulate_covariates).

C - empirical

Annual-block resampling of observed covariate years from data, preserving within-year temporal structure.

For each approach and return period T, the return level solves G(z) = 1 - 1/T, where G is the annual non-exceedance probability averaged over trajectories. For approaches B and C the conditional Monte Carlo standard error of the root, \hat s_T / (\sqrt{n}\,\hat g_T), is reported in se_mc, with \hat s_T the standard deviation of the pathwise probabilities at the root and \hat g_T the mean pathwise derivative with respect to the level. It measures only finite-ensemble integration error, conditional on the fitted model and trajectory law.

Value

A data frame of class nhpp_rl with columns approach, scenario, TR, RL, and se_mc (NA for approach A). See plot.nhpp_rl and rl_table.

Examples

set.seed(1)
n_years <- 30; n_obs <- 50
df <- data.frame(x = rnorm(n_years * n_obs),
                 year = rep(seq_len(n_years), each = n_obs))
df$y <- 10 + 2 * df$x + rexp(nrow(df), 1 / 3)
fit <- fit_nhpp(df, threshold = 18, loc_vars = "x", penalty = "none",
                obs_per_year = n_obs, verbose = FALSE)
rl <- marginalize(fit, df, TRs = c(5, 10, 50),
                  scenarios = list(low = list(x = -1), high = list(x = 1)),
                  n_boot = 200)
rl
rl_table(rl)

Extract the number of exceedances from a fitted model

Description

Extract the number of exceedances from a fitted model

Usage

n_exceedances(fit, y)

Arguments

Value

Integer. Number of observations exceeding the threshold.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)))
fit <- fit_nhpp(df, threshold = 4, penalty = "none", verbose = FALSE)
n_exceedances(fit)

Log-likelihood, number of observations, and prediction for nhpp_fit

Description

logLik() returns the unpenalized log-likelihood at the operational estimator, with degrees of freedom equal to the number of active parameters and nobs equal to the number of exceedances, so that AIC() and BIC() work; BIC(fit) equals bic_nhpp(fit). predict() is a thin wrapper around predict_params.

Usage

## S3 method for class 'nhpp_fit'
logLik(object, ...)

## S3 method for class 'nhpp_fit'
nobs(object, ...)

## S3 method for class 'nhpp_fit'
predict(object, newdata = NULL, operational = TRUE, ...)

Arguments

Value

logLik(): an object of class logLik. nobs(): the number of exceedances. predict(): a list with components mu, sigma, xi.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)), x = rnorm(330))
fit <- fit_nhpp(df, threshold = 4, loc_vars = "x", penalty = "none",
                verbose = FALSE)
logLik(fit)
c(AIC = AIC(fit), BIC = BIC(fit), nobs = nobs(fit))
str(predict(fit, newdata = data.frame(x = c(-1, 0, 1))))

Plot diagnostics for a fitted nhpp_fit object

Description

Plot diagnostics for a fitted nhpp_fit object

Usage

## S3 method for class 'nhpp_fit'
plot(
  x,
  type = c("intensity", "fitted", "path", "qq"),
  residuals = c("interarrival", "magnitude"),
  ...
)

Arguments

Value

Invisibly returns NULL, called for side effects.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)), x = rnorm(330))
fit <- fit_nhpp(df, threshold = 4, loc_vars = "x", penalty = "none",
                verbose = FALSE)
plot(fit, type = "intensity")
plot(fit, type = "qq")

Plot return-level curves from marginalize()

Description

Draws return level against return period (logarithmic axis), one curve per approach and scenario. For approaches with a Monte Carlo standard error, dashed lines show \pm 2 standard errors.

Usage

## S3 method for class 'nhpp_rl'
plot(x, ...)

Arguments

Value

Invisibly returns NULL, called for side effects.

Examples

set.seed(1)
df <- data.frame(x = rnorm(1500), year = rep(1:30, each = 50))
df$y <- 10 + 2 * df$x + rexp(1500, 1 / 3)
fit <- fit_nhpp(df, threshold = 18, loc_vars = "x", penalty = "none",
                obs_per_year = 50, verbose = FALSE)
plot(marginalize(fit, df, TRs = c(2, 5, 10, 50, 100), n_boot = 200))

Analytical gradient of the penalized negative log-likelihood

Description

Analytical gradient of the penalized negative log-likelihood

Usage

pp_grad(
  par,
  dm,
  y,
  threshold,
  lambda = 0.1,
  alpha = 0.5,
  pen_xi = TRUE,
  eps_l1 = 1e-10,
  obs_per_year = 365.25,
  penalty_weights = NULL
)

Arguments

Value

Numeric vector of length length(par).

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)), x = rnorm(330))
dm <- build_design_matrices(df, loc_vars = "x")
pp_grad(c(6, 0.1, log(1.5), 0.1), dm, df$y, threshold = 4, lambda = 0)

Penalized negative log-likelihood for the NHPP model

Description

Evaluates the negative log-likelihood of the non-homogeneous Poisson process (point-process) representation of threshold exceedances, plus a smoothed elastic-net penalty on the penalized coefficients.

Usage

pp_nllh(
  par,
  dm,
  y,
  threshold,
  lambda = 0.1,
  alpha = 0.5,
  pen_xi = TRUE,
  eps_l1 = 1e-10,
  obs_per_year = 365.25,
  penalty_weights = NULL
)

Arguments

Value

Scalar. The penalized negative log-likelihood value. Parameter vectors outside the point-process support return 1e9.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)), x = rnorm(330))
dm <- build_design_matrices(df, loc_vars = "x")
par <- c(6, 0.1, log(1.5), 0.1)
pp_nllh(par, dm, df$y, threshold = 4, lambda = 0)
pp_nllh(par, dm, df$y, threshold = 4, lambda = 1)

Compute time-varying GEV parameters from a fitted model

Description

Given an nhpp_fit object and a data frame of covariate values, returns the implied location, scale, and shape parameters at each row.

Usage

predict_params(fit, newdata = NULL, operational = TRUE)

Arguments

Value

A list with three numeric vectors of length nrow(newdata):

Examples

set.seed(1)
df <- data.frame(x = rnorm(2000))
df$y <- 10 + 2 * df$x + rexp(2000, 1 / 3)
fit <- fit_nhpp(df, threshold = 18, loc_vars = "x", penalty = "none",
                verbose = FALSE)
predict_params(fit, data.frame(x = c(-1, 0, 1)))

Print a fitted nhpp_fit model

Description

Print a fitted nhpp_fit model

Usage

## S3 method for class 'nhpp_fit'
print(x, ...)

Arguments

Value

Invisibly returns x.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)), x = rnorm(330))
print(fit_nhpp(df, threshold = 4, loc_vars = "x", penalty = "none",
               verbose = FALSE))

Print a summary of a fitted nhpp_fit model

Description

Print a summary of a fitted nhpp_fit model

Usage

## S3 method for class 'summary.nhpp_fit'
print(x, digits = 5L, ...)

Arguments

Value

Invisibly returns x.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)))
print(summary(fit_nhpp(df, threshold = 4, penalty = "none",
                       verbose = FALSE)))

Transformed residuals of a fitted nhpp_fit model

Description

Two residual types are available, both standard exponential under a correctly specified model:

"interarrival"

Time-rescaled inter-exceedance increments Z_k = n_y^{-1}\sum_{j: t_{k-1} \le j < t_k} \hat\tau_u(t_j), k = 2, \ldots, m, where \hat\tau_u(t) is the fitted exceedance rate at the threshold (Ogata, 1988; Smith and Shively, 1995).

"magnitude"

Exponential transforms of the threshold excesses under the implied generalized Pareto law, \xi^{-1}\log\{1 + \xi (y - u)/\beta(t)\} with \beta(t) = \sigma(t) + \xi(t)\{u - \mu(t)\}.

Usage

## S3 method for class 'nhpp_fit'
residuals(
  object,
  type = c("interarrival", "magnitude"),
  operational = TRUE,
  ...
)

Arguments

Value

Numeric vector of residuals.

References

Ogata, Y. (1988). Statistical models for earthquake occurrences and residual analysis for point processes. Journal of the American Statistical Association, 83(401), 9–27.

Smith, R. L. and Shively, T. S. (1995). Point process approach to modeling trends in tropospheric ozone based on exceedances of a high threshold. Atmospheric Environment, 29(23), 3489–3499.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)), x = rnorm(330))
fit <- fit_nhpp(df, threshold = 4, loc_vars = "x", penalty = "none",
                verbose = FALSE)
z <- residuals(fit)
c(mean = mean(z), var = var(z))

Wide-format return level table from marginalize() output

Description

Converts the long-format output of marginalize into a wide table with one row per scenario and one column per return period.

Usage

rl_table(marg_result)

Arguments

Value

A wide data frame.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)),
                 year = rep(1991:2020, each = 11))
fit <- fit_nhpp(df, threshold = 4, penalty = "none", obs_per_year = 11,
                verbose = FALSE)
rl_table(marginalize(fit, df, TRs = c(10, 50), approaches = "A"))

Seasonal Fourier harmonics

Description

Builds the seasonal harmonic columns used by margEVT. The k-th pair is \cos(2\pi k t / P) and \sin(2\pi k t / P), where t is the position of the observation within the year and P is the seasonal period. Using this function to build the fitting data guarantees that build_cov_annual reconstructs the same seasonal cycle during marginalization, where positions run from 1 to n_obs within each simulated or resampled year.

Usage

seasonal_harmonics(
  t,
  period = 365.25,
  n_harmonics = 2L,
  sine = c("sen", "sin")
)

Arguments

Value

A data frame with columns cos1, sen1 (or sin1), ..., one row per element of t.

Examples

h <- seasonal_harmonics(1:365, period = 365.25)
head(h)
names(seasonal_harmonics(1:52, period = 52, n_harmonics = 3, sine = "sin"))

Simulate stationary covariate trajectories from a VAR generator

Description

Simulates n_mc independent annual trajectories. All trajectories are propagated simultaneously, so the cost grows with the trajectory length rather than with n_mc times the length. Each trajectory is initialized from the stationary Gaussian law of the companion state (default) or from zero followed by burn_in transitions. Seasonality is re-added at positions 1, ..., n_obs; if the generator has a derive function, model covariates are recomputed from the simulated base variables, using derive_lead additional leading observations that are discarded afterwards.

Usage

simulate_covariates(
  generator,
  n_mc,
  n_obs = NULL,
  burn_in = 300L,
  initialization = c("stationary", "burnin"),
  seed = NULL
)

Arguments

Value

A list of length n_mc, each element an n_obs-row data frame of covariates, directly usable as mc_sample in marginalize (approach B).

Examples

set.seed(1)
n <- 3650
t <- rep(1:365, 10)
df <- cbind(data.frame(x = as.numeric(stats::arima.sim(list(ar = 0.6), n))),
            seasonal_harmonics(t))
df$y <- 10 + df$x + rexp(n, 1 / 3)
fit <- fit_nhpp(df, threshold = 20, loc_vars = "x", penalty = "none",
                verbose = FALSE)
mc <- simulate_covariates(fit_var_generator(fit, df), n_mc = 3, seed = 1)
length(mc); dim(mc[[1]])

Summarise a fitted nhpp_fit model

Description

Builds a structured summary of a fitted nhpp_fit object, including the penalty used, the active and shrunk coefficients, information criteria, convergence diagnostics, and (when available) Hessian conditioning. Coefficient standard errors are deliberately not reported: after penalization and data-driven selection they lack a valid interpretation (see bootstrap_coef for conditional intervals).

Usage

## S3 method for class 'nhpp_fit'
summary(object, tol = NULL, ...)

Arguments

Value

An object of class summary.nhpp_fit, printed by print.summary.nhpp_fit.

Examples

set.seed(1)
df <- data.frame(y = c(rexp(300), runif(30, 5, 10)), x = rnorm(330))
fit <- fit_nhpp(df, threshold = 4, loc_vars = "x", penalty = "none",
                verbose = FALSE)
summary(fit)