R/utils.R

In exdqlm: Extended Dynamic Quantile Linear Models

#
log_g<-function(gam){	base::log(2)+stats::pnorm(-abs(gam),log=TRUE)+0.5*gam^2 }
L.fn<-function(p0){ stats::uniroot(function(gam) exp(log_g(gam))-(1-p0), c(-1000,0))$root }
U.fn<-function(p0){ stats::uniroot(function(gam) exp(log_g(gam))-p0, c(0,1000))$root }
p.fn<-function(p0,gam){ (p0-as.numeric(gam<0))/exp(log_g(gam))+as.numeric(gam<0)}
A.fn<-function(p0,gam){ temp.p = p.fn(p0,gam); return((1-2*temp.p)/(temp.p*(1-temp.p))) }
B.fn<-function(p0,gam){ temp.p = p.fn(p0,gam); return((2)/(temp.p*(1-temp.p))) }
C.fn<-function(p0,gam){ temp.p = p.fn(p0,gam); return((as.numeric(gam>0)-temp.p)^(-1)) }
# Internal helper: validate bounds and fall back to R reference if needed.
.gamma_bounds_ok_basic <- function(L, U) {
  if (!is.numeric(L) || !is.numeric(U) || length(L) != 1L || length(U) != 1L) return(FALSE)
  if (!is.finite(L) || !is.finite(U)) return(FALSE)
  if (L >= U) return(FALSE)
  # Gamma bounds should always straddle 0 for p0 in (0, 1).
  if (L > 0 || U < 0) return(FALSE)
  TRUE
}

.gamma_bounds_ok_cpp <- function(L, U, p0, tol_log = 1e-4) {
  if (!.gamma_bounds_ok_basic(L, U)) return(FALSE)
  log_target_L <- base::log1p(-p0)
  log_target_U <- base::log(p0)
  if (!is.finite(log_target_L) || !is.finite(log_target_U)) return(FALSE)

  # Validate against the defining equations on the log scale:
  #   g_gamma(L) = 1 - p0,  g_gamma(U) = p0
  logL <- log_g(L)
  logU <- log_g(U)
  if (!is.finite(logL) || !is.finite(logU)) return(FALSE)

  (abs(logL - log_target_L) <= tol_log) && (abs(logU - log_target_U) <= tol_log)
}

.gamma_bounds_ref <- function(p0) {
  c(L = L.fn(p0), U = U.fn(p0))
}

.gamma_bounds <- function(p0) {
  stopifnot(is.numeric(p0), length(p0) == 1L, is.finite(p0), p0 > 0, p0 < 1)
  out_cpp <- try(get_gamma_bounds_cpp(p0), silent = TRUE)
  if (!inherits(out_cpp, "try-error") && length(out_cpp) == 2L) {
    if (is.null(names(out_cpp))) names(out_cpp) <- c("L", "U")
    if (.gamma_bounds_ok_cpp(out_cpp[1], out_cpp[2], p0)) return(out_cpp)
  }
  out_ref <- try(.gamma_bounds_ref(p0), silent = TRUE)
  if (!inherits(out_ref, "try-error") && .gamma_bounds_ok_basic(out_ref[1], out_ref[2])) {
    return(out_ref)
  }
  stop("Unable to compute valid gamma bounds for p0 = ", p0)
}

.gig_b_floor <- function() {
  1e-10
}

.sample_gig_devroye_required <- function(n_samples, p, a, b_vec, context = "gig") {
  if (!exists("sample_gig_devroye_vector", mode = "function")) {
    stop(sprintf("%s requires sample_gig_devroye_vector(), but it is not available", context))
  }

  eps_gig <- sqrt(.Machine$double.eps)
  b_floor <- .gig_b_floor()
  p <- as.numeric(p)[1]
  a <- as.numeric(a)[1]
  b_vec <- as.numeric(b_vec)

  if (!is.finite(p)) {
    stop(sprintf("%s requires a finite lambda; got %.6g", context, p))
  }

  if (!is.finite(a) || a <= 0) a <- eps_gig
  b_vec[!is.finite(b_vec)] <- b_floor
  b_vec <- pmax(b_vec, b_floor)

  draws <- sample_gig_devroye_vector(
    as.integer(n_samples)[1],
    p = p,
    a = a,
    b_vec = b_vec
  )

  bad <- which(!is.finite(draws) | draws <= 0)
  if (length(bad)) {
    first <- bad[1]
    stop(sprintf("%s returned %d invalid draws (first index=%d, value=%.6g)",
                 context, length(bad), first, draws[first]))
  }

  pmax(draws, eps_gig)
}

.sample_gig_devroye_pairs_required <- function(n_samples, p, a_vec, b_vec, context = "gig") {
  if (!exists("sample_gig_devroye_pairs", mode = "function")) {
    stop(sprintf("%s requires sample_gig_devroye_pairs(), but it is not available", context))
  }

  eps_gig <- sqrt(.Machine$double.eps)
  b_floor <- .gig_b_floor()
  p <- as.numeric(p)[1]
  a_vec <- as.numeric(a_vec)
  b_vec <- as.numeric(b_vec)

  if (!is.finite(p)) {
    stop(sprintf("%s requires a finite lambda; got %.6g", context, p))
  }
  if (length(a_vec) != length(b_vec)) {
    stop(sprintf("%s requires a_vec and b_vec to have the same length", context))
  }

  a_vec[!is.finite(a_vec) | a_vec <= 0] <- eps_gig
  b_vec[!is.finite(b_vec)] <- b_floor
  b_vec <- pmax(b_vec, b_floor)

  draws <- sample_gig_devroye_pairs(
    as.integer(n_samples)[1],
    p = p,
    a_vec = a_vec,
    b_vec = b_vec
  )

  bad <- which(!is.finite(draws) | draws <= 0)
  if (length(bad)) {
    first <- bad[1]
    stop(sprintf("%s returned %d invalid draws (first index=%d, value=%.6g)",
                 context, length(bad), first, draws[first]))
  }

  pmax(draws, eps_gig)
}

.exdqlm_format_progress_value <- function(x) {
  if (is.null(x) || length(x) == 0L) return(NULL)
  if (length(x) > 1L) {
    return(paste(vapply(as.list(x), .exdqlm_format_progress_value, character(1)), collapse = ","))
  }
  if (is.logical(x)) {
    if (is.na(x)) return("NA")
    return(if (isTRUE(x)) "yes" else "no")
  }
  if (is.character(x)) {
    return(as.character(x)[1])
  }
  if (is.numeric(x)) {
    x <- as.numeric(x)[1]
    if (!is.finite(x)) return("NA")
    ax <- abs(x)
    if (ax >= 1e4 || (ax > 0 && ax < 1e-3)) return(sprintf("%.3e", x))
    if (ax >= 100) return(sprintf("%.2f", x))
    if (ax >= 1) return(sprintf("%.3f", x))
    return(sprintf("%.3g", x))
  }
  as.character(x)[1]
}

.exdqlm_progress <- function(label, ..., .verbose = TRUE) {
  if (!isTRUE(.verbose)) return(invisible(NULL))
  fields <- list(...)
  pieces <- character(0)
  for (nm in names(fields)) {
    val <- .exdqlm_format_progress_value(fields[[nm]])
    if (is.null(val)) next
    if (nzchar(nm)) {
      pieces <- c(pieces, sprintf("%s=%s", nm, val))
    } else {
      pieces <- c(pieces, val)
    }
  }
  cat(paste(c(label, pieces), collapse = " | "), "\n")
  utils::flush.console()
  try(flush(stdout()), silent = TRUE)
  invisible(NULL)
}

.exdqlm_static_model_label <- function(dqlm.ind) {
  if (isTRUE(dqlm.ind)) "AL special case (gamma = 0)" else "exAL"
}


.exdqlm_unwrap_fit_bundle <- function(obj, max_depth = 32L) {
  depth <- 0L
  fit <- obj
  normalized <- NULL
  meta <- NULL
  is_bundle <- function(x) {
    is.list(x) && all(c("fit", "normalized", "meta") %in% names(x))
  }
  while (is_bundle(fit) && depth < max_depth) {
    if (is.null(normalized) && !is.null(fit$normalized)) normalized <- fit$normalized
    if (is.null(meta) && !is.null(fit$meta)) meta <- fit$meta
    fit <- fit$fit
    depth <- depth + 1L
  }
  list(fit = fit, normalized = normalized, meta = meta, depth = depth)
}

.normalize_gamma_prior_trunc_t <- function(PriorGamma = NULL) {
  if (is.null(PriorGamma)) {
    PriorGamma <- list(m_gam = 0, s_gam = 1, df_gam = 1)
  } else {
    need <- c("m_gam", "s_gam", "df_gam")
    if (!is.list(PriorGamma) || any(is.na(match(need, names(PriorGamma))))) {
      stop("`PriorGamma` must be a list containing `m_gam`, `s_gam`, and `df_gam`")
    }
  }

  PriorGamma$m_gam <- as.numeric(PriorGamma$m_gam)[1]
  PriorGamma$s_gam <- as.numeric(PriorGamma$s_gam)[1]
  PriorGamma$df_gam <- as.numeric(PriorGamma$df_gam)[1]

  if (!is.finite(PriorGamma$m_gam)) stop("`PriorGamma$m_gam` must be finite")
  if (!is.finite(PriorGamma$s_gam) || PriorGamma$s_gam <= 0) stop("`PriorGamma$s_gam` must be > 0")
  if (!is.finite(PriorGamma$df_gam) || PriorGamma$df_gam <= 0) stop("`PriorGamma$df_gam` must be > 0")

  PriorGamma
}

