Extended Dynamic Quantile Linear Models

Bayesian quantile-regression routines for dynamic state-space models and static regression under the extended asymmetric Laplace (exAL) error distribution. The dynamic state-space models are extended dynamic quantile linear models (exDQLMs). The package combines dynamic exDQLM inference via Laplace-delta variational Bayes (LDVB), Markov chain Monte Carlo (MCMC), and legacy importance-sampling variational Bayes (ISVB) with static exAL regression via LDVB and MCMC, reduced asymmetric Laplace/dynamic quantile linear model (AL/DQLM) paths through fixed skewness, component builders for trend/seasonality/regression blocks, static shrinkage priors including ridge, regularized horseshoe, and 'rhs_ns', evidence lower bound (ELBO) diagnostics, optional C++ accelerators, and posterior predictive synthesis across separately fitted quantiles through 'quantileSynthesis()'. Dynamic exDQLM methods are described in Barata et al. (2020) .

exdqlm — Extended Dynamic Quantile Linear Models

exdqlm is a Bayesian quantile-regression package that combines dynamic state-space quantile models with static extended asymmetric Laplace (exAL) regression through a common application programming interface (API). It is built for problems where quantiles, rather than means, are the main object of interest, but the user still wants familiar state-space/model-matrix inputs and explicit posterior inference.

The current release provides:

• extended dynamic quantile linear model (exDQLM) inference via Laplace-delta variational Bayes (LDVB) and Markov chain Monte Carlo (MCMC), with legacy importance-sampling variational Bayes (ISVB) retained for backward compatibility and historical comparisons
• static Bayesian exAL regression with the same LDVB and MCMC engines
• modular model builders for trend, seasonality, and regression components
• reduced asymmetric Laplace (AL)/dynamic quantile linear model (DQLM) paths through dqlm.ind = TRUE (dynamic and static), with a static convenience alias al.ind = TRUE
• standardized variational Bayes (VB) trace diagnostics at fit$diagnostics$vb_trace for the evidence lower bound (ELBO), sigma, gamma, and convergence deltas
• structured exAL scale-skewness inference: LDVB uses a q(gamma) q(sigma | gamma) block by default, and MCMC uses an exact scale-collapsed gamma slice transition by default
• static shrinkage priors beta_prior = "ridge", "rhs", and "rhs_ns"
• post hoc posterior-predictive synthesis via quantileSynthesis()

The package collects Bayesian dynamic quantile state-space modeling, static Bayesian quantile regression, multiple inference engines, shrinkage priors for static coefficients, and post hoc synthesis across quantiles in one interface.

Terminology. The exAL family generalizes the standard AL by adding a skewness parameter, allowing for asymmetric tails. The standard AL is a special case with zero skewness. We refer to the generalized AL from Kotz et al. as Kotz-GAL to avoid confusion.

Installation

CRAN:

install.packages("exdqlm")

Development (GitHub):

# install.packages("pak")
pak::pak("AntonioAPDL/exdqlm")

Current release note

Version 1.1.2 is a narrow reproducibility update to the fast C++ dynamic MCMC backend. The FFBS state simulation step now uses a Cholesky covariance square root for stochastic multivariate-normal state draws, avoiding platform-dependent SVD bases under fixed seeds. The exported API and model specification workflow are unchanged.

Why exdqlm is distinctive

• Dynamic Bayesian quantile state-space modeling is the core use case: exdqlmLDVB() is the main VB path, exdqlmMCMC() provides posterior simulation, and legacy exdqlmISVB() remains available for backward compatibility and historical comparisons.
• Static Bayesian quantile regression is now part of the same package via exalStaticLDVB() and exalStaticMCMC(), rather than living in a separate code path or companion repository.
• Multiple inference engines are available depending on the use case: structured Laplace-delta variational Bayes (LDVB) as the main dynamic VB engine, posterior simulation (MCMC), and legacy fast approximate dynamic VB (ISVB) when older workflows need to be reproduced.
• User-facing VB diagnostics are standardized through fit$diagnostics$vb_trace, so plotting and monitoring code can use the same iteration-wise API across VB engines.
• Static shrinkage priors go beyond ridge: rhs and rhs_ns provide a horseshoe-family regularization story for sparse or weakly identified static coefficient problems.
• Post hoc multi-quantile synthesis is built in through quantileSynthesis(), which combines separately fitted quantile models into a unified posterior predictive distribution using isotonic correction and optional monotone rearrangement.

