Fix reported inference, forecasting, serialization, time-index handling and numerical robustness

Package/DESCRIPTION

@@ -1,8 +1,8 @@
 Package: MSGARCH
 Type: Package
 Title: Markov-Switching GARCH Models
-Version: 2.51
-Date: 2022-12-05
+Version: 2.52
+Date: 2026-08-10
 Authors@R: c(person("David", "Ardia", role = c("aut"),
              email = "david.ardia.ch@gmail.com",
 			 comment = c(ORCID = "0000-0003-2823-782X")),

Package/NEWS

@@ -1,3 +1,29 @@
+Changes in Version 2.52
+  o The Hamilton filter now stabilises the regime log densities on their largest entry; anchoring on the smallest one overflowed to NaN when the regimes differed by more than about 1400 in log density
+  o The stationary distribution is obtained with a guarded solve: a transition matrix outside the simplex, reachable whenever fixed or regime-constant parameters are used with three or more regimes, no longer aborts the fit with an Armadillo exception
+  o Fitted-object methods no longer recycle observations when newdata is supplied with ts or zoo data, and forecast indexes advance by the series' own time step
+  o simulate() accepts nburn = 0
+  o FitML reports optimizer failures instead of failing on '$ operator is invalid for atomic vectors'
+  o UncVol averages the horizons after the burn-in rather than a descending range inside it
+  o The data must now contain at least two observations and no NA, NaN or infinite values; parameters are checked likewise
+  o Risk() validates alpha, nahead and ctr$nmesh, and warns when the evaluation grid does not cover enough of the predictive distribution to resolve the requested tail
+  o prior$sd is validated with the standard-deviation checker and must be finite and strictly positive
+  o Transition probabilities can no longer be passed to constraint.spec$fixed, which never worked
+  o FitMCMC drops fixed parameters from a user-supplied ctr$par0, as FitML already did
+  o The identification sort is skipped when constraint.spec$fixed is used, since relabelling regimes would move the fixed value into another regime
+  o CreateSpec validates switch.spec$K and rejects an explicitly heterogeneous distribution vector when expanding a single regime
+  o The mixture transition matrix built for Viterbi decoding is now row-stochastic (the decoded path is unchanged)
+  o Fixed the standard errors reported by summary(): the delta-method sandwich was transposed
+  o Pr(>|t|) in summary() is now the two-sided p-value, as its label states
+  o AIC/BIC now drop K-1 degrees of freedom per regime-constant parameter, not one
+  o Saved specifications and fits (saveRDS/readRDS) are usable again: the Rcpp modules are rebuilt
+  o predict() on an MCMC fit now averages over the posterior draws instead of returning the first draw
+  o Fixed a transposed cube write in the multi-regime in-sample CDF, which corrupted the first row of PIT(do.its = TRUE) and could overrun the buffer
+  o The is_log branches of the native mixture density and CDF now return the mixture instead of the last regime
+  o Two-sided p-values are evaluated on the lower tail, so they no longer cancel to zero for large statistics
+  o Removed the unreachable Sim methods for fit objects, which R CMD check flagged as unregistered
+  o Dropped the obsolete Rcpp:::LdFlags() call from src/Makevars
+  o Added regression tests for all of the above
 Changes in Version 2.51
   o Fix warning: use of bitwise '|' with boolean operands
 Changes in Version 2.5

Package/R/CondVol.R

@@ -29,7 +29,7 @@ f_CondVol <- function(object, par, data, do.its = FALSE, nahead = 1L, do.cumulat
   vol <- sqrt(vol)
   draw <- NULL
   if (!isTRUE(do.its)) {
-    tmp    <- mean(vol[dim(PredProb)[1]])
+    tmp    <- mean(vol[dim(PredProb)[1], ])
     vol    <- vector(mode = "numeric", length = nahead)
     vol[1] <- tmp
     if (nahead > 1) {
@@ -41,13 +41,7 @@ f_CondVol <- function(object, par, data, do.its = FALSE, nahead = 1L, do.cumulat
       vol[2:nahead] = apply(draw[2:nahead,, drop = FALSE], 1, sd)
     }
     names(vol) <- paste0("h=", 1:nahead)
-    if(zoo::is.zoo(data)){
-      vol = zoo::zooreg(vol, order.by = zoo::index(data)[length(data)]+(1:nahead))
-    }
-    if(is.ts(data)){
-      vol = zoo::zooreg(vol, order.by = zoo::index(data)[length(data)]+(1:nahead))
-      vol = as.ts(vol)
-    }
+    vol <- f_index_result(vol, data, nahead)
   } else {
     draw <- NULL
     if (nrow(par.check) > 1) {
@@ -58,13 +52,7 @@ f_CondVol <- function(object, par, data, do.its = FALSE, nahead = 1L, do.cumulat
       vol <- vol[1:length(data_)]
     }
     names(vol) <- paste0("t=", 1:(length(data_)))
-    if(zoo::is.zoo(data)){
-      vol = zoo::zooreg(vol, order.by = zoo::index(data))
-    }
-    if(is.ts(data)){
-      vol = zoo::zooreg(vol, order.by = zoo::index(data))
-      vol = as.ts(vol)
-    }
+    vol <- f_index_result(vol, data)
   }
   out = list()
   class(vol) <- c("MSGARCH_CONDVOL",class(vol))

Package/R/CreateSpec.R

@@ -200,6 +200,12 @@ CreateSpec <- function(variance.spec = list(model = c("sGARCH", "sGARCH")),
                        constraint.spec  = list(fixed = list(), regime.const = NULL),
                        prior = list(mean = list(), sd = list())) {
 
+  # whether the caller supplied these, as opposed to falling back on the
+  # two-regime defaults in the signature; only an explicit vector conflicts
+  # with expanding a single regime through switch.spec$K
+  bVarGiven  <- !missing(variance.spec)
+  bDistGiven <- !missing(distribution.spec)
+
   ## check
   variance.spec     <- f_check_variance_spec(variance.spec)
   distribution.spec <- f_check_distribution_spec(distribution.spec, length(variance.spec$model))
@@ -210,12 +216,18 @@ CreateSpec <- function(variance.spec = list(model = c("sGARCH", "sGARCH")),
   prior.sd          <- prior$sd
 
   if (!is.null(switch.spec$K)) {
-    if (length(variance.spec$model) > 1 | length(distribution.spec$model) > 1) {
+    if (length(switch.spec$K) != 1L || !is.numeric(switch.spec$K) ||
+        !is.finite(switch.spec$K) || switch.spec$K < 1 ||
+        switch.spec$K != round(switch.spec$K)) {
+      stop("switch.spec$K has to be a single positive whole number.")
+    }
+    if ((bVarGiven && length(variance.spec$model) > 1L) ||
+        (bDistGiven && length(distribution.spec$distribution) > 1L)) {
       stop("you can only use the variable K if you specified one
            regime in variance.spec and distribution.spec")
     } else {
-      variance.spec$model = rep(variance.spec$model, switch.spec$K)
-      distribution.spec$distribution = rep(distribution.spec$distribution, switch.spec$K)
+      variance.spec$model = rep(variance.spec$model[1L], switch.spec$K)
+      distribution.spec$distribution = rep(distribution.spec$distribution[1L], switch.spec$K)
     }
   }
 
@@ -309,7 +321,7 @@ CreateSpec <- function(variance.spec = list(model = c("sGARCH", "sGARCH")),
   }
   ## prior Sd
   if (length(prior.sd) >= 1) {
-    prior.sd <- f_check_parameterPriorMean(prior.sd, out$label)
+    prior.sd <- f_check_parameterPriorSd(prior.sd, out$label)
     out$prior.sd <- f_substitute_fixedpar(out$prior.sd, prior.sd)
     out$rcpp.func$set_sd(out$prior.sd)
   }

Package/R/FitMCMC.R

@@ -148,6 +148,10 @@ FitMCMC.MSGARCH_SPEC <- function(spec, data, ctr = list()) {
       par0 <- f_substitute_fixedpar(par0, spec$fixed.pars)
     }
     par0 <- f_unmapPar(par0, spec, do.plm = TRUE)
+    if (isTRUE(spec$fixed.pars.bool)) {
+      # as in FitML: the sampler works on the free parameters only
+      par0 <- f_remove_fixedpar(par0, spec$fixed.pars)
+    }
   }
   par    <- ctr$SamplerFUN(f_posterior = f_posterior, data = data_, spec = spec, par0 = par0, ctr = ctr)
   np     <- length(par0)
@@ -193,8 +197,15 @@ FitMCMC.MSGARCH_SPEC <- function(spec, data, ctr = list()) {
       par <- f_add_regimeconstpar_matrix(par, spec$K, spec$label)
     }
   }
-  if(isTRUE(ctr$do.sort)){
-    par <- f_sort_par(spec, par)
+  if (isTRUE(ctr$do.sort)) {
+    if (isTRUE(spec$fixed.pars.bool)) {
+      # sorting relabels the regimes by unconditional variance, which would move
+      # a parameter fixed in one regime into another; the constraint wins
+      message("do.sort is ignored: constraint.spec$fixed ties parameters to ",
+              "specific regimes, which the identification sort would relabel.")
+    } else {
+      par <- f_sort_par(spec, par)
+    }
   }
   par <- coda::mcmc(par)
   ctr$par0 <- par0

Package/R/FitML.R

@@ -128,12 +128,20 @@ FitML.MSGARCH_SPEC <- function(spec, data, ctr = list()) {
   }
   optimizer <- ctr$OptimFUN(vPw, f_nll, spec, data_, ctr$do.plm)
 
+  if (inherits(optimizer, "try-error")) {
+    stop("FitML: the optimizer failed with: ", as.character(optimizer))
+  }
+  if (is.null(optimizer$value) || is.null(optimizer$par)) {
+    stop("FitML: OptimFUN must return a list with elements 'value' and 'par'.")
+  }
+
   llk <- -optimizer$value
-
-  if (llk == 1e+10) {
-    str <- "FitML -> Error during optimization"
-    f_error(str)
-    stop()
+
+  # f_nll returns +1e10 when the likelihood cannot be evaluated, so a failed
+  # optimization comes back as llk = -1e10, not +1e10 as this test once assumed
+  if (!is.finite(llk) || llk <= -1e+10) {
+    stop("FitML: optimization failed; the log-likelihood could not be evaluated ",
+         "away from the starting values. Check the data and the specification.")
   }
 
   vPw <- optimizer$par
@@ -150,7 +158,11 @@ FitML.MSGARCH_SPEC <- function(spec, data, ctr = list()) {
   }
 
   par <- matrix(vPn, nrow = 1L, dimnames = list(NULL, names(vPn)))
-  par <- f_sort_par(spec, par)
+  if (!isTRUE(spec$fixed.pars.bool)) {
+    # see FitMCMC: the identification sort would relabel regimes and so break a
+    # parameter fixed in a particular one
+    par <- f_sort_par(spec, par)
+  }
   par <- as.vector(par)
   names(par) <- spec$label
   vPww <- f_unmapPar(par, spec, ctr$do.plm)

Package/R/Inference.R

@@ -35,12 +35,12 @@ f_InferenceFun <- function(vPw, data, spec, do.plm, mNegHessian = NULL) {
 
   mJacob      <- numDeriv::jacobian(f_mapPar, vPw_mod, spec = spec, do.plm = do.plm)
   mInvHessian <- MASS::ginv(mNegHessian)
-  mSandwitch  <- t(mJacob) %*% mInvHessian %*% mJacob
+  mSandwitch  <- mJacob %*% mInvHessian %*% t(mJacob)
 
   vSE   <- sqrt(diag(mSandwitch))
   vTest <- vPn/vSE
 
-  vPvalues <- 1 - pnorm(abs(vTest))
+  vPvalues <- 2 * pnorm(-abs(vTest))
 
   out[, "Estimate"]   <- vPn
   out[, "Std. Error"] <- vSE

Package/R/Likelihood.R

@@ -15,6 +15,12 @@ f_nll <- function(vPw, data, spec, do.plm) {
     vPn <- f_add_regimeconstpar(vPn, spec$K, spec$label)
   }
 
+  # the working-to-natural map can overflow for extreme trial values; treat that
+  # as an infeasible point rather than letting the strict parameter check throw
+  if (anyNA(vPn) || any(!is.finite(vPn))) {
+    return(1e+10)
+  }
+
   dLLK <- Kernel(spec, vPn, data, log = TRUE, do.prior = FALSE)
 
   if (!is.finite(dLLK)) {

Package/R/ParameterConstraints.R

@@ -1,14 +1,30 @@
 #################################################### fixed.pars ####
 
 f_check_parameterConstraints <- function(fixed.pars, vParNames) {
-  
+
   if (any(!names(fixed.pars) %in% vParNames)) {
     vWrongPars <- names(fixed.pars)[!names(fixed.pars) %in% vParNames]
-    stop(cat(paste("Wrong name in fixed.pars:", vWrongPars)))
+    stop("Wrong name in fixed.pars: ", paste(vWrongPars, collapse = ", "))
   }
-
+
+  # Transition probabilities are named in vParNames but cannot be fixed: the
+  # starting-value routine splits fixed parameters by regime and hands them to a
+  # single-regime specification that has no P_i_j, and the prior correction in
+  # Kernel() indexes prior.mean, which only covers the within-regime
+  # coefficients, so the log-posterior would silently become NA.
+  vIsP <- grepl("^P_", names(fixed.pars))
+  if (any(vIsP)) {
+    stop("Transition probabilities cannot be fixed through constraint.spec$fixed: ",
+         paste(names(fixed.pars)[vIsP], collapse = ", "))
+  }
+
+  vFixed <- unlist(fixed.pars)
+  if (length(vFixed) > 0L && (!is.numeric(vFixed) || any(!is.finite(vFixed)))) {
+    stop("Every entry of constraint.spec$fixed must be a finite number.")
+  }
+
   return(fixed.pars)
-  
+
 }
 
 f_remove_fixedpar <- function(vPar, fixed.pars) {

Package/R/Pit.R

@@ -132,14 +132,7 @@ PIT.MSGARCH_SPEC <- function(object, x = NULL, par = NULL, data = NULL,
     }
     tmp <- tmp/nrow(par)
     rownames(tmp) =  paste0("t=",1:length(data_))
-    if(zoo::is.zoo(data)){
-      tmp = zoo::zooreg(tmp, order.by = zoo::index(data))
-    }
-    if(is.ts(data)){
-      tmp = zoo::zooreg(tmp, order.by = zoo::index(data))
-      tmp = as.ts(tmp)
-      colnames(tmp) = rep("",ncol(tmp)) 
-    }
+    tmp <- f_index_result(tmp, data)
   } else {
     x <- matrix(x)
     if (ncol(x) != 1L) {
@@ -160,14 +153,7 @@ PIT.MSGARCH_SPEC <- function(object, x = NULL, par = NULL, data = NULL,
       }
     }
     rownames(tmp) <- paste0("h=",1:nahead)
-    if(zoo::is.zoo(data)){
-      tmp = zoo::zooreg(tmp, order.by = zoo::index(data)[length(data)]+(1:nahead))
-    }
-    if(is.ts(data)){
-      tmp = zoo::zooreg(tmp, order.by = zoo::index(data)[length(data)]+(1:nahead))
-      tmp = as.ts(tmp)
-      colnames(tmp) = rep("",ncol(tmp)) 
-    }
+    tmp <- f_index_result(tmp, data, nahead)
   }
   if (!isTRUE(ctr$do.return.draw)) {
     draw <- NULL
@@ -191,15 +177,7 @@ PIT.MSGARCH_SPEC <- function(object, x = NULL, par = NULL, data = NULL,
 #' @export
 PIT.MSGARCH_ML_FIT <- function(object, x = NULL, newdata = NULL,
                                do.norm = TRUE, do.its = FALSE, nahead = 1L, do.cumulative = FALSE, ctr = list(), ...) {
-  data = c(object$data, newdata)
-  if(is.ts(object$data)){
-    if(is.null(newdata)){
-      data = zoo::zooreg(data, order.by =  c(zoo::index(data)))
-    } else {
-      data = zoo::zooreg(data, order.by =  c(zoo::index(data),zoo::index(data)[length(data)]+(1:length(newdata))))
-    }
-    data = as.ts(data)
-  }
+  data <- f_combine_data(object$data, newdata)
   out <- PIT(object = object$spec, x = x, par = object$par, data = data,
              do.norm = do.norm, do.its = do.its, nahead = nahead, do.cumulative = do.cumulative,  ctr = ctr)
   return(out)
@@ -209,15 +187,7 @@ PIT.MSGARCH_ML_FIT <- function(object, x = NULL, newdata = NULL,
 #' @export
 PIT.MSGARCH_MCMC_FIT <- function(object, x = NULL, newdata = NULL,
                                  do.norm = TRUE, do.its = FALSE, nahead = 1L, do.cumulative = FALSE, ctr = list(), ...) {
-  data = c(object$data, newdata)
-  if(is.ts(object$data)){
-    if(is.null(newdata)){
-      data = zoo::zooreg(data, order.by =  c(zoo::index(data)))
-    } else {
-      data = zoo::zooreg(data, order.by =  c(zoo::index(data),zoo::index(data)[length(data)]+(1:length(newdata))))
-    }
-    data = as.ts(data)
-  }
+  data <- f_combine_data(object$data, newdata)
   out <- PIT(object = object$spec, x = x, par = object$par, data = data,
              do.norm = do.norm, do.its = do.its, nahead = nahead, do.cumulative = do.cumulative,  ctr = ctr)
   return(out)

Package/R/Posterior.R

@@ -14,6 +14,10 @@ f_posterior <- function(vPw, data, spec, PriorFun) {
     vPn <- f_add_regimeconstpar(vPn, spec$K, spec$label)
   }
 
+  if (anyNA(vPn) || any(!is.finite(vPn))) {
+    return(-1e10)
+  }
+
   dLLK <- Kernel(spec, vPn, data, log = TRUE, do.prior = TRUE) + sum(log(diag(abs(mJacob))))
 
   if (!is.finite(dLLK)) {

Package/R/PredPDF.R

@@ -120,14 +120,7 @@ PredPdf.MSGARCH_SPEC <- function(object, x = NULL, par = NULL, data = NULL,
     }
     tmp <- tmp/nrow(par)
     rownames(tmp) <-  paste0("t=",1:length(data_))
-    if(zoo::is.zoo(data)){
-      tmp = zoo::zooreg(tmp, order.by = zoo::index(data))
-    }
-    if(is.ts(data)){
-      tmp = zoo::zooreg(tmp, order.by = zoo::index(data))
-      tmp = as.ts(tmp)
-      colnames(tmp) = rep("",ncol(tmp)) 
-    }
+    tmp <- f_index_result(tmp, data)
   } else {
     if (is.null(x)) {
       stop("x is NULL: x must be a vector or a matrix of size N x 1")
@@ -151,14 +144,7 @@ PredPdf.MSGARCH_SPEC <- function(object, x = NULL, par = NULL, data = NULL,
       }
     }
     rownames(tmp) <- paste0("h=",1:nahead)
-    if(zoo::is.zoo(data)){
-      tmp = zoo::zooreg(tmp, order.by = zoo::index(data)[length(data)]+(1:nahead))
-    }
-    if(is.ts(data)){
-      tmp = zoo::zooreg(tmp, order.by = zoo::index(data)[length(data)]+(1:nahead))
-      tmp = as.ts(tmp)
-      colnames(tmp) = rep("",ncol(tmp)) 
-    }
+    tmp <- f_index_result(tmp, data, nahead)
   }
 
   if (!isTRUE(ctr$do.return.draw)) {
@@ -183,15 +169,7 @@ PredPdf.MSGARCH_SPEC <- function(object, x = NULL, par = NULL, data = NULL,
 #' @export
 PredPdf.MSGARCH_ML_FIT <- function(object, x = NULL, newdata = NULL,
                                 log = FALSE, do.its = FALSE, nahead = 1L, do.cumulative = FALSE, ctr = list(), ...) {
-  data <- c(object$data, newdata)
-  if(is.ts(object$data)){
-    if(is.null(newdata)){
-      data = zoo::zooreg(data, order.by =  c(zoo::index(data)))
-    } else {
-      data = zoo::zooreg(data, order.by =  c(zoo::index(data),zoo::index(data)[length(data)]+(1:length(newdata))))
-    }
-    data = as.ts(data)
-  }
+  data <- f_combine_data(object$data, newdata)
   out  <- PredPdf(object = object$spec, x = x, par = object$par, data = data,
                log = log, do.its = do.its, nahead = nahead, do.cumulative = do.cumulative, ctr = ctr)
   return(out)
@@ -201,15 +179,7 @@ PredPdf.MSGARCH_ML_FIT <- function(object, x = NULL, newdata = NULL,
 #' @export
 PredPdf.MSGARCH_MCMC_FIT <- function(object, x = NULL, newdata = NULL,
                                   log = FALSE, do.its = FALSE, nahead = 1L, do.cumulative = FALSE, ctr = list(), ...) {
-  data <- c(object$data, newdata)
-  if(is.ts(object$data)){
-    if(is.null(newdata)){
-      data = zoo::zooreg(data, order.by =  c(zoo::index(data)))
-    } else {
-      data = zoo::zooreg(data, order.by =  c(zoo::index(data),zoo::index(data)[length(data)]+(1:length(newdata))))
-    }
-    data = as.ts(data)
-  }
+  data <- f_combine_data(object$data, newdata)
   out  <- PredPdf(object = object$spec, x = x, par = object$par, data = data,
                log = log, do.its = do.its, nahead = nahead, do.cumulative = do.cumulative, ctr = ctr)
   return(out)