.gamma_log_prior_trunc_t <- function(gamma, bounds, PriorGamma = NULL) {
  bounds <- as.numeric(bounds)
  if (length(bounds) != 2L || !all(is.finite(bounds)) || bounds[1] >= bounds[2]) {
    stop("`bounds` must be a finite length-2 vector with bounds[1] < bounds[2]")
  }
  prior <- .normalize_gamma_prior_trunc_t(PriorGamma)
  crch::dtt(
    gamma,
    location = prior$m_gam,
    scale = prior$s_gam,
    df = prior$df_gam,
    left = bounds[1],
    right = bounds[2],
    log = TRUE
  )
}

.gamma_prior_density_trunc_t <- function(gamma, bounds, PriorGamma = NULL, log = FALSE) {
  lp <- .gamma_log_prior_trunc_t(gamma, bounds = bounds, PriorGamma = PriorGamma)
  if (isTRUE(log)) lp else exp(lp)
}

.exdqlm_uni_slice_bounded <- function(
  x0,
  log_density,
  w = 0.1,
  m = Inf,
  lower = -Inf,
  upper = Inf,
  ...
) {
  if (!is.numeric(x0) || length(x0) != 1L || !is.finite(x0)) {
    stop("x0 must be a single finite numeric value.")
  }
  if (!is.function(log_density)) stop("log_density must be a function.")
  if (!is.numeric(w) || length(w) != 1L || !is.finite(w) || w <= 0) {
    stop("w must be a single positive finite numeric value.")
  }
  if (!is.numeric(m) || length(m) != 1L ||
      (!is.infinite(m) && (!is.finite(m) || m <= 0 || floor(m) != m))) {
    stop("m must be Inf or a single positive integer.")
  }
  if (!is.numeric(lower) || length(lower) != 1L || !is.finite(lower)) {
    stop("lower must be a single finite numeric value.")
  }
  if (!is.numeric(upper) || length(upper) != 1L || !is.finite(upper)) {
    stop("upper must be a single finite numeric value.")
  }
  if (lower >= upper) stop("lower must be strictly less than upper.")
  if (x0 < lower || x0 > upper) stop("x0 must lie within [lower, upper].")

  n_eval <- 0L
  logf <- function(x) {
    n_eval <<- n_eval + 1L
    as.numeric(log_density(x, ...))[1]
  }

  gx0 <- logf(x0)
  if (!is.finite(gx0)) {
    stop("log_density(x0) must be finite for slice sampling.")
  }

  logy <- gx0 - stats::rexp(1)
  u <- stats::runif(1, 0, w)
  L <- x0 - u
  R <- x0 + (w - u)

  if (is.infinite(m)) {
    repeat {
      if (L <= lower) break
      if (logf(L) <= logy) break
      L <- L - w
    }
    repeat {
      if (R >= upper) break
      if (logf(R) <= logy) break
      R <- R + w
    }
  } else if (m > 1) {
    J <- floor(stats::runif(1, 0, m))
    K <- (m - 1) - J
    while (J > 0) {
      if (L <= lower) break
      if (logf(L) <= logy) break
      L <- L - w
      J <- J - 1L
    }
    while (K > 0) {
      if (R >= upper) break
      if (logf(R) <= logy) break
      R <- R + w
      K <- K - 1L
    }
  }

  L <- max(L, lower)
  R <- min(R, upper)

  repeat {
    x1 <- stats::runif(1, L, R)
    gx1 <- logf(x1)
    if (gx1 >= logy) {
      return(list(
        value = x1,
        log_density = gx1,
        evals = n_eval,
        interval = c(lower = L, upper = R)
      ))
    }
    if (x1 > x0) {
      R <- x1
    } else {
      L <- x1
    }
  }
}
#
CheckLossFn = function(p0,diff){diff*p0 - diff*as.numeric(diff<0)}
#
dlm_df = function(y, model, df, dim.df, s.priors = list(l0=1,S0=10), just.lik=FALSE){
  ### Gets the Time Series Length / Replicate number
  y = check_ts(y)
  TT = nrow(y)
  ### Gets the State Parameter dimension and Prior Distribution Parameters
  m0 = model$m0
  C0 = model$C0
  l0 = s.priors$l0
  S0 = s.priors$S0
  n = length(m0)
  ### Constructs F and G
  FF = model$FF
  GG = model$GG
  ### Variable Saving
  ### Posterior Distribution
  m = matrix(0,TT,n)
  C = array(0,c(TT,n,n))
  ### Predictive State Distribution
  a = matrix(0,TT,n)
  R = array(0,dim = c(TT,n,n))
  P = array(0,dim = c(TT,n,n))
  W = array(0,dim = c(TT,n,n))
  ### One-Step Ahead Forecast
  f = matrix(0,TT,1)
  Q = array(0,c(TT,1,1))
  inv.Q = array(0,c(TT,1,1))
  ### Regression Variables
  e = matrix(0,TT,1)
  A = array(0,c(TT,n,1))
  ### Sample Variance
  S = vector("numeric",TT)
  l = vector("numeric",TT)

  # Prior Dim Check
  m0 = matrix(m0,n,1)
  C0 = .exdqlm_regularize_cov(matrix(C0,n,n), context = "dlm_df_C0")
  ### Discount Factor Blocking
  df.mat = make_df_mat(df,dim.df,n)

  ### First Update
  ### One-step state forecast
  a[1,]  = GG[,,1] %*% m0
  P[1,,] = .exdqlm_regularize_cov(GG[,,1] %*% C0 %*% t(GG[,,1]), context = "dlm_df_P_t1")
  W[1,,] = df.mat * P[1,,]
  R[1,,] = .exdqlm_regularize_cov(P[1,,] + W[1,,], context = "dlm_df_R_t1")
  ### One-step ahead forecast
  f[1,] = t(FF[,1]) %*% a[1,]
  Q[1,,] = matrix(.exdqlm_regularize_var(1 + t(FF[,1]) %*% R[1,,] %*% FF[,1], context = "dlm_df_Q_t1"),1,1)
  inv.Q[1,,] = matrix(1 / Q[1,,], 1, 1)
  ### Auxilary Variables
  e[1,]  = as.matrix(y[1,] - f[1,],1,1)
  A[1,,] = R[1,,] %*% FF[,1] %*% inv.Q[1,,]
  ### Variance update
  l[1] = l0 + 1
  S[1] = l0 * S0 / l[1] + (t(e[1,]) %*% inv.Q[1,,] %*% e[1,] / l[1])
  ### Posterior Distribution
  m[1,]  = a[1,] + as.matrix(A[1,,],n,1) %*% e[1,]
  C[1,,] = .exdqlm_regularize_cov(
    R[1,,] - as.matrix(A[1,,],n,1) %*% Q[1,,] %*% t(A[1,,]),
    context = "dlm_df_C_t1"
  )

  for(i in 2:TT){
    ### One-step state forecast
    a[i,]  = GG[,,i] %*% m[i-1,]
    P[i,,] = .exdqlm_regularize_cov(GG[,,i] %*% C[i-1,,] %*% t(GG[,,i]), context = sprintf("dlm_df_P_t%d", i))
    W[i,,] = df.mat * P[i,,]
    R[i,,] = .exdqlm_regularize_cov(P[i,,] + W[i,,], context = sprintf("dlm_df_R_t%d", i))
    ### One-step ahead forecast
    f[i,] = t(FF[,i]) %*% a[i,]
    Q[i,,] = matrix(.exdqlm_regularize_var(1 + t(FF[,i])%*% R[i,,]%*% FF[,i], context = sprintf("dlm_df_Q_t%d", i)),1,1)
    inv.Q[i,,] = matrix(1 / Q[i,,], 1, 1)
    ### Auxilary Variables
    e[i,]  = as.matrix(y[i,] - f[i,],1,1)
    A[i,,] = as.matrix(R[i,,] %*% FF[,i] %*% inv.Q[i,,],n,1)
    ### Variance update
    l[i] = l[i-1] + 1
    S[i] = l[i-1] * S[i-1] / l[i] + (t(e[i,]) %*% inv.Q[i,,] %*% e[i,] / l[i])
    ### Posterior Distribution
    m[i,]  = a[i,] + as.matrix(A[i,,],n,1) %*% e[i,]
    C[i,,] = .exdqlm_regularize_cov(
      R[i,,] - as.matrix(A[i,,],n,1) %*% Q[i,,] %*% t(as.matrix(A[i,,],n,1)),
      context = sprintf("dlm_df_C_t%d", i)
    )
  }

  ### Adjust By Variance
  R[1,,] = S0 * R[1,,]
  Q[1,,]   = S0 * Q[1,,]
  C[1,,]   = S[1] * C[1,,]
  for(i in 2:TT){
    R[i,,] = S[i-1] * R[i,,]
    Q[i,,]   = S[i-1] * Q[i,,]
    C[i,,]   = S[i] * C[i,,]
  }

  # Calculate Log-Likelihood
  det.Q = log(abs(Q[1,,])) ; llik = lgamma((l0+1)/2)-lgamma(l0/2)-log(pi*l0)/2-det.Q/2-(l0+1)*log(1+t(e[1,])%*%inv.Q[1,,]%*%e[1,]/l0)/2
  for(t in 2:TT){
    det.Q = log(abs(Q[t,,]))
    llik = llik + lgamma((l[t-1]+1)/2)-lgamma(l[t-1]/2)-log(pi*l[t-1])/2-det.Q/2-(l[t-1]+1)*log(1+t(e[t,])%*%inv.Q[t,,]%*%e[t,]/l[t-1])/2
  }
  if(just.lik){
    return(list(llik = llik))
  }

  ## SMOOTHING
  ### Initializes recursive relations
  sa = matrix(0,TT,n)
  sR = array(0, dim = c(TT,n,n))
  ### Runs the recursive equations
  sa[TT,]  = m[TT,]
  sR[TT,,] = C[TT,,]
  for(k in 1:(TT-1)){
  ### Computes the Auxilary recursion Variable B
    t_next <- TT - k + 1L
    R_next_info <- .exdqlm_cov_inverse(R[t_next,,], context = sprintf("dlm_df_smoother_R_t%d", t_next))
    B = C[TT-k,,] %*% t(GG[,,t_next]) %*% R_next_info$inverse
    sa[TT-k,] = m[TT-k,] + B %*% (sa[TT-k+1,] - a[TT-k+1,])
    sR[TT-k,,] = .exdqlm_regularize_cov(
      C[TT-k,,] + B %*% (sR[TT-k+1,,] - R_next_info$Sigma) %*% t(B),
      context = sprintf("dlm_df_smoother_sR_t%d", TT - k)
    )
  }
  ### Adjusts the variance update
  for(k in 1:TT){
    sR[TT-k,,] = S[TT] * sR[TT-k,,] / S[TT-k]
  }
  return(list(fm = m, fC = C, m = sa, C = sR,model = model, s = S, n = l))
}
#
make_df_mat = function(df,dim.df,n){
  if(sum(dim.df)!=n){ stop("sum of component dimensions given in dim.df does not match m0") }
  if(length(df)!=length(dim.df)){ stop("length of component discount factors does not match length of component dimensions") }
  dfs = rep(df,dim.df)
  n.dfs = length(dim.df)
  ind.dfs = c(0,sapply(1:length(dim.df),function(x){sum(dim.df[1:x])}),n)
  df.mat = matrix(0,n,n)
  for(j in 1:n.dfs){
    df.mat[(ind.dfs[j]+1):ind.dfs[(j+1)],(ind.dfs[j]+1):ind.dfs[(j+1)]] = (1-dfs[ind.dfs[(j+1)]])/dfs[ind.dfs[(j+1)]]
  }
  return(df.mat)
}