Workflow map

Object system

exdqlm uses ordinary S3 classes for fitted models and post-processing objects. Dynamic fits keep their engine-specific first class (exdqlmLDVB, exdqlmMCMC, or legacy exdqlmISVB) and also inherit from the shared exdqlmFit family. Static fits similarly inherit from exalStaticFit. Diagnostic, forecast, synthesis, and static-diagnostic functions return explicit objects with print(), summary(), and plot() methods where appropriate. Dynamic fits also support predict() for forecasting, and fitted/forecast/static objects support diagnostics() where diagnostic summaries are defined. Named helpers such as exdqlmPlot(), compPlot(), exdqlmForecast(), exdqlmDiagnostics(), exdqlmForecastDiagnostics(), and exalStaticDiagnostics() remain available for explicit or backward-compatible workflows.

Which engine should I use?

Default warmup behavior

The package now applies a conservative automatic warmup profile for the most numerically delicate shared blocks, so ordinary users do not have to assemble nested warmup lists just to get a stable first fit.

• beta_prior = "rhs" and beta_prior = "rhs_ns" keep the package's shared tau warmup schedule on by default.
• exAL VB entry points (exalStaticLDVB(), exdqlmLDVB()) apply a light automatic warmup for the (sigma, gamma) block. For unrestricted exAL fits, the default VB scale-skewness block uses a structured q(gamma) q(sigma | gamma) approximation with conditional GIG moments for sigma; the previous two-dimensional Laplace-delta block remains available through exal_make_vb_sigmagam_control(factorization = "laplace_delta").
• exAL MCMC entry points (exalStaticMCMC(), exdqlmMCMC()) apply a light automatic (sigma, gamma) warmup and keep VB warm starts available for the harder cases where they help. Their default unrestricted exAL kernel is mh.proposal = "collapsed_slice", an exact-target update that slices gamma after integrating out sigma and then redraws sigma from its exact conditional distribution.
• Advanced warmup control remains available through vb_control and mcmc_control, but those controls are now intended as the override path, not the default user workflow.

In practice, the recommended workflow is:

1. fit the model with the default API;
2. inspect diagnostics if the fit still looks unstable;
3. only then override the warmup controls explicitly.

Quick start (≤ 10 lines)

Local-level model at a single quantile (the median) using the package's main dynamic VB engine. We keep the pure-R path for CRAN-style reproducibility and use the reduced DQLM path to keep the example small.

set.seed(1)
library(exdqlm)

T      <- 120
state  <- cumsum(rnorm(T, sd = 0.2))
y      <- state + rnorm(T, sd = 1.0)

model  <- list(FF = matrix(1), GG = matrix(1), m0 = 0, C0 = 100)
options(exdqlm.use_cpp_kf = FALSE, exdqlm.use_cpp_samplers = FALSE)

fit <- exdqlmLDVB(
  y = y, p0 = 0.5, model = model, df = 0.98, dim.df = 1,
  dqlm.ind = TRUE, sig.init = 1.0
)

tail(fit$diagnostics$elbo, 3)

For plotting or monitoring VB convergence, use the standardized trace table:

head(fit$diagnostics$vb_trace[, c("iter", "elbo", "sigma", "gamma")])

Core concepts (at a glance)

• State-space skeleton: design (FF) and evolution (GG) matrices with a prior for the state vector (m0, C0).
• Quantile of interest: p0 (e.g., 0.1, 0.5, 0.9).
• exAL errors: controlled by scale and skewness; ordinary workflows should usually let LDVB or MCMC update them, while fixed-parameter paths remain available for explicit baseline or compatibility runs.
• Discount factors: df and dim.df control evolution per block (e.g., trend vs seasonality).
• VB traces: fit$diagnostics$vb_trace provides a standardized iteration-by-iteration table for the evidence lower bound (ELBO), sigma, gamma, and convergence deltas across VB engines.
• ELBO: retained at fit$diagnostics$elbo and mirrored in fit$diagnostics$vb_trace$elbo. LDVB traces are deterministic conditional on the data and tuning inputs, up to numerical tolerance; legacy ISVB traces may also reflect Monte Carlo variability from the importance-sampling step.

Current release highlights