.exdqlm_dynamic_cov_controls <- function() {
  eig_floor <- as.numeric(getOption("exdqlm.dynamic.cov_eig_floor", 1e-10))[1]
  if (!is.finite(eig_floor) || eig_floor <= 0) eig_floor <- 1e-10
  eig_cap <- as.numeric(getOption("exdqlm.dynamic.cov_eig_cap", 1e8))[1]
  if (!is.finite(eig_cap) || eig_cap <= eig_floor) eig_cap <- 1e8
  q_floor <- as.numeric(getOption("exdqlm.dynamic.q_floor", 1e-10))[1]
  if (!is.finite(q_floor) || q_floor <= 0) q_floor <- 1e-10
  q_cap <- as.numeric(getOption("exdqlm.dynamic.q_cap", 1e12))[1]
  if (!is.finite(q_cap) || q_cap <= q_floor) q_cap <- 1e12
  list(eig_floor = eig_floor, eig_cap = eig_cap, q_floor = q_floor, q_cap = q_cap)
}

.exdqlm_regularize_cov <- function(Sigma, context = "dynamic_covariance",
                                   eig_floor = NULL, eig_cap = NULL) {
  ctrl <- .exdqlm_dynamic_cov_controls()
  if (is.null(eig_floor)) eig_floor <- ctrl$eig_floor
  if (is.null(eig_cap)) eig_cap <- ctrl$eig_cap

  M <- as.matrix(Sigma)
  M <- (M + t(M)) / 2
  if (any(!is.finite(M))) {
    bad <- which(!is.finite(M), arr.ind = TRUE)[1, ]
    stop(sprintf("%s has non-finite entry at (%d,%d)", context, bad[1], bad[2]))
  }

  sv <- svd(M)
  vals <- pmin(pmax(sv$d, eig_floor), eig_cap)
  M_reg <- sv$u %*% diag(vals, nrow(M)) %*% t(sv$u)
  M_reg <- (M_reg + t(M_reg)) / 2
  attr(M_reg, "svd_u") <- sv$u
  attr(M_reg, "svd_d") <- vals
  M_reg
}

.exdqlm_cov_inverse <- function(Sigma, context = "dynamic_covariance",
                                eig_floor = NULL, eig_cap = NULL) {
  M_reg <- .exdqlm_regularize_cov(Sigma, context = context, eig_floor = eig_floor, eig_cap = eig_cap)
  u <- attr(M_reg, "svd_u")
  d <- attr(M_reg, "svd_d")
  inv <- u %*% diag(1 / d, nrow(M_reg)) %*% t(u)
  list(Sigma = M_reg, inverse = inv, values = d)
}

.exdqlm_regularize_var <- function(q, context = "dynamic_q",
                                   q_floor = NULL, q_cap = NULL) {
  ctrl <- .exdqlm_dynamic_cov_controls()
  if (is.null(q_floor)) q_floor <- ctrl$q_floor
  if (is.null(q_cap)) q_cap <- ctrl$q_cap
  qv <- as.numeric(q)[1]
  if (is.na(qv)) {
    stop(sprintf("%s is NA", context))
  }
  if (!is.finite(qv)) {
    qv <- sign(qv) * q_cap
  }
  qv <- min(max(qv, q_floor), q_cap)
  qv
}
#
check_mod = function(model){
  if(!is.exdqlm(model)){
    stop("Model must be an 'exdqlm' object. To create an 'exdqlm', use functions as.exdqlm(), seasMod(), or polytrendMod().")
  }
  
  ## check all dimensions
  # m0
  if(!is.vector(model$m0)){
    if(nrow(model$m0) != 1 & ncol(model$m0) != 1){
      stop("m0 must be a vector, or a matrix with 1 column or 1 row")
    }
  }
  model$m0 = as.matrix(c(model$m0))
  p = nrow(model$m0)
  # C0
  model$C0 = as.matrix(model$C0)
  if(p != dim(model$C0)[1] & p != dim(model$C0)[2]){
    stop("C0 must be a square matrix matching the dimension of m0")
    }
  if(!all.equal(model$C0, t(model$C0)) | !all(eigen(model$C0)$values >= 0)){
    stop("C0 must be a covariance matrix")
  }
  # FF
  if(!is.vector(model$FF)){
    if(nrow(model$FF) != p & ncol(model$FF) != p){
      stop("FF must be a vector of length matching the dimension of m0, or a matrix with number of rows matching the dimension of m0")
    }
    if(ncol(model$FF) == p){
      model$FF = t(model$FF)
    }
  }else{
    if(length(model$FF) != p){
      stop("FF must be a vector of length matching the dimension of m0, or a matrix with number of rows matching the dimension of m0")
    }else{
      model$FF = matrix(model$FF,p,1)
    }
  }
  # GG
  if(is.null(dim(model$GG)[3])){
    model$GG = as.matrix(model$GG)
  }else{
    if(is.na(dim(model$GG)[3])){
      model$GG = as.matrix(model$GG)
    }else{
      model$GG = as.array(model$GG)
    }
  }
  if(p != dim(model$GG)[1] & p != dim(model$GG)[2]){
    stop("GG must be a square matrix matching the dimension of m0, or an array with first two dimensions matching the dimension of m0")
  }
  
  return(model)
}
#
check_logics = function(gam.init,sig.init,fix.gamma,fix.sigma,dqlm.ind){
  retval <- NULL
  retval$gam.init = gam.init
  retval$fix.gamma = fix.gamma
  retval$dqlm.ind = dqlm.ind
  if(dqlm.ind){
    if(gam.init!=0 | !fix.gamma){
      retval$gam.init <- gam.init <- 0
      retval$fix.gamma <- fix.gamma <- TRUE
    }
  }else{
    if(gam.init==0 && fix.gamma==TRUE){
      retval$dqlm.ind = TRUE
    }
  }
  if(fix.gamma & is.na(gam.init)){ stop("when fix.gamma = TRUE, gam.init must be specified") }
  if(fix.sigma & is.na(sig.init)){ stop("when fix.sigma = TRUE, sig.init must be specified") }
  return(retval)
}

.resolve_static_al_alias <- function(dqlm.ind, al.ind, dqlm_supplied, where) {
  if (!is.logical(dqlm.ind) || length(dqlm.ind) != 1L || is.na(dqlm.ind)) {
    stop("dqlm.ind must be a single non-missing logical.")
  }
  if (!is.null(al.ind)) {
    if (!is.logical(al.ind) || length(al.ind) != 1L || is.na(al.ind)) {
      stop("al.ind must be NULL or a single non-missing logical.")
    }
    if (isTRUE(dqlm_supplied) && !identical(isTRUE(dqlm.ind), isTRUE(al.ind))) {
      stop(
        sprintf(
          "%s received conflicting inputs: dqlm.ind = %s and al.ind = %s. Use one flag or set both equal.",
          where,
          as.character(dqlm.ind),
          as.character(al.ind)
        )
      )
    }
    dqlm.ind <- isTRUE(al.ind)
  }
  dqlm.ind
}
#
check_ts = function(dat){
  if(is.null(dim(dat))){
    dim(dat) <- c(length(dat),1)
  }
  if(all(dim(dat)>1) || all(dim(dat)==1)){
    stop("data must be univariate time-series")
  }
  if(dim(dat)[1]<dim(dat)[2]){
    dat = t(dat)
  }
  return(invisible(dat))
}