• Dynamic LDVB algorithm via exdqlmLDVB() as the main VB routine for exDQLM fitting, with legacy exdqlmISVB() retained for backward compatibility.
• Synthesis helper quantileSynthesis() to combine posterior predictive draws from separately fitted quantile models.
• Dynamic regression blocks via regMod() and static exAL regression via exalStaticLDVB() and exalStaticMCMC().
• Reduced AL/DQLM paths across dynamic and static APIs via dqlm.ind = TRUE, with al.ind = TRUE available as a static convenience alias.
• Static shrinkage priors in both static LDVB and MCMC via beta_prior = "ridge", "rhs", or "rhs_ns".
• Transfer-function helpers exdqlmTransferLDVB() and exdqlmTransferMCMC(), with legacy exdqlmTransferISVB() retained for backward compatibility, plus expanded static object generics (exalStaticMCMC, exalStaticLDVB).
• Standardized user-facing naming: the primary API now uses exalStatic..., exdqlmTransfer..., and quantileSynthesis(), while documented legacy ISVB entry points remain available for backward-compatible workflows.
• C++ backend controls retained as optional: Kalman bridge default TRUE; builders, samplers, and post-predictive C++ paths default FALSE.
• Standardized VB trace diagnostics via fit$diagnostics$vb_trace, giving plot-ready iteration histories for ELBO, sigma, gamma, and convergence deltas across VB fits.
• Deterministic dynamic diagnostics using a single primary Kullback--Leibler (KL) normality diagnostic, continuous ranked probability score (CRPS) through an integrated quantile-score approximation, and held-out forecast scoring through diagnostics(). Advanced KL sensitivity details are available under diagnostics$kl.details.
• Standard S3 post-processing methods for dynamic fits: plot(fit) draws the fitted quantile, plot(fit, type = "component", index = ...) and plot(fit, type = "state", index = ...) display component/state summaries, predict(fit, ...) returns an exdqlmForecast object, and diagnostics(fit) returns an exdqlmDiagnostic object. Forecast and static-fit diagnostics use the same generic where defined. The named helpers remain available for explicit workflows and backward compatibility.

For CI/CRAN-style runs, keep optional C++ builders/samplers/post-pred FALSE and set exdqlm.use_cpp_kf = FALSE for strict R-path reproducibility.

Runtime options (summary)

Set with:

options(exdqlm.use_cpp_kf = TRUE)
options(exdqlm.use_cpp_builders = FALSE)
options(exdqlm.use_cpp_samplers = TRUE)
options(exdqlm.use_cpp_postpred = FALSE)
options(exdqlm.use_cpp_mcmc = TRUE)
options(exdqlm.cpp_mcmc_mode = "fast")

Backend control (minimal):

• Force pure-R backend: set options(exdqlm.use_cpp_kf = FALSE, exdqlm.use_cpp_builders = FALSE).
• Keep builder calls explicit with backend = "R" or backend = "cpp" in polytrendMod() and seasMod().

Minimal examples (CRAN-safe)

1) Single-quantile fit on built-in data (tiny slice)

Trend + seasonality + one climate-index regressor. Note: FF for the regressor is 1 × T. Combine components pairwise.

data("BTflow", package = "exdqlm")
data("climateIndices", package = "exdqlm")

set.seed(2)
T <- 150
y <- log(BTflow[seq_len(T)])
bt_dates <- seq(as.Date("1987-01-01"), by = "month", length.out = T)
idx <- match(bt_dates, climateIndices$date)
x <- scale(climateIndices$noi[idx])[, 1]

trend.comp <- polytrendMod(order = 1, m0 = 0, C0 = 1)
seas.comp  <- seasMod(p = 12, h = 1, C0 = diag(1, 2))

# 1-d regressor block (explicit 1 x T design)
reg.comp <- list(m0 = 0, C0 = 1, FF = matrix(x, nrow = 1), GG = matrix(1))

# combine via +.exdqlm
reg.comp <- as.exdqlm(reg.comp)
model    <- trend.comp + seas.comp + reg.comp

# one discount per block: (trend, seasonal[2-d], reg)
df     <- c(1.00, 0.98, 1.00)
dim.df <- c(1,       2,   1)

options(exdqlm.use_cpp_kf = FALSE, exdqlm.use_cpp_samplers = FALSE)

fit <- exdqlmLDVB(
  y = y, p0 = 0.5, model = model,
  df = df, dim.df = dim.df,
  dqlm.ind = TRUE, sig.init = 0.2
)

# quick checks
tail(fit$diagnostics$elbo, 2)
dim(fit$theta.out$sm)  # state-dimension x time

2) exAL helper sanity check (CDF ↔ quantile)