# Internal helpers for reduced AL (DQLM) variational updates.
.dqlm_gig_moment <- function(k, chi, psi, r) {
  chi <- pmax(as.numeric(chi), 1e-14)
  psi <- pmax(as.numeric(psi), 1e-14)
  z <- sqrt(chi * psi)
  num <- besselK(z, nu = k + r, expon.scaled = TRUE)
  den <- besselK(z, nu = k, expon.scaled = TRUE)
  ratio <- num / den
  ratio[!is.finite(ratio)] <- 1
  (sqrt(chi / psi)^r) * ratio
}

.dqlm_gig_elog <- function(k, chi, psi) {
  chi <- pmax(as.numeric(chi), 1e-14)
  psi <- pmax(as.numeric(psi), 1e-14)
  z <- sqrt(chi * psi)
  eps <- 1e-6
  logK <- function(nu) {
    val <- besselK(z, nu = nu, expon.scaled = TRUE)
    log(pmax(val, 1e-300)) - z
  }
  dlogK <- (logK(k + eps) - logK(k - eps)) / (2 * eps)
  0.5 * (log(chi) - log(psi)) + dlogK
}

.dqlm_gig_entropy <- function(k, chi, psi, E_inv_v, E_v, E_log_v) {
  chi <- pmax(as.numeric(chi), 1e-14)
  psi <- pmax(as.numeric(psi), 1e-14)
  z <- sqrt(chi * psi)
  logK <- log(pmax(besselK(z, nu = k, expon.scaled = TRUE), 1e-300)) - z
  logc <- (k / 2) * (log(psi) - log(chi)) - log(2) - logK
  sum(-logc - (k - 1) * E_log_v + 0.5 * (chi * E_inv_v + psi * E_v))
}

.vb_joint_controls <- function(tol_state, has_gamma = TRUE) {
  tol_state <- as.numeric(tol_state)[1]
  if (!is.finite(tol_state) || tol_state <= 0) tol_state <- 1e-3

  tol_sigma <- as.numeric(getOption("exdqlm.tol_sigma", tol_state))[1]
  if (!is.finite(tol_sigma) || tol_sigma <= 0) tol_sigma <- tol_state

  tol_gamma <- as.numeric(getOption("exdqlm.tol_gamma", tol_state))[1]
  if (!is.finite(tol_gamma) || tol_gamma <= 0) tol_gamma <- tol_state

  tol_elbo <- as.numeric(getOption("exdqlm.tol_elbo", pmax(1e-5, tol_state / 10)))[1]
  if (!is.finite(tol_elbo) || tol_elbo <= 0) tol_elbo <- pmax(1e-5, tol_state / 10)

  min_iter <- suppressWarnings(as.integer(getOption("exdqlm.vb.min_iter", 10L))[1])
  if (!is.finite(min_iter) || min_iter < 1L) min_iter <- 10L

  patience <- suppressWarnings(as.integer(getOption("exdqlm.vb.patience", 3L))[1])
  if (!is.finite(patience) || patience < 1L) patience <- 3L

  allow_elbo_drop <- as.numeric(getOption("exdqlm.vb.allow_elbo_drop", tol_elbo))[1]
  if (!is.finite(allow_elbo_drop) || allow_elbo_drop < 0) allow_elbo_drop <- tol_elbo

  list(
    tol_state = tol_state,
    tol_sigma = tol_sigma,
    tol_gamma = if (isTRUE(has_gamma)) tol_gamma else NA_real_,
    tol_elbo = tol_elbo,
    min_iter = min_iter,
    patience = patience,
    allow_elbo_drop = allow_elbo_drop,
    has_gamma = isTRUE(has_gamma)
  )
}

.vb_joint_step <- function(
  iter,
  d_state,
  d_sigma,
  d_gamma = NA_real_,
  d_elbo = NA_real_,
  controls,
  compute_elbo = TRUE,
  stable_count = 0L
) {
  state_ok <- is.finite(d_state) && (d_state <= controls$tol_state)
  sigma_ok <- is.finite(d_sigma) && (d_sigma <= controls$tol_sigma)
  gamma_ok <- if (isTRUE(controls$has_gamma)) {
    is.finite(d_gamma) && (d_gamma <= controls$tol_gamma)
  } else {
    TRUE
  }

  if (!isTRUE(compute_elbo)) {
    elbo_ok <- TRUE
  } else if (is.finite(d_elbo)) {
    elbo_ok <- (d_elbo <= controls$tol_elbo) && (d_elbo >= -controls$allow_elbo_drop)
  } else {
    elbo_ok <- FALSE
  }

  joint_ok <- state_ok && sigma_ok && gamma_ok && elbo_ok
  stable_next <- if (iter >= controls$min_iter && joint_ok) stable_count + 1L else 0L
  stop_now <- stable_next >= controls$patience

  list(
    state_ok = state_ok,
    sigma_ok = sigma_ok,
    gamma_ok = gamma_ok,
    elbo_ok = elbo_ok,
    joint_ok = joint_ok,
    stable_count = stable_next,
    stop_now = stop_now
  )
}

.exdqlm_pos_truncnorm_moments <- function(mu, tau2) {
  mu <- as.numeric(mu)
  tau2 <- pmax(as.numeric(tau2), 1e-14)
  tau <- sqrt(tau2)
  alpha <- mu / tau
  log_Phi <- stats::pnorm(alpha, log.p = TRUE)
  log_phi <- stats::dnorm(alpha, log = TRUE)
  Lambda <- exp(log_phi - log_Phi)
  Phi <- exp(log_Phi)
  E_pos <- mu + tau * Lambda
  E2_pos <- tau2 + mu^2 + tau * mu * Lambda

  extreme_left <- is.finite(alpha) & alpha < -8
  if (any(extreme_left)) {
    x <- -alpha[extreme_left]
    mean_shift_std <- 1 / x - 2 / (x^3) + 10 / (x^5)
    Lambda[extreme_left] <- x + mean_shift_std
    E_pos[extreme_left] <- tau[extreme_left] * mean_shift_std
    E2_std <- 1 - x * mean_shift_std
    E2_pos[extreme_left] <- tau2[extreme_left] * pmax(E2_std, mean_shift_std^2)
  }

  E2_pos <- pmax(E2_pos, E_pos^2)
  list(
    mean = E_pos,
    second = E2_pos,
    sd = sqrt(pmax(E2_pos - E_pos^2, 0)),
    tau = tau,
    alpha = alpha,
    Phi = Phi,
    Lambda = Lambda
  )
}

.exdqlm_pos_truncnorm_entropy <- function(mu, tau2, moments = NULL, total = TRUE) {
  mu <- as.numeric(mu)
  tau2 <- pmax(as.numeric(tau2), 1e-14)
  tau <- sqrt(tau2)
  if (is.null(moments)) {
    moments <- .exdqlm_pos_truncnorm_moments(mu, tau2)
  }
  alpha <- mu / tau
  H <- 0.5 * log(2 * pi * tau2) + log(pmax(moments$Phi, 1e-12)) +
    0.5 * (1 - alpha * moments$Lambda)
  if (isTRUE(total)) sum(H) else H
}

.exdqlm_ldvb_sample_gaussian <- function(mu, Sigma, n_samp) {
  ns <- suppressWarnings(as.integer(n_samp)[1])
  if (!is.finite(ns) || ns < 1L) return(NULL)

  mu <- as.numeric(mu)
  p <- length(mu)
  S <- as.matrix(Sigma)
  if (!all(dim(S) == c(p, p))) {
    stop("Sigma must be a square covariance matrix matching length(mu).", call. = FALSE)
  }
  S <- (S + t(S)) / 2

  if (p == 1L) {
    sd1 <- sqrt(max(S[1, 1], 0))
    out <- matrix(mu[1] + stats::rnorm(ns, sd = sd1), ncol = 1L)
    colnames(out) <- names(mu)
    return(out)
  }

  L <- try(chol(S), silent = TRUE)
  if (inherits(L, "try-error")) {
    eig <- eigen(S, symmetric = TRUE)
    vals <- pmax(eig$values, .Machine$double.eps)
    L <- eig$vectors %*% diag(sqrt(vals), p) %*% t(eig$vectors)
  }

  Z <- matrix(stats::rnorm(ns * p), nrow = ns, ncol = p)
  out <- sweep(Z %*% L, 2L, mu, `+`)
  colnames(out) <- names(mu)
  out
}