set.seed(3)
x      <- seq(-2, 2, length.out = 5)
p0     <- 0.25
mu     <- 0
sigma  <- 1
gamma  <- 0.0

# CDF then invert with QF — should approximately return x
cdf_vals <- pexal(x,  p0 = p0, mu = mu, sigma = sigma, gamma = gamma)
x_back   <- qexal(cdf_vals, p0 = p0, mu = mu, sigma = sigma, gamma = gamma)

round(cbind(x, x_back), 4)
#>       x x_back
#> [1,] -2     -2
#> [2,] -1     -1
#> [3,]  0      0
#> [4,]  1      1
#> [5,]  2      2

# A few random draws
rexal(5, p0 = p0, mu = mu, sigma = sigma, gamma = gamma)
#> [1] -0.5296664  5.4402490  0.7934288  0.4376783  2.5354967

CRAN-safety. All examples set a seed, use tiny data, finish in a few seconds, and explicitly keep the pure-R path.

3) Static Bayesian regression with reduced AL and RHS-family priors

set.seed(4)
n <- 80
p <- 5
X <- matrix(rnorm(n * p), n, p)
beta <- c(1, -1, 0, 0, 0.5)
y <- as.numeric(X %*% beta + rnorm(n))

# Reduced AL fit (gamma fixed at zero)
fit_al <- exalStaticLDVB(
  y = y, X = X, p0 = 0.5,
  al.ind = TRUE,
  max_iter = 150, tol = 1e-4, verbose = FALSE
)

# exAL fit with regularized horseshoe prior on coefficients
fit_rhs <- exalStaticMCMC(
  y = y, X = X, p0 = 0.5,
  beta_prior = "rhs",
  n.burn = 200, n.mcmc = 200, thin = 1,
  trace.diagnostics = FALSE,
  verbose = FALSE
)

# exAL fit with rhs_ns controls (same API family, additive option)
fit_rhs_ns <- exalStaticMCMC(
  y = y, X = X, p0 = 0.5,
  beta_prior = "rhs_ns",
  beta_prior_controls = list(
    tau0 = 0.5,
    a_zeta = 2,
    b_zeta = 1,
    shrink_intercept = FALSE
  ),
  n.burn = 200, n.mcmc = 200, thin = 1,
  trace.diagnostics = FALSE,
  verbose = FALSE
)

fit_al$dqlm.ind
fit_rhs$beta_prior$type
fit_rhs_ns$beta_prior$type

4) Multi-quantile synthesis (conceptual sketch)

Fit several quantiles separately, then combine their posterior predictive draws into a single unified posterior predictive distribution.

p_grid <- c(0.1, 0.5, 0.9)
fits <- lapply(p_grid, function(tau) {
  exdqlmLDVB(
    y = y, p0 = tau, model = model,
    df = df, dim.df = dim.df,
    sig.init = 0.2, gam.init = 0
  )
})

draws <- lapply(fits, function(m) m$samp.post.pred)

syn <- quantileSynthesis(
  draws_list = draws,
  p = p_grid,
  T_expected = length(y)
)

names(syn)

FAQ / Troubleshooting

• It runs slowly. Use short series (≤ 200), fix scale/skewness, and keep discount factors near but below one (≈ 0.96–0.99). Enable C++ bridges only if your toolchain supports them.

• Which VB diagnostic object should I plot? Start with fit$diagnostics$vb_trace. It provides one iteration-wise table across VB engines; use engine-specific internals only when you need lower-level block traces.

• ELBO dips slightly—bug? Small downward blips in fit$diagnostics$vb_trace$elbo can occur from numerical tolerance, approximation updates, or, for legacy ISVB, Monte Carlo variability. Look for an overall upward trend; if not, simplify the model or adjust variance/discounts.

• OpenMP not available. That’s fine. It is optional. Everything runs serially; examples here use the pure-R path.

• Numerical stability tips. Avoid extremely tight C0; start with moderate priors (e.g., C0 around 1–100 for simple models), and fix scale/skewness for initial runs.

How to cite

Barata, R., Prado, R., & Sansó, B. (2022). Fast inference for time-varying quantiles via flexible dynamic models with application to the characterization of atmospheric rivers. Annals of Applied Statistics, 16(1), 247–271. https://doi.org/10.1214/21-AOAS1497

License

MIT © The authors. See LICENSE.

Getting help

Open an issue: https://github.com/AntonioAPDL/exdqlm/issues