.exal_static_ldvb_sample_posterior <- function(fit, n_samp) {
  ns <- suppressWarnings(as.integer(n_samp)[1])
  if (!is.finite(ns) || ns < 1L) return(NULL)

  beta_mu <- as.numeric(fit$qbeta$m)
  names(beta_mu) <- colnames(fit$X)
  beta_draws <- .exdqlm_ldvb_sample_gaussian(beta_mu, fit$qbeta$V, ns)
  if (is.null(colnames(beta_draws))) {
    colnames(beta_draws) <- colnames(fit$X)
  }

  out <- list(samp.beta = beta_draws)

  if (isTRUE(fit$dqlm.ind)) {
    a <- as.numeric(fit$qsig$a)[1]
    b <- as.numeric(fit$qsig$b)[1]
    out$samp.sigma <- 1 / stats::rgamma(ns, shape = a, rate = b)
    out$samp.gamma <- rep(0, ns)
    return(out)
  }

  ld_mu <- c(
    eta = as.numeric(fit$qsiggam$eta_hat),
    ell = as.numeric(fit$qsiggam$ell_hat)
  )
  ld_draws <- .exdqlm_ldvb_sample_gaussian(ld_mu, fit$qsiggam$Sigma, ns)
  bounds <- fit$misc$bounds
  L <- as.numeric(bounds["L"])
  U <- as.numeric(bounds["U"])
  if (!is.finite(L)) L <- as.numeric(bounds[1])
  if (!is.finite(U)) U <- as.numeric(bounds[2])
  out$samp.sigma <- exp(ld_draws[, 2L])
  out$samp.gamma <- L + (U - L) * stats::plogis(ld_draws[, 1L])
  out
}

.exdqlm_trace_summary <- function(x) {
  z <- as.numeric(x)
  z <- z[is.finite(z)]
  if (!length(z)) {
    return(list(
      mean = NA_real_,
      sd = NA_real_,
      q05 = NA_real_,
      median = NA_real_,
      q95 = NA_real_,
      min = NA_real_,
      max = NA_real_
    ))
  }
  qs <- stats::quantile(z, probs = c(0.05, 0.5, 0.95), na.rm = TRUE, names = FALSE, type = 8)
  list(
    mean = mean(z),
    sd = stats::sd(z),
    q05 = qs[1],
    median = qs[2],
    q95 = qs[3],
    min = min(z),
    max = max(z)
  )
}

.exdqlm_validate_crps_probs <- function(crps_probs) {
  crps_probs <- as.numeric(crps_probs)
  if (!length(crps_probs) ||
      any(!is.finite(crps_probs)) ||
      any(crps_probs <= 0 | crps_probs >= 1)) {
    stop("`crps_probs` must be a non-empty numeric vector with values strictly between 0 and 1.", call. = FALSE)
  }
  if (anyDuplicated(crps_probs)) {
    stop("`crps_probs` must not contain duplicate values.", call. = FALSE)
  }
  sort(crps_probs)
}

.exdqlm_validate_crps_weights <- function(crps_weights, n_probs) {
  if (is.null(crps_weights)) {
    return(rep(1 / n_probs, n_probs))
  }
  crps_weights <- as.numeric(crps_weights)
  if (length(crps_weights) != n_probs ||
      any(!is.finite(crps_weights)) ||
      any(crps_weights < 0) ||
      sum(crps_weights) <= 0) {
    stop("`crps_weights` must be a non-negative numeric vector with length equal to `crps_probs` and positive sum.", call. = FALSE)
  }
  crps_weights / sum(crps_weights)
}

.exdqlm_empirical_quantile_type7 <- function(z, probs) {
  z <- sort(as.numeric(z))
  z <- z[is.finite(z)]
  m <- length(z)
  if (!m) {
    return(rep(NA_real_, length(probs)))
  }
  h <- 1 + (m - 1) * probs
  lo <- floor(h)
  hi <- ceiling(h)
  z[lo] + (h - lo) * (z[hi] - z[lo])
}

.exdqlm_crps_sample_row <- function(y_true, draws_vec) {
  z <- sort(as.numeric(draws_vec))
  z <- z[is.finite(z)]
  m <- length(z)
  if (m < 2L || !is.finite(y_true)) {
    return(NA_real_)
  }
  mean(abs(z - y_true)) - sum((2 * seq_len(m) - m - 1) * z) / (m^2)
}

.exdqlm_crps_sample_vec <- function(y_true, draws_mat) {
  draws_mat <- as.matrix(draws_mat)
  stopifnot(length(y_true) == nrow(draws_mat))
  vapply(seq_len(nrow(draws_mat)), function(i) {
    .exdqlm_crps_sample_row(y_true[[i]], draws_mat[i, ])
  }, numeric(1))
}

.exdqlm_crps_row_validated <- function(y_true, draws_vec, probs, weights) {
  if (!is.finite(y_true)) {
    return(NA_real_)
  }
  qhat <- .exdqlm_empirical_quantile_type7(draws_vec, probs)
  if (!any(is.finite(qhat))) {
    return(NA_real_)
  }
  loss <- CheckLossFn(probs, y_true - qhat)
  2 * sum(weights * loss)
}

.exdqlm_crps_row <- function(y_true, draws_vec,
                             probs = seq(0.01, 0.99, by = 0.01),
                             weights = NULL) {
  probs <- .exdqlm_validate_crps_probs(probs)
  weights <- .exdqlm_validate_crps_weights(weights, length(probs))
  .exdqlm_crps_row_validated(y_true, draws_vec, probs, weights)
}

.exdqlm_crps_vec <- function(y_true, draws_mat,
                             probs = seq(0.01, 0.99, by = 0.01),
                             weights = NULL) {
  probs <- .exdqlm_validate_crps_probs(probs)
  weights <- .exdqlm_validate_crps_weights(weights, length(probs))
  draws_mat <- as.matrix(draws_mat)
  stopifnot(length(y_true) == nrow(draws_mat))
  vapply(seq_len(nrow(draws_mat)), function(i) {
    .exdqlm_crps_row_validated(y_true[[i]], draws_mat[i, ], probs, weights)
  }, numeric(1))
}

.exdqlm_normalize_vb_trace_vector <- function(x, n_iter, name, drop_init = FALSE, fill = NA_real_) {
  n_iter <- suppressWarnings(as.integer(n_iter)[1])
  if (!is.finite(n_iter) || n_iter < 0L) {
    stop("`n_iter` must be a non-negative integer.", call. = FALSE)
  }
  if (!n_iter) {
    return(rep(fill, 0L))
  }
  if (is.null(x) || !length(x)) {
    return(rep(fill, n_iter))
  }

  x <- as.numeric(x)
  if (isTRUE(drop_init) && length(x) == (n_iter + 1L)) {
    return(x[-1L])
  }
  if (length(x) == n_iter) {
    return(x)
  }
  if (length(x) < n_iter) {
    return(c(rep(fill, n_iter - length(x)), x))
  }

  expected <- if (isTRUE(drop_init)) {
    sprintf("%d or %d (including initialization)", n_iter, n_iter + 1L)
  } else {
    sprintf("at most %d", n_iter)
  }
  stop(
    sprintf("`%s` trace must have length %s; got %d.", name, expected, length(x)),
    call. = FALSE
  )
}

.exdqlm_make_vb_trace <- function(
  iter,
  engine = "VB",
  dqlm.ind = FALSE,
  elbo = NULL,
  sigma = NULL,
  gamma = NULL,
  delta_state = NULL,
  delta_sigma = NULL,
  delta_gamma = NULL,
  delta_s = NULL,
  delta_elbo = NULL,
  sigma_has_init = FALSE,
  gamma_has_init = FALSE
) {
  n_iter <- suppressWarnings(as.integer(iter)[1])
  if (!is.finite(n_iter) || n_iter < 0L) {
    stop("`iter` must be a non-negative integer.", call. = FALSE)
  }

  data.frame(
    iter = seq_len(n_iter),
    engine = rep(as.character(engine)[1], n_iter),
    dqlm.ind = rep(isTRUE(dqlm.ind), n_iter),
    elbo = .exdqlm_normalize_vb_trace_vector(elbo, n_iter, "elbo"),
    sigma = .exdqlm_normalize_vb_trace_vector(
      sigma, n_iter, "sigma", drop_init = sigma_has_init
    ),
    gamma = .exdqlm_normalize_vb_trace_vector(
      gamma, n_iter, "gamma", drop_init = gamma_has_init
    ),
    delta_state = .exdqlm_normalize_vb_trace_vector(delta_state, n_iter, "delta_state"),
    delta_sigma = .exdqlm_normalize_vb_trace_vector(delta_sigma, n_iter, "delta_sigma"),
    delta_gamma = .exdqlm_normalize_vb_trace_vector(delta_gamma, n_iter, "delta_gamma"),
    delta_s = .exdqlm_normalize_vb_trace_vector(delta_s, n_iter, "delta_s"),
    delta_elbo = .exdqlm_normalize_vb_trace_vector(delta_elbo, n_iter, "delta_elbo"),
    stringsAsFactors = FALSE
  )
}

.exdqlm_chain_health_metrics <- function(x, n_keep = length(x)) {
  z <- as.numeric(x)
  z <- z[is.finite(z)]
  n_keep <- suppressWarnings(as.numeric(n_keep)[1])
  if (!is.finite(n_keep) || n_keep <= 0) n_keep <- length(z)

  ess <- if (length(z) >= 10L) {
    tryCatch(as.numeric(coda::effectiveSize(coda::as.mcmc(z))), error = function(e) NA_real_)
  } else {
    NA_real_
  }

  acf1 <- if (length(z) >= 10L) {
    ac <- tryCatch(stats::acf(z, lag.max = 1L, plot = FALSE)$acf, error = function(e) NULL)
    if (is.null(ac) || length(ac) < 2L) NA_real_ else as.numeric(ac[2L])
  } else {
    NA_real_
  }

  geweke_absz <- if (length(z) >= 20L) {
    gz <- tryCatch(coda::geweke.diag(coda::as.mcmc(z))$z, error = function(e) NA_real_)
    as.numeric(abs(gz[1]))
  } else {
    NA_real_
  }

  half_drift <- if (length(z) >= 20L) {
    i <- floor(length(z) / 2L)
    s <- stats::sd(z)
    if (!is.finite(s) || s <= 0 || i < 5L || (length(z) - i) < 5L) {
      NA_real_
    } else {
      as.numeric(abs(mean(z[(i + 1L):length(z)]) - mean(z[seq_len(i)])) / s)
    }
  } else {
    NA_real_
  }

  list(
    n = as.integer(length(z)),
    ess = ess,
    ess_per1k = if (is.finite(ess) && is.finite(n_keep) && n_keep > 0) as.numeric(ess / n_keep * 1000) else NA_real_,
    acf1 = acf1,
    geweke_absz = geweke_absz,
    half_drift = half_drift
  )
}

# Reduced dynamic DQLM CAVI core (no gamma / no s_t block).
.run_dynamic_dqlm_cavi <- function(
  y, p0, model, df, dim.df,
  fix.sigma = FALSE, sig.init = NA_real_,
  tol = 0.1, n.samp = 200L,
  PriorSigma = NULL,
  verbose = TRUE,
  exps0 = NULL,
  max_iter = 200L,
  engine = "VB"
) {
  y <- y
  TT <- length(y)
  p <- length(model$m0)

  GG <- array(model$GG, c(p, p, TT))
  FF <- matrix(model$FF, p, TT)
  m0 <- as.numeric(model$m0)
  C0 <- as.matrix(model$C0)
  df.mat <- make_df_mat(df, dim.df, p)

  # Fixed AL constants at gamma = 0.
  A_tau <- (1 - 2 * p0) / (p0 * (1 - p0))
  B_tau <- 2 / (p0 * (1 - p0))

  if (is.null(PriorSigma)) {
    m_sigma <- 1
    v_sigma <- 10
    PriorSigma <- list(
      a_sig = (m_sigma^2) / v_sigma + 2,
      b_sig = (m_sigma^3) / v_sigma + m_sigma
    )
  }
  a0 <- as.numeric(PriorSigma$a_sig)[1]
  b0 <- as.numeric(PriorSigma$b_sig)[1]
  if (!is.finite(a0) || !is.finite(b0) || a0 <= 0 || b0 <= 0) {
    stop("PriorSigma must define positive finite a_sig and b_sig.")
  }

  sig0 <- if (!is.na(sig.init)) as.numeric(sig.init)[1] else 1
  if (!is.finite(sig0) || sig0 <= 0) sig0 <- 1
  if (isTRUE(fix.sigma) && is.na(sig.init)) {
    stop("fix.sigma=TRUE requires a finite sig.init in reduced DQLM CAVI.")
  }

  # Initialize q(v) moments and q(sigma) moments.
  E_v <- rep(sig0, TT)
  E_inv_v <- rep(1 / sig0, TT)
  E_log_v <- rep(0, TT)
  kappa <- 1 / sig0
  E_sigma <- sig0
  E_inv_sigma <- 1 / sig0
  E_log_sigma <- log(sig0)
  shape_sigma <- NA_real_
  scale_sigma <- NA_real_

  # Local R smoother used for both updates and dynamic ELBO state block.
  update_theta_reduced <- function(ex.f, ex.q) {
    m <- sm <- matrix(NA_real_, p, TT)
    C <- sC <- array(NA_real_, c(p, p, TT))
    a_store <- matrix(NA_real_, p, TT)
    P_store <- array(NA_real_, c(p, p, TT))
    f_vec <- rep(NA_real_, TT)
    q_vec <- rep(NA_real_, TT)
    sfe <- rep(NA_real_, TT)

    # Forward filter
    a <- as.vector(GG[, , 1] %*% m0)
    P <- .exdqlm_regularize_cov(GG[, , 1] %*% C0 %*% t(GG[, , 1]), context = "dqlm_cavi_P_t1")
    R <- .exdqlm_regularize_cov(P + df.mat * P, context = "dqlm_cavi_R_t1")
    f <- as.numeric(t(FF[, 1]) %*% a + ex.f[1])
    q <- .exdqlm_regularize_var(t(FF[, 1]) %*% R %*% FF[, 1] + ex.q[1], context = "dqlm_cavi_q_t1")
    m[, 1] <- a + as.vector(t(R) %*% FF[, 1]) * (y[1] - f) / q
    C[, , 1] <- .exdqlm_regularize_cov(
      R - (t(R) %*% FF[, 1] %*% t(FF[, 1]) %*% R) / q,
      context = "dqlm_cavi_C_t1"
    )
    a_store[, 1] <- a
    P_store[, , 1] <- R
    f_vec[1] <- f
    q_vec[1] <- q
    sfe[1] <- (y[1] - f) / sqrt(q)

    if (TT >= 2) {
      for (t in 2:TT) {
        a <- as.vector(GG[, , t] %*% m[, (t - 1)])
        P <- .exdqlm_regularize_cov(GG[, , t] %*% C[, , (t - 1)] %*% t(GG[, , t]), context = sprintf("dqlm_cavi_P_t%d", t))
        R <- .exdqlm_regularize_cov(P + df.mat * P, context = sprintf("dqlm_cavi_R_t%d", t))
        f <- as.numeric(t(FF[, t]) %*% a + ex.f[t])
        # Keep matrix shape (1 x p) so covariance update uses a p x p outer product.
        fB <- t(FF[, t]) %*% R
        q <- .exdqlm_regularize_var(fB %*% FF[, t] + ex.q[t], context = sprintf("dqlm_cavi_q_t%d", t))
        m[, t] <- a + as.vector(t(fB)) * (y[t] - f) / q
        C[, , t] <- .exdqlm_regularize_cov(
          R - (t(fB) %*% fB) / q,
          context = sprintf("dqlm_cavi_C_t%d", t)
        )
        a_store[, t] <- a
        P_store[, , t] <- R
        f_vec[t] <- f
        q_vec[t] <- q
        sfe[t] <- (y[t] - f) / sqrt(q)
      }
    }

    # Backward smoothing
    sC[, , TT] <- C[, , TT]
    sm[, TT] <- m[, TT]
    if (TT >= 2) {
      for (t in (TT - 1):1) {
        Pn_info <- .exdqlm_cov_inverse(P_store[, , (t + 1)], context = sprintf("dqlm_cavi_smoother_R_t%d", t + 1L))
        J <- C[, , t] %*% t(GG[, , (t + 1)]) %*% Pn_info$inverse
        sm[, t] <- m[, t] + J %*% (sm[, (t + 1)] - a_store[, (t + 1)])
        sC[, , t] <- .exdqlm_regularize_cov(
          C[, , t] + J %*% (sC[, , (t + 1)] - Pn_info$Sigma) %*% t(J),
          context = sprintf("dqlm_cavi_smoother_C_t%d", t)
        )
      }
    }

    exps <- apply(FF * sm, 2, sum)
    vars <- vapply(seq_len(TT), function(t) {
      as.numeric(t(FF[, t]) %*% sC[, , t] %*% FF[, t])
    }, numeric(1))
    exps2 <- exps^2 + vars

    # Dynamic ELBO state block via pseudo-model identity.
    y_star <- y - ex.f
    log_py_star <- -0.5 * sum(log(2 * pi * q_vec) + ((y - f_vec)^2) / q_vec)
    E_log_pseudo <- -0.5 * sum(log(2 * pi * ex.q) + (vars + (y_star - exps)^2) / ex.q)
    elbo_alpha <- as.numeric(log_py_star - E_log_pseudo)

    list(
      exps = exps,
      vars = vars,
      exps2 = exps2,
      standard.forecast.errors = sfe,
      sm = sm,
      sC = sC,
      fm = m,
      fC = C,
      elbo_alpha = elbo_alpha
    )
  }

  # Initial theta moments
  if (!is.null(exps0)) {
    if (length(exps0) != TT) stop("exps0 must have same length as y.")
    exps_init <- as.numeric(exps0)
  } else {
    init_dlm <- dlm_df(y, model, df, dim.df, s.priors = list(l0 = 1, S0 = sig0), just.lik = FALSE)
    exps_init <- apply(FF * t(init_dlm$m), 2, sum)
  }

  prev_exps <- exps_init
  prev_sigma <- E_sigma
  iter <- 0L
  stable_count <- 0L
  seq.sigma <- E_sigma
  elbo.seq <- numeric(0)
  delta_state <- numeric(0)
  delta_sigma <- numeric(0)
  delta_elbo <- numeric(0)
  controls <- .vb_joint_controls(tol_state = tol, has_gamma = FALSE)
  stop_reason <- "max_iter"

  .exdqlm_progress(
    "LDVB start",
    tol = tol,
    .verbose = verbose
  )
  tictoc::tic("run time")
  while (iter < max_iter) {
    iter <- iter + 1L

    # q(alpha_{0:T}) update using reduced pseudo-observation model
    ex.f <- A_tau / pmax(E_inv_v, 1e-12)
    ex.q <- B_tau / pmax(kappa * E_inv_v, 1e-12)
    theta.out <- update_theta_reduced(ex.f, ex.q)

    # Residual moments under q(alpha_t)
    E_r <- y - theta.out$exps
    E_r2 <- y^2 - 2 * y * theta.out$exps + theta.out$exps2

    # q(v_t) update (lambda = 1/2)
    chi <- (kappa / B_tau) * E_r2
    psi <- kappa * (2 + (A_tau^2) / B_tau)
    chi <- pmax(chi, 1e-12)
    psi <- pmax(psi, 1e-12)

    E_inv_v <- sqrt(psi / chi)
    E_v <- sqrt(chi / psi) * (1 + 1 / sqrt(chi * psi))
    E_log_v <- .dqlm_gig_elog(0.5, chi, psi)

    # q(sigma) update (or fixed sigma)
    if (!isTRUE(fix.sigma)) {
      shape_sigma <- a0 + 1.5 * TT
      scale_sigma <- b0 + sum(E_v) + (1 / (2 * B_tau)) * sum(E_inv_v * E_r2 - 2 * A_tau * E_r + (A_tau^2) * E_v)
      scale_sigma <- pmax(scale_sigma, 1e-12)

      E_inv_sigma <- shape_sigma / scale_sigma
      E_sigma <- if (shape_sigma > 1) scale_sigma / (shape_sigma - 1) else scale_sigma / shape_sigma
      E_log_sigma <- log(scale_sigma) - digamma(shape_sigma)
      kappa <- E_inv_sigma
    } else {
      E_sigma <- as.numeric(sig.init)[1]
      E_inv_sigma <- 1 / E_sigma
      E_log_sigma <- log(E_sigma)
      kappa <- E_inv_sigma
      shape_sigma <- NA_real_
      scale_sigma <- NA_real_
    }

    # Dynamic ELBO (reduced model)
    E_log_p_sigma <- a0 * log(b0) - lgamma(a0) - (a0 + 1) * E_log_sigma - b0 * E_inv_sigma
    E_log_p_v <- -TT * E_log_sigma - E_inv_sigma * sum(E_v)
    E_log_p_y <- - (TT / 2) * log(2 * pi) - (TT / 2) * log(B_tau) - (TT / 2) * E_log_sigma -
      0.5 * sum(E_log_v) -
      (E_inv_sigma / (2 * B_tau)) * sum(E_inv_v * E_r2 - 2 * A_tau * E_r + (A_tau^2) * E_v)
    H_sigma <- if (!isTRUE(fix.sigma)) {
      shape_sigma + log(scale_sigma) + lgamma(shape_sigma) - (shape_sigma + 1) * digamma(shape_sigma)
    } else {
      0
    }
    H_v <- .dqlm_gig_entropy(0.5, chi, psi, E_inv_v, E_v, E_log_v)
    elbo <- as.numeric(theta.out$elbo_alpha + E_log_p_sigma + E_log_p_v + E_log_p_y + H_sigma + H_v)
    elbo.seq <- c(elbo.seq, elbo)

    # Convergence diagnostics (joint: state + sigma + ELBO)
    d_state <- max(abs(theta.out$exps - prev_exps))
    d_sigma <- abs(E_sigma - prev_sigma)
    d_elbo <- if (length(elbo.seq) >= 2L) {
      elbo.seq[length(elbo.seq)] - elbo.seq[length(elbo.seq) - 1L]
    } else {
      NA_real_
    }
    step <- .vb_joint_step(
      iter = iter,
      d_state = d_state,
      d_sigma = d_sigma,
      d_elbo = d_elbo,
      controls = controls,
      compute_elbo = TRUE,
      stable_count = stable_count
    )
    stable_count <- step$stable_count
    delta_state <- c(delta_state, d_state)
    delta_sigma <- c(delta_sigma, d_sigma)
    delta_elbo <- c(delta_elbo, d_elbo)

    prev_exps <- theta.out$exps
    prev_sigma <- E_sigma
    seq.sigma <- c(seq.sigma, E_sigma)

    if (iter %% 5L == 0L) {
      .exdqlm_progress(
        "LDVB progress",
        model = "DQLM",
        iter = iter,
        elbo = elbo,
        .verbose = verbose
      )
    }

    if (step$stop_now) {
      stop_reason <- "joint_converged"
      break
    }
  }
  run.time <- tictoc::toc(quiet = TRUE)
  .exdqlm_progress(
    "LDVB done",
    model = "DQLM",
    status = if (identical(stop_reason, "joint_converged")) "converged" else "stopped",
    iter = iter,
    runtime_sec = run.time$toc - run.time$tic,
    .verbose = verbose
  )

  # Posterior sampling from variational factors
  ns <- as.integer(n.samp)
  if (!isTRUE(fix.sigma) && is.finite(shape_sigma) && is.finite(scale_sigma)) {
    samp.sigma <- 1 / stats::rgamma(ns, shape = shape_sigma, rate = scale_sigma)
  } else {
    samp.sigma <- rep(E_sigma, ns)
  }

  # q(v_t): require the package C++ Devroye sampler.
  samp.v <- t(.sample_gig_devroye_required(
    ns, p = 0.5, a = psi, b_vec = chi, context = "q(v_t) sampling"
  ))

  # q(theta_t): independent Gaussian draws per t from smoothed marginals.
  samp.theta <- array(NA_real_, dim = c(p, TT, ns))
  for (t in seq_len(TT)) {
    mu_t <- as.numeric(theta.out$sm[, t])
    S_t <- matrix(theta.out$sC[, , t], nrow = p, ncol = p)
    S_t <- (S_t + t(S_t)) / 2
    chol_t <- tryCatch(chol(S_t), error = function(e) NULL)
    if (is.null(chol_t)) {
      eig <- eigen(S_t, symmetric = TRUE)
      vals <- pmax(eig$values, 1e-12)
      chol_t <- eig$vectors %*% diag(sqrt(vals), p) %*% t(eig$vectors)
    }
    Z <- matrix(stats::rnorm(p * ns), nrow = p, ncol = ns)
    samp.theta[, t, ] <- mu_t + chol_t %*% Z
  }

  # Posterior predictive draws at gamma = 0 (AL model).
  samp.post.pred <- matrix(NA_real_, nrow = TT, ncol = ns)
  for (t in seq_len(TT)) {
    theta_t <- matrix(samp.theta[, t, ], nrow = p, ncol = ns)
    xb <- as.numeric(crossprod(FF[, t], theta_t))
    samp.post.pred[t, ] <- rexal(ns, p0, xb, samp.sigma, 0)
  }

  sig.out <- list(
    E.sigma = E_sigma,
    E.inv.sigma = E_inv_sigma,
    E.log.sigma = E_log_sigma,
    shape = shape_sigma,
    scale = scale_sigma
  )
  vts.out <- list(
    E.uts = E_v,
    E.inv.uts = E_inv_v,
    E.log.uts = sum(E_log_v),
    uts.chi = chi,
    uts.psi = psi,
    uts.lambda = 0.5
  )

  list(
    y = y,
    run.time = (run.time$toc - run.time$tic),
    iter = iter,
    dqlm.ind = TRUE,
    model = model,
    p0 = p0,
    df = df,
    dim.df = dim.df,
    sig.init = sig.init,
    seq.sigma = seq.sigma,
    samp.theta = samp.theta,
    samp.post.pred = samp.post.pred,
    map.standard.forecast.errors = theta.out$standard.forecast.errors,
    samp.sigma = samp.sigma,
    samp.vts = samp.v,
    theta.out = theta.out,
    sig.out = sig.out,
    vts.out = vts.out,
    fix.sigma = fix.sigma,
    fix.gamma = TRUE,
    diagnostics = list(
      elbo = elbo.seq,
      vb_trace = .exdqlm_make_vb_trace(
        iter = iter,
        engine = engine,
        dqlm.ind = TRUE,
        elbo = elbo.seq,
        sigma = seq.sigma,
        gamma = NULL,
        delta_state = delta_state,
        delta_sigma = delta_sigma,
        delta_gamma = NULL,
        delta_s = NULL,
        delta_elbo = delta_elbo,
        sigma_has_init = TRUE
      ),
      convergence = list(
        converged = identical(stop_reason, "joint_converged"),
        stop_reason = stop_reason,
        iter = iter,
        stable_count = stable_count,
        criteria = controls,
        final = list(
          delta_state = if (length(delta_state)) utils::tail(delta_state, 1L) else NA_real_,
          delta_sigma = if (length(delta_sigma)) utils::tail(delta_sigma, 1L) else NA_real_,
          delta_gamma = NA_real_,
          delta_elbo = if (length(delta_elbo)) utils::tail(delta_elbo, 1L) else NA_real_
        )
      ),
      deltas = list(
        state = delta_state,
        sigma = delta_sigma,
        gamma = rep(NA_real_, length(delta_state)),
        elbo = delta_elbo
      )
    )
  )
}


# Reduced static DQLM CAVI core (no gamma / no s block).
.run_static_dqlm_cavi <- function(
  y, X, p0,
  max_iter = 1000L,
  tol = 1e-4,
  b0 = NULL,
  V0 = NULL,
  beta_prior_obj = NULL,
  a_sigma = 1,
  b_sigma = 1,
  init = NULL,
  n.samp = 200L,
  verbose = TRUE
) {
  y <- as.numeric(y)
  X <- as.matrix(X)
  storage.mode(X) <- "double"
  n <- length(y)
  p <- ncol(X)
  if (nrow(X) != n) stop("nrow(X) must match length(y).")
  if (!(p0 > 0 && p0 < 1)) stop("p0 must be in (0,1).")

  if (is.null(b0)) b0 <- rep(0, p)
  if (is.null(V0)) V0 <- diag(1e6, p)
  V0 <- as.matrix(V0)
  if (!all(dim(V0) == c(p, p))) stop("V0 must be p x p.")
  if (is.null(beta_prior_obj)) {
    beta_prior_obj <- .static_beta_prior_make(
      beta_prior = "ridge",
      p = p,
      b0 = b0,
      V0 = V0,
      beta_prior_controls = NULL,
      warn_rhs_b0 = FALSE,
      warn_rhs_V0 = FALSE
    )
  }

  A_tau <- (1 - 2 * p0) / (p0 * (1 - p0))
  B_tau <- 2 / (p0 * (1 - p0))

  # Initialization
  m_beta <- if (!is.null(init$beta)) as.numeric(init$beta) else rep(0, p)
  if (length(m_beta) != p) m_beta <- rep(m_beta[1], p)
  V_beta <- V0
  beta_state <- beta_prior_obj$init_vb()

  sigma0 <- if (!is.null(init$sigma)) as.numeric(init$sigma)[1] else 1
  if (!is.finite(sigma0) || sigma0 <= 0) sigma0 <- 1
  a_q <- a_sigma + 1.5 * n
  b_q <- a_q * sigma0
  kappa <- a_q / b_q

  ell <- rep(1, n) # E[1 / v_t]
  nu <- rep(1, n)  # E[v_t]
  mlogv <- rep(0, n)
  chi <- rep(1, n)
  psi <- 1

  converged <- FALSE
  elbo_trace <- numeric(0)
  iter <- 0L
  stable_count <- 0L
  seq.sigma <- sigma0
  delta_beta <- numeric(0)
  delta_sigma <- numeric(0)
  delta_elbo <- numeric(0)
  controls <- .vb_joint_controls(tol_state = tol, has_gamma = FALSE)
  stop_reason <- "max_iter"

  .exdqlm_progress(
    "LDVB start",
    max_iter = as.integer(max_iter),
    tol = tol,
    .verbose = verbose
  )

  t0 <- proc.time()[3]
  for (iter in seq_len(as.integer(max_iter))) {
    prev_m_beta <- m_beta
    prev_sigma <- if (a_q > 1) b_q / (a_q - 1) else b_q / a_q

    # (1) q(beta): Normal
    W <- (kappa / B_tau) * ell
    Xw <- X * sqrt(W)
    prior_sys <- beta_prior_obj$beta_system_vb(beta_state)
    V_inv <- crossprod(Xw) + prior_sys$Prec
    Uc <- tryCatch(chol(V_inv), error = function(e) NULL)
    if (is.null(Uc)) Uc <- chol(V_inv + 1e-10 * diag(p))
    V_beta <- chol2inv(Uc)

    rhs <- prior_sys$h + (kappa / B_tau) * (
      crossprod(X, ell * y) - A_tau * colSums(X)
    )
    m_beta <- as.numeric(V_beta %*% rhs)
    beta_state <- beta_prior_obj$update_vb(
      beta_state,
      list(m = m_beta, V = V_beta)
    )

    # Residual moments under q(beta)
    m_eta <- as.numeric(X %*% m_beta)
    s_eta <- rowSums((X %*% V_beta) * X)
    E_r <- y - m_eta
    E_r2 <- s_eta + E_r^2

    # (2) q(v_t): GIG(lambda=1/2)
    chi <- (kappa / B_tau) * E_r2
    psi <- kappa * (2 + (A_tau^2) / B_tau)
    chi <- pmax(chi, 1e-12)
    psi <- pmax(psi, 1e-12)

    ell <- sqrt(psi / chi)
    nu <- sqrt(chi / psi) * (1 + 1 / sqrt(chi * psi))
    mlogv <- .dqlm_gig_elog(0.5, chi, psi)

    # (3) q(sigma): Inverse-gamma
    a_q <- a_sigma + 1.5 * n
    b_q <- b_sigma + sum(nu) + (1 / (2 * B_tau)) *
      sum(ell * E_r2 - 2 * A_tau * E_r + (A_tau^2) * nu)
    b_q <- pmax(b_q, 1e-12)
    kappa <- a_q / b_q

    E_sigma <- if (a_q > 1) b_q / (a_q - 1) else b_q / a_q
    E_inv_sigma <- kappa
    E_log_sigma <- log(b_q) - digamma(a_q)

    # ELBO
    E_log_p_beta <- beta_prior_obj$elbo_vb(
      beta_state,
      list(m = m_beta, V = V_beta)
    )

    E_log_p_sigma <- a_sigma * log(b_sigma) - lgamma(a_sigma) -
      (a_sigma + 1) * E_log_sigma - b_sigma * E_inv_sigma

    E_log_p_v <- -n * E_log_sigma - E_inv_sigma * sum(nu)

    E_log_p_y <- -(n / 2) * log(2 * pi) - (n / 2) * log(B_tau) -
      (n / 2) * E_log_sigma - 0.5 * sum(mlogv) -
      (E_inv_sigma / (2 * B_tau)) *
      sum(ell * E_r2 - 2 * A_tau * E_r + (A_tau^2) * nu)

    logdetVb <- as.numeric(determinant(V_beta, logarithm = TRUE)$modulus)
    H_beta <- 0.5 * (p * (1 + log(2 * pi)) + logdetVb)
    H_sigma <- a_q + log(b_q) + lgamma(a_q) - (a_q + 1) * digamma(a_q)
    H_v <- .dqlm_gig_entropy(0.5, chi, rep(psi, n), ell, nu, mlogv)

    elbo <- as.numeric(E_log_p_beta + E_log_p_sigma + E_log_p_v + E_log_p_y + H_beta + H_sigma + H_v)
    elbo_trace <- c(elbo_trace, elbo)

    d_beta <- max(abs(m_beta - prev_m_beta))
    d_sigma <- abs(E_sigma - prev_sigma)
    d_elbo <- if (length(elbo_trace) >= 2) {
      elbo_trace[length(elbo_trace)] - elbo_trace[length(elbo_trace) - 1]
    } else {
      NA_real_
    }
    step <- .vb_joint_step(
      iter = iter,
      d_state = d_beta,
      d_sigma = d_sigma,
      d_elbo = d_elbo,
      controls = controls,
      compute_elbo = TRUE,
      stable_count = stable_count
    )
    stable_count <- step$stable_count
    seq.sigma <- c(seq.sigma, E_sigma)
    delta_beta <- c(delta_beta, d_beta)
    delta_sigma <- c(delta_sigma, d_sigma)
    delta_elbo <- c(delta_elbo, d_elbo)

    if (iter %% 25L == 0L) {
      .exdqlm_progress(
        "LDVB progress",
        model = "AL special case",
        iter = iter,
        elbo = elbo,
        .verbose = verbose
      )
    }

    if (step$stop_now) {
      converged <- TRUE
      stop_reason <- "joint_converged"
      break
    }
  }

  t1 <- proc.time()[3]

  ret <- list(
    y = y,
    X = X,
    p0 = p0,
    dqlm.ind = TRUE,
    qbeta = list(m = m_beta, V = V_beta),
    qv = list(
      chi = chi,
      psi = psi,
      E_v = nu,
      E_inv_v = ell,
      E_log_v = mlogv
    ),
    qsig = list(
      a = a_q,
      b = b_q,
      E_sigma = if (a_q > 1) b_q / (a_q - 1) else b_q / a_q,
      E_inv_sigma = kappa,
      E_log_sigma = log(b_q) - digamma(a_q)
    ),
    converged = converged,
    iter = iter,
    run.time = as.numeric(t1 - t0),
    beta_prior = list(
      type = beta_prior_obj$type,
      controls = beta_prior_obj$controls,
      summary = beta_prior_obj$summary_vb(beta_state),
      state = if (.static_is_rhs_family(beta_prior_obj$type)) beta_state else NULL
    ),
    misc = list(
      p0 = p0,
      n = n,
      p = p,
      A = A_tau,
      B = B_tau,
      elbo = elbo_trace
    ),
    diagnostics = list(
      elbo = elbo_trace,
      vb_trace = .exdqlm_make_vb_trace(
        iter = iter,
        engine = "LDVB",
        dqlm.ind = TRUE,
        elbo = elbo_trace,
        sigma = seq.sigma,
        gamma = NULL,
        delta_state = delta_beta,
        delta_sigma = delta_sigma,
        delta_gamma = NULL,
        delta_s = NULL,
        delta_elbo = delta_elbo,
        sigma_has_init = TRUE
      ),
      convergence = list(
        converged = converged,
        stop_reason = stop_reason,
        iter = iter,
        stable_count = stable_count,
        criteria = controls,
        final = list(
          delta_state = if (length(delta_beta)) utils::tail(delta_beta, 1L) else NA_real_,
          delta_sigma = if (length(delta_sigma)) utils::tail(delta_sigma, 1L) else NA_real_,
          delta_gamma = NA_real_,
          delta_elbo = if (length(delta_elbo)) utils::tail(delta_elbo, 1L) else NA_real_
        )
      ),
      deltas = list(
        state = delta_beta,
        sigma = delta_sigma,
        gamma = rep(NA_real_, length(delta_beta)),
        elbo = delta_elbo
      )
    )
  )
  draws <- .exal_static_ldvb_sample_posterior(ret, n.samp)
  if (!is.null(draws)) {
    ret[names(draws)] <- draws
  }
  ret$run.time <- as.numeric(proc.time()[3] - t0)
  if (.static_is_rhs_family(beta_prior_obj$type)) {
    .static_rhs_maybe_warn_collapse(ret$beta_prior$summary, beta_prior_obj$controls)
  }
  .exdqlm_progress(
    "LDVB done",
    model = "AL special case",
    status = if (isTRUE(converged)) "converged" else "stopped",
    iter = iter,
    runtime_sec = ret$run.time,
    .verbose = verbose
  )
  ret
}