From 0d3fad6906ce7ea05b9a9f76063ba8a6d67c3b64 Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Mon, 10 Aug 2026 20:42:39 -0400
Subject: [PATCH 01/20] Fix the standard errors and p-values reported by
 summary()

f_InferenceFun built the delta-method sandwich as t(J) V J. numDeriv::jacobian
returns d f_i / d x_j, so the variance of the natural parameters is J V t(J).
The two agree only when J is symmetric, which it is not: the working-to-natural
map is triangular inside each regime (the sGARCH bound on beta is 0.9999 -
alpha1; the gjrGARCH and tGARCH bounds on beta also involve alpha2 and the
shape/skew parameters) and the transition-probability block is anti-diagonal.

On the MS(2)-GARCH(1,1)-Normal fit to SMI, six of the eight reported standard
errors were wrong, by factors from 0.06 to 16, and the two transition
probabilities had their standard errors exchanged. The corrected values agree
with the observed information computed directly in the natural parameterisation
to within 2%, where the old ones were off by up to 85%.

The Pr(>|t|) column reported 1 - pnorm(abs(t)), i.e. half a two-sided p-value
under a two-sided label.
---
 Package/R/Inference.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Package/R/Inference.R b/Package/R/Inference.R
index 1196a0c..1bea314 100644
--- a/Package/R/Inference.R
+++ b/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 * (1 - pnorm(abs(vTest)))
   
   out[, "Estimate"]   <- vPn
   out[, "Std. Error"] <- vSE

From 34d3f31c7dce69f9dd336406d5174f0fd9dc635d Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Mon, 10 Aug 2026 20:43:08 -0400
Subject: [PATCH 02/20] Rebuild the Rcpp modules when a saved spec or fit is
 reloaded

A MSGARCH_SPEC holds Rcpp module objects. R serializes external pointers as
NULL, so a spec or fit written with saveRDS comes back with dead pointers and
every method on it fails with "NULL value passed as symbol address". That is
what f_check_spec is for, but its recovery branch called
spec$rcpp.func$get_mean() and get_sd() -- the very pointer whose failure had
just triggered the branch -- so the rebuild could never run.

The two values it read there are already kept on the R side, in
spec$prior.mean and spec$prior.sd, which the next two lines were using anyway;
the C++ round trip was dead code. Dropping it makes the rebuild work: a
reloaded spec, ML fit or MCMC fit now returns exactly the same volatilities,
state probabilities, forecasts and information criteria as before it was
saved, with user-supplied priors preserved.

This affects the ordinary workflow of fitting a model, saving it, and
analysing it in a later session, and of sending a spec to a parallel worker.
---
 Package/R/Utils.R | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/Package/R/Utils.R b/Package/R/Utils.R
index 8ffd8d8..ea1d1cd 100644
--- a/Package/R/Utils.R
+++ b/Package/R/Utils.R
@@ -345,15 +345,18 @@ f_check_spec <- function(spec) {
     FALSE
   })
   if (!isTRUE(is.OK)) {
+    # The Rcpp module objects behind spec$rcpp.func do not survive
+    # serialization (saveRDS / parallel workers), so rebuild them. The prior
+    # mean and sd must be read from the R-side copies kept in the spec: the
+    # C++ getters go through the very pointer that is already dead.
     spec.new = f_spec(models = spec$name, do.mix = spec$is.mix)
-    prior.mean = spec$rcpp.func$get_mean()
-    prior.sd = spec$rcpp.func$get_sd()
-    names(prior.mean) = names(prior.sd) = spec$label[1:length(prior.mean)]
-    prior.mean[names(spec$prior.mean)] = spec$prior.mean
-    prior.sd[names(spec$prior.sd)] = spec$prior.sd
     spec$rcpp.func = spec.new$rcpp.func
-    spec$rcpp.func$set_mean(spec$prior.mean)
-    spec$rcpp.func$set_sd(spec$prior.sd)
+    if (!is.null(spec$prior.mean)) {
+      spec$rcpp.func$set_mean(spec$prior.mean)
+    }
+    if (!is.null(spec$prior.sd)) {
+      spec$rcpp.func$set_sd(spec$prior.sd)
+    }
   }
   return(spec)
 }

From 3a706565af6ca7dee5a30811eccc7de1892a89d4 Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Mon, 10 Aug 2026 20:43:08 -0400
Subject: [PATCH 03/20] Count regime-constant parameters correctly in AIC and
 BIC

Setting constraint.spec$regime.const ties a parameter across all K regimes, so
it leaves one free value where there were K and removes K - 1 degrees of
freedom (see f_rename_par, which strips name_2 ... name_K). dofMSGARCH
subtracted one per constrained name, which is right only at K = 2.

With a regime-constant shape parameter the reported df was 17 instead of 16 at
K = 3 and 27 instead of 25 at K = 4, so stats::AIC and stats::BIC over-
penalised the constrained model in every K >= 3 model-selection table.
---
 Package/R/Utils.R | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/Package/R/Utils.R b/Package/R/Utils.R
index ea1d1cd..fd7bc92 100644
--- a/Package/R/Utils.R
+++ b/Package/R/Utils.R
@@ -421,5 +421,8 @@ f_check_parameterPriorSd <- function(prior.sd, vParNames) {
 
 
 dofMSGARCH = function(object){
-  return(length(object$spec$par0) - length(object$spec[["regime.const.pars"]]) - length(object$spec[["fixed.pars"]]))
+  # each regime-constant parameter removes K - 1 free parameters, not one
+  return(length(object$spec$par0)
+         - length(object$spec[["regime.const.pars"]]) * (object$spec$K - 1L)
+         - length(object$spec[["fixed.pars"]]))
 }

From 742674eaa333364ecc567ceec7d5ae05dd86932f Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Mon, 10 Aug 2026 20:43:32 -0400
Subject: [PATCH 04/20] Average over the posterior when forecasting from an
 MCMC fit

In f_CondVol, vol is a (T + 1) x ndraw matrix and the one-step-ahead value was
taken as vol[dim(PredProb)[1]]. A single index into a matrix is linear
indexing, so that is the last row of the first column: predict() on a
MSGARCH_MCMC_FIT reported the forecast of MCMC draw #1 rather than the
posterior mean, and mean() of the resulting scalar was a no-op.

On a 100-draw chain fitted to SMI the reported value was 1.019753, the first
draw, against a posterior mean of 1.039776 and a range across draws of
0.999837 to 1.097337. Volatility() already averaged across draws correctly, so
the two methods disagreed on the same fit. Adding the missing comma leaves the
single-parameter (ML) path untouched.
---
 Package/R/CondVol.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Package/R/CondVol.R b/Package/R/CondVol.R
index fa0318d..ccdd60e 100644
--- a/Package/R/CondVol.R
+++ b/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) {

From 3a562fd795a9457df0c6288f4513f410dd78813d Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Mon, 10 Aug 2026 20:43:32 -0400
Subject: [PATCH 05/20] Add regression tests for the five fixes

test_Inference.R (new) pins the delta method against an independently computed
observed information. It is anchored on a single-regime GARCH(1,1)-Normal
because every parameter there is interior, so a central-difference Hessian of
the natural-scale negative log-likelihood is well conditioned: it agrees with
J V t(J) to 1e-4 in relative terms while the transposed sandwich is off by 44%
and 85%. It also checks the exact sandwich on a two-regime fit, where the
transition-probability block is the part that gets exchanged, that Pr(>|t|) is
two-sided, and that dofMSGARCH and the AIC/BIC arithmetic drop K - 1 values per
regime-constant parameter for K = 2, 3 and 4.

test_Serialization.R (new) round-trips a spec, an ML fit and an MCMC fit
through saveRDS/readRDS. R restores external pointers as NULL, so this
reproduces the cross-session failure inside a single session, and the tests
require Volatility, State, predict, AIC, summary, DIC and the log-kernel to
return values identical to the originals, with user priors preserved.

test_Volatility.R gains one block that recomputes the one-step-ahead
volatility draw by draw through the public interface -- each such call carries
a single parameter vector, so it cannot depend on how the draws are pooled --
and requires predict() on the MCMC fit to equal their mean.

Every block opens with a guard asserting its own precondition (the two
sandwich orientations really differ; the round trip really did invalidate the
pointers; the posterior mean really differs from the first draw) so that none
of them can pass vacuously if the surrounding code changes. All eleven blocks
fail on the unpatched package.
---
 Package/tests/testthat/test_Inference.R     | 167 ++++++++++++++++++++
 Package/tests/testthat/test_Serialization.R |  97 ++++++++++++
 Package/tests/testthat/test_Volatility.R    |  30 +++-
 3 files changed, 291 insertions(+), 3 deletions(-)
 create mode 100644 Package/tests/testthat/test_Inference.R
 create mode 100644 Package/tests/testthat/test_Serialization.R

diff --git a/Package/tests/testthat/test_Inference.R b/Package/tests/testthat/test_Inference.R
new file mode 100644
index 0000000..98624c6
--- /dev/null
+++ b/Package/tests/testthat/test_Inference.R
@@ -0,0 +1,167 @@
+testthat::context("Test Inference (standard errors, p-values, degrees of freedom)")
+
+data("SMI", package = "MSGARCH")
+
+# Single-regime GARCH(1,1)-Normal is the cleanest probe for the orientation of the
+# delta method: the working -> natural map is triangular (the upper bound on beta is
+# 0.9999 - alpha1, so d beta / d alpha1_tilde != 0 while d alpha1 / d beta_tilde == 0)
+# and no parameter sits near a bound, so the observed information in the natural
+# parameterisation is well conditioned and can be used as an independent reference.
+spec.sr <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                               distribution.spec = list(distribution = c("norm")),
+                               switch.spec = list(do.mix = FALSE, K = 1))
+fit.sr <- MSGARCH::FitML(spec.sr, data = SMI)
+
+f_jacob <- function(fit) {
+  vPw <- MSGARCH:::f_unmapPar(fit$par, fit$spec, fit$ctr$do.plm)
+  numDeriv::jacobian(MSGARCH:::f_mapPar, vPw, spec = fit$spec, do.plm = fit$ctr$do.plm)
+}
+
+testthat::test_that("Standard errors use J V J', not its transpose", {
+
+  J  <- f_jacob(fit.sr)
+  V  <- MASS::ginv(fit.sr$Inference$Hessian)
+  se <- unname(fit.sr$Inference$MatCoef[, "Std. Error"])
+
+  se.delta      <- sqrt(diag(J %*% V %*% t(J)))
+  se.transposed <- sqrt(diag(t(J) %*% V %*% J))
+
+  # guard against a vacuous test: the two orientations must actually differ here
+  testthat::expect_true(max(abs(se.delta - se.transposed)) > 1e-4)
+
+  testthat::expect_true(max(abs(se - se.delta)) < 1e-8)
+
+})
+
+testthat::test_that("Standard errors match the observed information in the natural scale", {
+
+  f_nll_natural <- function(vPn) {
+    -MSGARCH:::Kernel(fit.sr$spec, vPn, SMI, log = TRUE, do.prior = FALSE)
+  }
+
+  # central-difference Hessian with a small step, to stay away from the constraint
+  # boundaries at which the kernel is floored to -1e10
+  vPn  <- fit.sr$par
+  d    <- length(vPn)
+  step <- pmax(abs(vPn), 1) * 1e-4
+  mH   <- matrix(data = 0, nrow = d, ncol = d)
+  for (i in 1:d) {
+    for (j in i:d) {
+      ei <- rep(0, d); ei[i] <- step[i]
+      ej <- rep(0, d); ej[j] <- step[j]
+      mH[i, j] <- mH[j, i] <- (f_nll_natural(vPn + ei + ej) - f_nll_natural(vPn + ei - ej)
+                               - f_nll_natural(vPn - ei + ej) + f_nll_natural(vPn - ei - ej)) /
+                              (4 * step[i] * step[j])
+    }
+  }
+  se.natural <- sqrt(diag(solve(mH)))
+  se         <- unname(fit.sr$Inference$MatCoef[, "Std. Error"])
+
+  # agreement is ~1e-4 in relative terms; the transposed sandwich is off by 44% and 85%
+  testthat::expect_true(max(abs(se / se.natural - 1)) < 0.02)
+
+})
+
+testthat::test_that("Transition-probability standard errors are not swapped across regimes", {
+
+  # The P-block of the Jacobian is anti-diagonal (f_mapGamma enumerates the
+  # off-diagonal entries column-major but carries the row-major parameter names),
+  # so transposing the sandwich exchanges the two transition probabilities' errors.
+  spec.ms <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                                 distribution.spec = list(distribution = c("norm")),
+                                 switch.spec = list(do.mix = FALSE, K = 2))
+  par <- c(0.021631876185, 0.087024443479, 0.881493722371, 0.020659831566,
+           0.005396009353, 0.994040728662, 0.978348086740, 0.998703301894)
+  names(par) <- spec.ms$label
+
+  vPw <- MSGARCH:::f_unmapPar(par, spec.ms, FALSE)
+  inf <- MSGARCH:::f_InferenceFun(vPw, SMI, spec.ms, do.plm = FALSE)
+
+  J  <- numDeriv::jacobian(MSGARCH:::f_mapPar, vPw, spec = spec.ms, do.plm = FALSE)
+  V  <- MASS::ginv(inf$Hessian)
+  se <- unname(inf$MatCoef[, "Std. Error"])
+
+  se.delta      <- sqrt(diag(J %*% V %*% t(J)))
+  se.transposed <- sqrt(diag(t(J) %*% V %*% J))
+  iP            <- match(c("P_1_1", "P_2_1"), spec.ms$label)
+
+  # the two orientations must disagree on the P block, otherwise this proves nothing
+  testthat::expect_true(max(abs(se.delta[iP] - se.transposed[iP])) > 1e-4)
+
+  testthat::expect_true(max(abs(se - se.delta)) < 1e-8)
+  testthat::expect_true(all(is.finite(se)) && all(se > 0))
+
+})
+
+testthat::test_that("Pr(>|t|) is a two-sided p-value", {
+
+  mCoef <- fit.sr$Inference$MatCoef
+  testthat::expect_true(max(abs(mCoef[, "Pr(>|t|)"] -
+                                2 * (1 - stats::pnorm(abs(mCoef[, "t value"]))))) < 1e-12)
+  testthat::expect_true(all(mCoef[, "Pr(>|t|)"] >= 0 & mCoef[, "Pr(>|t|)"] <= 1))
+
+})
+
+testthat::test_that("Degrees of freedom drop K - 1 values per regime-constant parameter", {
+
+  # a regime-constant parameter leaves one free value where there were K, so it
+  # removes K - 1 degrees of freedom, not one. K = 2 cannot tell the two apart.
+  for (K in 2:4) {
+    spec.rc <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                                   distribution.spec = list(distribution = c("std")),
+                                   switch.spec = list(do.mix = FALSE, K = K),
+                                   constraint.spec = list(regime.const = c("nu")))
+    exp.dof <- length(spec.rc$label) - (K - 1L)
+    testthat::expect_equal(MSGARCH:::dofMSGARCH(list(spec = spec.rc)), exp.dof)
+  }
+
+  # unconstrained and fixed-parameter specifications must be unaffected
+  spec.free <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                                   distribution.spec = list(distribution = c("norm")),
+                                   switch.spec = list(do.mix = FALSE, K = 2))
+  testthat::expect_equal(MSGARCH:::dofMSGARCH(list(spec = spec.free)),
+                         length(spec.free$label))
+
+  spec.fix <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                                  distribution.spec = list(distribution = c("norm")),
+                                  switch.spec = list(do.mix = FALSE, K = 2),
+                                  constraint.spec = list(fixed = list(beta_1 = 0.8)))
+  testthat::expect_equal(MSGARCH:::dofMSGARCH(list(spec = spec.fix)),
+                         length(spec.fix$label) - 1L)
+
+})
+
+testthat::test_that("AIC and BIC use the free-parameter count of a constrained fit", {
+
+  # end-to-end on a fitted K = 2 model: logLik()'s df attribute and the AIC/BIC
+  # arithmetic must both follow dofMSGARCH
+  spec.rc <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                                 distribution.spec = list(distribution = c("std")),
+                                 switch.spec = list(do.mix = FALSE, K = 2),
+                                 constraint.spec = list(regime.const = c("nu")))
+  fit.rc  <- MSGARCH::FitML(spec.rc, data = SMI, ctr = list(do.se = FALSE))
+  exp.dof <- length(spec.rc$label) - 1L
+
+  testthat::expect_equal(as.integer(attr(stats::logLik(fit.rc), "df")), exp.dof)
+  testthat::expect_true(abs(AIC(fit.rc) - (-2 * fit.rc$loglik + 2 * exp.dof)) < 1e-8)
+  testthat::expect_true(abs(BIC(fit.rc) -
+                            (-2 * fit.rc$loglik + log(length(SMI)) * exp.dof)) < 1e-8)
+
+  # K = 2 cannot separate "one" from "K - 1", so repeat for K = 3 and K = 4 on a
+  # fit-shaped object: logLik.MSGARCH_ML_FIT only reads $loglik, $data and $spec
+  for (K in 3:4) {
+    spec.K <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                                  distribution.spec = list(distribution = c("std")),
+                                  switch.spec = list(do.mix = FALSE, K = K),
+                                  constraint.spec = list(regime.const = c("nu")))
+    fit.K  <- structure(list(loglik = -3000, data = SMI, spec = spec.K),
+                        class = "MSGARCH_ML_FIT")
+    exp.K  <- length(spec.K$label) - (K - 1L)
+
+    testthat::expect_equal(as.integer(attr(stats::logLik(fit.K), "df")), exp.K)
+    testthat::expect_true(abs(AIC(fit.K) - (-2 * fit.K$loglik + 2 * exp.K)) < 1e-8)
+    testthat::expect_true(abs(BIC(fit.K) -
+                              (-2 * fit.K$loglik + log(length(SMI)) * exp.K)) < 1e-8)
+  }
+
+})
diff --git a/Package/tests/testthat/test_Serialization.R b/Package/tests/testthat/test_Serialization.R
new file mode 100644
index 0000000..70cb5ea
--- /dev/null
+++ b/Package/tests/testthat/test_Serialization.R
@@ -0,0 +1,97 @@
+testthat::context("Test Serialization (saveRDS / readRDS round trip)")
+
+# A MSGARCH_SPEC carries Rcpp module objects. R serializes external pointers as NULL,
+# so any spec or fit written with saveRDS comes back with dead pointers and has to be
+# rebuilt by f_check_spec. That happens in one session too, which is what these tests
+# exploit; it is the same failure a user hits when reloading an overnight fit or when
+# shipping a spec to a parallel worker.
+
+data("SMI", package = "MSGARCH")
+
+f_roundtrip <- function(object) {
+  sFile <- tempfile(fileext = ".rds")
+  on.exit(unlink(sFile))
+  saveRDS(object, file = sFile)
+  return(readRDS(sFile))
+}
+
+spec <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                            distribution.spec = list(distribution = c("norm")),
+                            switch.spec = list(do.mix = FALSE, K = 2))
+par <- c(0.021631876185, 0.087024443479, 0.881493722371, 0.020659831566,
+         0.005396009353, 0.994040728662, 0.978348086740, 0.998703301894)
+
+testthat::test_that("The round trip really does invalidate the Rcpp pointers", {
+
+  # if this ever stops holding the tests below become vacuous
+  spec.rt <- f_roundtrip(spec)
+  testthat::expect_error(spec.rt$rcpp.func$get_sd())
+
+})
+
+testthat::test_that("A reloaded spec is usable and gives identical results", {
+
+  spec.rt <- f_roundtrip(spec)
+
+  exp.vol <- Volatility(object = spec,    par = par, data = SMI)
+  est.vol <- Volatility(object = spec.rt, par = par, data = SMI)
+  testthat::expect_true(max(abs(as.numeric(est.vol) - as.numeric(exp.vol))) < 1e-12)
+
+  exp.llk <- MSGARCH:::Kernel(spec,    par, SMI, log = TRUE, do.prior = FALSE)
+  est.llk <- MSGARCH:::Kernel(spec.rt, par, SMI, log = TRUE, do.prior = FALSE)
+  testthat::expect_true(abs(est.llk - exp.llk) < 1e-12)
+
+  exp.state <- State(object = spec,    par = par, data = SMI)$SmoothProb
+  est.state <- State(object = spec.rt, par = par, data = SMI)$SmoothProb
+  testthat::expect_true(max(abs(est.state - exp.state)) < 1e-12)
+
+})
+
+testthat::test_that("A reloaded spec keeps its user-supplied priors", {
+
+  spec.prior <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                                    distribution.spec = list(distribution = c("norm")),
+                                    switch.spec = list(do.mix = FALSE, K = 2),
+                                    prior = list(mean = list(beta_1 = 0.7),
+                                                 sd   = list(beta_1 = 0.1)))
+  exp.mean <- spec.prior$rcpp.func$get_mean()
+  exp.sd   <- spec.prior$rcpp.func$get_sd()
+
+  spec.rt <- MSGARCH:::f_check_spec(f_roundtrip(spec.prior))
+
+  testthat::expect_true(max(abs(spec.rt$rcpp.func$get_mean() - exp.mean)) < 1e-12)
+  testthat::expect_true(max(abs(spec.rt$rcpp.func$get_sd()   - exp.sd))   < 1e-12)
+
+})
+
+testthat::test_that("A reloaded MSGARCH_ML_FIT is usable", {
+
+  fit    <- MSGARCH::FitML(spec, data = SMI[1:500])
+  fit.rt <- f_roundtrip(fit)
+
+  testthat::expect_true(max(abs(as.numeric(Volatility(fit.rt)) -
+                                as.numeric(Volatility(fit)))) < 1e-12)
+
+  set.seed(1234)
+  exp.pred <- predict(object = fit,    nahead = 1L)$vol
+  set.seed(1234)
+  est.pred <- predict(object = fit.rt, nahead = 1L)$vol
+  testthat::expect_true(abs(as.numeric(est.pred) - as.numeric(exp.pred)) < 1e-12)
+
+  testthat::expect_true(abs(AIC(fit.rt) - AIC(fit)) < 1e-12)
+  testthat::expect_silent(invisible(capture.output(summary(fit.rt))))
+
+})
+
+testthat::test_that("A reloaded MSGARCH_MCMC_FIT is usable", {
+
+  set.seed(1234)
+  fit    <- MSGARCH::FitMCMC(spec, data = SMI[1:500],
+                             ctr = list(nburn = 100L, nmcmc = 100L, nthin = 1L))
+  fit.rt <- f_roundtrip(fit)
+
+  testthat::expect_true(max(abs(as.numeric(Volatility(fit.rt)) -
+                                as.numeric(Volatility(fit)))) < 1e-12)
+  testthat::expect_true(abs(DIC(fit.rt)$DIC - DIC(fit)$DIC) < 1e-12)
+
+})
diff --git a/Package/tests/testthat/test_Volatility.R b/Package/tests/testthat/test_Volatility.R
index 37f407b..e3f2b0f 100644
--- a/Package/tests/testthat/test_Volatility.R
+++ b/Package/tests/testthat/test_Volatility.R
@@ -18,12 +18,36 @@ testthat::test_that("Forecast", {
   
 })
 
+testthat::test_that("Forecast from an MCMC fit averages over the posterior draws", {
+
+  y <- SMI[1:500]
+  set.seed(1234)
+  fit  <- MSGARCH::FitMCMC(spec, data = y,
+                           ctr = list(nburn = 100L, nmcmc = 200L, nthin = 20L))
+  mPar <- as.matrix(fit$par)
+
+  # one-step-ahead volatility draw by draw, through the public interface; each of
+  # these calls carries a single parameter vector, so it cannot depend on how the
+  # draws are pooled afterwards
+  vVol <- vapply(seq_len(nrow(mPar)), function(i) {
+    as.numeric(predict(object = spec, par = mPar[i, ], newdata = y, nahead = 1L)$vol)
+  }, FUN.VALUE = numeric(1))
+
+  # guard: the posterior mean has to be distinguishable from the first draw,
+  # otherwise this test proves nothing
+  testthat::expect_true(abs(mean(vVol) - vVol[1]) > 1e-6)
+
+  est.forecast <- as.numeric(predict(object = fit, nahead = 1L)$vol)
+  testthat::expect_true(abs(est.forecast - mean(vVol)) < 1e-10)
+
+})
+
 testthat::test_that("Conditional Vol", {
-  
+
   tol <- 0.05
   est.Vol <- Volatility(object = spec, par = par, data = SMI)[2000]
   exp.Vol <- c(2.1321725800180471)
-  
+
   testthat::expect_true(max(abs(est.Vol - exp.Vol)) < tol)
-  
+
 })
\ No newline at end of file

From c93e21825ea82d640fe00304d414ba1198c3c4f7 Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Mon, 10 Aug 2026 20:43:56 -0400
Subject: [PATCH 06/20] Bump to 2.52 and record the fixes in NEWS

Also refreshes the Date field, which R CMD check --as-cran flags as stale
together with the unchanged version number. Drop this commit if the release
number should be decided separately; nothing else in the branch depends on it.
---
 Package/DESCRIPTION | 4 ++--
 Package/NEWS        | 7 +++++++
 2 files changed, 9 insertions(+), 2 deletions(-)

diff --git a/Package/DESCRIPTION b/Package/DESCRIPTION
index 5311c78..0fa3035 100644
--- a/Package/DESCRIPTION
+++ b/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")),
diff --git a/Package/NEWS b/Package/NEWS
index 97d5e9d..e04b454 100644
--- a/Package/NEWS
+++ b/Package/NEWS
@@ -1,3 +1,10 @@
+Changes in Version 2.52
+  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 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

From 6ce57bd76a7ebb23f65a8e02f60781061f1315c5 Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Mon, 10 Aug 2026 20:51:17 -0400
Subject: [PATCH 07/20] pr description

---
 PR_BODY.md      | 124 ++++++++++
 REVIEW.md       | 591 ++++++++++++++++++++++++++++++++++++++++++++++++
 REVIEW_codex.md | 148 ++++++++++++
 3 files changed, 863 insertions(+)
 create mode 100644 PR_BODY.md
 create mode 100644 REVIEW.md
 create mode 100644 REVIEW_codex.md

diff --git a/PR_BODY.md b/PR_BODY.md
new file mode 100644
index 0000000..d36a2e5
--- /dev/null
+++ b/PR_BODY.md
@@ -0,0 +1,124 @@
+Five bugs in the reporting and persistence layer around the likelihood, plus regression
+tests for each. The likelihood itself is not touched: no estimate, log-likelihood,
+conditional variance or state probability changes anywhere in this branch. What changes is
+what `summary()` prints, what `AIC`/`BIC` count, what `predict()` returns for a Bayesian
+fit, and whether a saved fit can be reloaded at all.
+
+Everything below is reproduced on `data("SMI")` with the shipped code. `R CMD check
+--as-cran` is unchanged by the branch (2 WARNINGs, 1 NOTE, `testthat` OK — and the version
+bump in the last commit clears one of the two WARNINGs).
+
+## 1. Standard errors used the delta method transposed — `R/Inference.R:38`
+
+```r
+mSandwitch <- t(mJacob) %*% mInvHessian %*% mJacob   # -> mJacob %*% mInvHessian %*% t(mJacob)
+```
+
+`numDeriv::jacobian` returns `∂f_i/∂x_j`, so `Var(g(θ̂)) = J V J'`. The two orientations
+agree only if `J` is symmetric, and it is not: the working→natural map is triangular inside
+each regime (the sGARCH bound on `beta` is `0.9999 − alpha1`; the gjrGARCH and tGARCH bounds
+on `beta` also involve `alpha2` and the shape/skew parameters), and the transition-probability
+block is *anti*-diagonal.
+
+On the default MS(2)-GARCH(1,1)-Normal fit to `SMI`, checked against the observed information
+computed directly in the natural parameterisation:
+
+| | natural-scale `H` | fixed (`J V J'`) | before (`J' V J`) |
+|---|---|---|---|
+| `alpha1_1` | 0.01606 | 0.01511 | **0.03413** |
+| `beta_1`   | 0.02197 | 0.02091 | **0.00958** |
+| `alpha1_2` | 0.00431 | 0.00437 | **0.00610** |
+| `beta_2`   | 0.00416 | 0.00426 | **0.00049** |
+| `P_1_1`    | 0.00982 | 0.00973 | **0.00059** |
+| `P_2_1`    | (at bound) | 0.03072 | **0.50250** |
+
+Six of eight were wrong, by factors from 0.06× to 16×. Because the P block is anti-diagonal,
+transposing **exchanged the two transition probabilities' standard errors**. Only `alpha0_k`,
+whose map is a plain `exp`, was unaffected.
+
+## 2. `Pr(>|t|)` was one-sided — `R/Inference.R:43`
+
+`1 - pnorm(abs(t))` under a two-sided label; now `2 * (1 - pnorm(abs(t)))`.
+
+## 3. `AIC`/`BIC` mis-counted `regime.const.pars` — `R/Utils.R:420`
+
+A regime-constant parameter leaves one free value where there were `K`, so it removes `K − 1`
+degrees of freedom (cf. `f_rename_par`, which strips `name_2 … name_K`). `dofMSGARCH`
+subtracted one, which is right only at `K = 2`: with a regime-constant shape parameter the df
+was 17 instead of 16 at `K = 3` and 27 instead of 25 at `K = 4`, always over-penalising the
+constrained model.
+
+## 4. Saved specs and fits could not be reloaded — `R/Utils.R:349-353`
+
+R serializes external pointers as `NULL`, so a `saveRDS`-ed spec or fit comes back with dead
+pointers. `f_check_spec` exists to rebuild them, but the rebuild branch called
+`spec$rcpp.func$get_mean()` / `get_sd()` — the very pointer whose failure had just triggered
+the branch:
+
+```
+Volatility(fit) : Error in .External(...): NULL value passed as symbol address
+State(fit)      : Error in .External(...): NULL value passed as symbol address
+predict(fit)    : Error in .External(...): NULL value passed as symbol address
+```
+
+The two values read there are already held on the R side in `spec$prior.mean` /
+`spec$prior.sd`, which the next two lines were using anyway, so the C++ round trip was dead
+code. Removing it makes the rebuild work; a reloaded spec, ML fit or MCMC fit now returns
+values identical to before it was saved, with user priors preserved.
+
+This is the ordinary workflow of fitting a model, saving it, and analysing it in a later
+session — or shipping a spec to a `parLapply` worker.
+
+## 5. `predict()` on an MCMC fit returned draw #1 — `R/CondVol.R:32`
+
+`vol` is `(T+1) × ndraw` and the one-step-ahead value was `vol[dim(PredProb)[1]]` — a single
+index into a matrix is linear indexing, i.e. the last row of the *first column*. On a
+100-draw chain fitted to `SMI` the reported value was 1.019753 (the first draw) against a
+posterior mean of 1.039776, with a spread of 0.999837–1.097337 across draws. `Volatility()`
+already averaged correctly, so the two methods disagreed on the same fit. The
+single-parameter (ML) path is unchanged.
+
+## Tests
+
+`test_Inference.R` and `test_Serialization.R` are new; `test_Volatility.R` gains one block.
+11 blocks, 34 assertions, +13s of check time. All 11 fail on the current code.
+
+Two things worth pointing out, since they are what makes the tests worth having:
+
+- The standard-error test is anchored on a **single-regime** GARCH(1,1)-Normal, not on the
+  MS(2) default. Every parameter there is interior, so a central-difference Hessian of the
+  natural-scale negative log-likelihood is well conditioned: it agrees with `J V J'` to 1e-4
+  in relative terms while the transposed sandwich is off by 44% and 85%. That is an
+  independent check on the *value*, not a restatement of the formula.
+- Every block opens with a guard asserting its own precondition — the two sandwich
+  orientations really differ for this model; the round trip really did invalidate the
+  pointers; the posterior mean really differs from the first draw — so none of them can pass
+  vacuously if the surrounding code changes.
+
+The `K = 3` and `K = 4` legs of the degrees-of-freedom test run against a fit-shaped list
+rather than a real fit, because a constrained `K ≥ 3` model cannot currently be fitted at
+all (see below).
+
+## Not in this branch
+
+The last commit (version bump + `NEWS`) is separable — drop it if the release number should
+be decided elsewhere; nothing depends on it.
+
+A review of the package turned up nine further issues that are **not** addressed here,
+several of them more serious than some of the above. The two worth flagging now:
+
+- **`src/MSgarch.h:283`** computes the ergodic distribution with a raw Armadillo `.i()` on
+  `I − P + U`, on every likelihood evaluation. When `do.plm = TRUE` — forced by `fixed.pars`
+  and `regime.const.pars`, and hard-coded in `FitMCMC` — the free transition entries are
+  mapped into `(0,1)` independently, so for `K ≥ 3` a row can leave the simplex, the matrix
+  can be exactly singular, and the uncaught exception aborts the entire run. Reproducible:
+  `FitML` on a `K = 3` spec with `regime.const = "nu"` dies at
+  `P = [[1, 0.999973, −0.999973], [0,1,0], [0,0,1]]`. This is why the tests above cannot fit
+  a constrained `K ≥ 3` model.
+- **`R/FitML.R:133`** guards optimisation failure with `if (llk == 1e+10)`, but `f_nll`
+  returns `+1e10` so `llk` is `−1e10`; and `f_OptimFUNDefault` wraps `optim` in `try()`, so
+  `optimizer$value` errors first. Every failure mode — including a single `NA` in the data,
+  which `f_check_y` lets through — therefore surfaces as
+  `Error in optimizer$value : $ operator is invalid for atomic vectors`.
+
+Happy to open these as separate issues or as a follow-up PR, whichever you prefer.
diff --git a/REVIEW.md b/REVIEW.md
new file mode 100644
index 0000000..deef0b7
--- /dev/null
+++ b/REVIEW.md
@@ -0,0 +1,591 @@
+# Review of the `MSGARCH` package (v2.51)
+
+Reviewed 2026-08-10 on R 4.5.2 (aarch64-apple-darwin20), Apple clang 21.
+Source: `package MSGARCH/Package/`, git `17017fa` (2022-12-05, maintainer K. Bluteau,
+David = aut). Fifth in the R-package sweep after `AdMit`, `bayesGARCH`, `DEoptim`, `GAS`.
+
+**Method.** (i) `R CMD check --as-cran` on a freshly built tarball; (ii) line-by-line read of
+all 31 `.R` files and all 30 `src/` files; (iii) analytic re-derivation of every conditional
+density, CDF, quantile function and truncated moment in `Normal.h` / `Student.h` / `Ged.h` /
+`Symmetric.h` / `Skewed.h`, and of the stationarity conditions in `sGARCH.h` / `gjrGARCH.h` /
+`tGARCH.h` / `eGARCH.h` / `sARCH.h`; (iv) numerical reproduction of every defect below on
+`data("SMI")` with the shipped binary. Everything reported here is reproduced, not inferred.
+
+> **Status (2026-08-10).** Items 1–4 of the work list at the bottom — **A1, A2, A3, A4 and
+> B1** — are **fixed** in `Package/` (7 changed lines across `R/Inference.R`, `R/Utils.R`,
+> `R/CondVol.R`). Each fix is verified below in the relevant section, `R CMD check --as-cran`
+> is byte-for-byte unchanged (2 WARNINGs, 1 NOTE; `testthat` OK), and all eight specifications
+> in a regression sweep (SR / MS-2 / MS-3 / mixture / heterogeneous / skew-t / GED / both
+> constraint paths) fit with finite, correct standard errors. Regression cover was added for
+> all five — `tests/testthat/test_Inference.R` (new), `tests/testthat/test_Serialization.R`
+> (new) and one block appended to `tests/testthat/test_Volatility.R`: 11 blocks / 34
+> assertions, all passing on the patched build and all 11 failing on stock 2.51.
+> `DESCRIPTION`'s version and `NEWS` are untouched — that is a release decision for the
+> maintainer. Everything from **B2** onwards is still open.
+
+> **Independent cross-check (Codex).** The whole of this document — the 14 claims, the five
+> fixes and the tests — was re-audited read-only by `codex-cli 0.139.0` working from source
+> alone (`codex exec --sandbox read-only`); its verbatim output is in `REVIEW_codex.md`.
+> It returned **CONFIRM on 13 of 14 claims, PARTIALLY CONFIRM on B3**, found no error in the
+> parts this review calls correct, and marked all five fixes OK. Its three pushbacks:
+>
+> 1. **B3's "trailing NA name" mechanism is overstated** — Codex expected `f_rename_par` to
+>    error on a too-long vector, or the fixed dimension to simply stay in the sampler.
+>    **Not upheld.** `names(x) <- <shorter>` pads with `NA` in R rather than erroring, and the
+>    trace is reproducible: `f_rename_par` returns
+>    `alpha0_1, alpha1_1, alpha0_2, alpha1_2, beta_2, P_1_1, P_2_1, NA` and `f_mapPar` then
+>    yields `0.1, 0.1, 80.008, 0.001, 0.1, 0.8001, 0.5, NA`. Codex's alternative also fails:
+>    in the working (no-`par0`) path the fixed parameter *is* respected — `beta_1` is exactly
+>    0.8 across all 300 draws. B3 stands as written.
+> 2. **The AIC/BIC test's only real fit is K = 2, where both dof formulas coincide** — fair,
+>    and already stated in that section; the K = 3/4 legs are deliberately synthetic because
+>    B2 makes a constrained K ≥ 3 fit impossible.
+> 3. **"All 11 blocks fail on stock 2.51" is not literally demonstrated** — it is, by
+>    measurement rather than by reading: the AIC/BIC block reports `pass=3 fail=6` on stock.
+>    Wording sharpened below.
+
+**`R CMD check --as-cran`: 2 WARNINGs, 1 NOTE — none of them is any of the defects below.**
+The WARNINGs are the usual pair (CRAN incoming feasibility: version not bumped + `Date`
+over a month old; and one spurious clang warning from R's own `R_ext/Boolean.h`). The NOTE is
+a genuine but harmless S3 signature mismatch on `Sim`/`Sim.MSGARCH_ML_FIT`. `tests/testthat/`
+(9 files as shipped) passes and covers none of this.
+
+**The conditional distributions are clean.** Unlike `GAS`, every density/CDF/quantile/moment
+checked out analytically: the standardised Student-*t* and GED constants, `E|z|` for all three
+families, the Fernández–Steel construction in `Skewed.h` (kernel, CDF, inverse-CDF and all four
+truncated moments `Eabsz`, `EzIpos`, `EzIneg`, `Ez2Ineg` — I re-derived each and they match,
+including both `xi >= 1` and `xi < 1` branches), the Hamilton filter and its over/underflow
+bookkeeping, Kim's smoother, the ergodic-distribution formula, the transition-matrix memory
+layout, and the tGARCH/gjrGARCH second-moment stationarity conditions. **The likelihood itself
+is correct.** The defects are in the layer around it: inference, persistence, forecasting for
+Bayesian fits, and the constrained-estimation paths.
+
+---
+
+## Tier 1 — wrong numbers reported to the user
+
+### A1. Every standard error is computed with the delta method transposed — `R/Inference.R:38` — **FIXED**
+
+```r
+mJacob     <- numDeriv::jacobian(f_mapPar, vPw_mod, spec = spec, do.plm = do.plm)
+mSandwitch <- t(mJacob) %*% mInvHessian %*% mJacob      # <- should be mJacob %*% ... %*% t(mJacob)
+```
+
+`numDeriv::jacobian(f, x)[i, j] = ∂f_i/∂x_j`, so with `θ_natural = g(θ_working)` the delta
+method is `Var(g) = J V J'`. The package computes `J' V J`. That is identical only when `J` is
+symmetric, and here it is not: the working→natural map is triangular within each regime
+(the sGARCH upper bound on `beta` is `0.9999 − alpha1`, gjr/tGARCH bounds on `beta` depend on
+`alpha1`, `alpha2` **and** on the shape/skew parameters), and the transition-probability block
+is *anti*-diagonal, because `f_mapGamma`/`f_unmapGamma`
+(`R/ParameterTransformation.R:180-213`) enumerate the off-diagonal entries in column-major
+order but attach the row-major parameter names to them. The Jacobian at the default
+MS(2)-GARCH(1,1)-Normal fit on `SMI`:
+
+```
+         alpha0_1 alpha1_1  beta_1 alpha0_2 alpha1_2 beta_2   P_1_1   P_2_1
+alpha1_1        0   0.0795  0.0000        0        0      0       0       0
+beta_1          0  -0.0767  0.0303        0        0      0       0       0
+beta_2          0        0  0.0000        0  -0.0054  5e-04       0       0
+P_1_1           0        0  0.0000        0        0      0  0.0000 -0.0212
+P_2_1           0        0  0.0000        0        0      0  0.0013  0.0000
+```
+
+**Evidence.** I computed the observed information directly in the natural parameterisation
+(central differences on `Kernel(spec, par, y, log = TRUE, do.prior = FALSE)`, which reproduces
+`fit$loglik` to machine precision) and compared:
+
+| | natural-scale `H` | delta method `J V J'` | **package `J' V J`** | package / correct |
+|---|---|---|---|---|
+| `alpha0_1` | 0.00744 | 0.00725 | 0.00725 | 1.00 |
+| `alpha1_1` | 0.01606 | 0.01511 | **0.03413** | **2.26×** |
+| `beta_1`   | 0.02197 | 0.02091 | **0.00958** | **0.46×** |
+| `alpha0_2` | 0.01688 | 0.01769 | 0.01769 | 1.00 |
+| `alpha1_2` | 0.00431 | 0.00437 | **0.00610** | **1.40×** |
+| `beta_2`   | 0.00416 | 0.00426 | **0.00049** | **0.12×** |
+| `P_1_1`    | 0.00982 | 0.00973 | **0.00059** | **0.06×** |
+| `P_2_1`    | (boundary, 0.9987) | 0.03072 | **0.50250** | **16×** |
+
+The natural-scale information agrees with `J V J'` to 2–6% on the seven interior parameters
+(`P_2_1` sits at 0.9987 where the natural-scale Hessian is not usable). The shipped numbers
+are wrong on **six of eight** parameters, by factors from 0.06× to 16×. Because the P-block of
+`J` is anti-diagonal, transposing literally **swaps the two transition probabilities' standard
+errors**: `summary(fit)` reports `P_1_1 = 0.978 (s.e. 0.0006)` and
+`P_2_1 = 0.9987 (s.e. 0.503)` when the correct values are ≈ 0.0097 and ≈ 0.031.
+
+Only `alpha0_k`, whose map is a plain `exp`, is unaffected. Every published table produced by
+`summary()` on a model with a `beta` or a transition probability is affected. **Fix**: one
+transpose. Verified.
+
+**After the fix** (`mSandwitch <- mJacob %*% mInvHessian %*% t(mJacob)`), the same fit:
+
+```
+         Estimate  Std. Error   t value  Pr(>|t|)   natural-scale H
+alpha0_1 0.021632    0.007245     2.986  0.002830          0.007444
+alpha1_1 0.087024    0.015111     5.759  0.000000          0.016063
+beta_1   0.881494    0.020912    42.152  0.000000          0.021975
+alpha0_2 0.020660    0.017687     1.168  0.242771          0.016884
+alpha1_2 0.005396    0.004373     1.234  0.217243          0.004309
+beta_2   0.994041    0.004256   233.588  0.000000          0.004161
+P_1_1    0.978348    0.009727   100.577  0.000000          0.009823
+P_2_1    0.998703    0.030720    32.509  0.000000    (boundary)
+```
+
+and across a sweep of eight specifications (single-regime; MS-2 sGARCH-Normal; MS-2
+gjrGARCH/tGARCH with sstd/std; MS-2 eGARCH-GED; MIX-2 sGARCH-std; MS-3; `fixed.pars`;
+`regime.const.pars`) every fit returns finite positive standard errors that reproduce
+`J V J'` to 1e-8 and p-values inside [0, 1].
+
+### A2. `Pr(>|t|)` is a one-sided p-value — `R/Inference.R:43` — **FIXED**
+
+```r
+vPvalues <- 1 - pnorm(abs(vTest))          # should be 2 * (1 - pnorm(abs(vTest)))
+```
+
+The column is labelled `Pr(>|t|)` and printed by `summary()`. Reported values are exactly half
+the two-sided p-value: on the `SMI` fit, `alpha1_1` shows `0.00539` where the two-sided value
+is `0.01077`, and `P_2_1` shows `0.0234` vs `0.0469`. Combined with A1 the two errors do not
+cancel — they compound.
+
+(Separately: t-tests against zero are meaningless for `beta`, for `nu`, and for the transition
+probabilities, all of which have non-zero-centred supports. A note in `?FitML` would help.)
+
+**After the fix**: the `Pr(>|t|)` column equals `2 * (1 - pnorm(abs(t)))` exactly
+(`all.equal` TRUE) on every specification in the sweep.
+
+### A3. `predict()` on an MCMC fit returns draw #1, not the posterior mean — `R/CondVol.R:32` — **FIXED**
+
+```r
+vol <- matrix(NA, nrow = dim(PredProb)[1], ncol = nrow(par.check))
+...
+tmp <- mean(vol[dim(PredProb)[1]])          # linear index -> row T+1 of COLUMN 1 only
+```
+
+`vol` is `(T+1) × ndraw`. A single index into a matrix is linear indexing, so `vol[T+1]` is the
+last row of the *first* column. `mean()` of a scalar is a no-op. The intended expression is
+`mean(vol[dim(PredProb)[1], ])`.
+
+Reproduced on `FitMCMC(CreateSpec(), SMI, ctr = list(nburn=500, nmcmc=1000, nthin=10))`,
+100 retained draws:
+
+```
+predict(mc, nahead = 1)$vol : 1.019753
+draw #1 only                : 1.019753   <- exact match
+posterior mean over draws   : 1.039776
+range across draws          : 0.999837 – 1.097337
+```
+
+The reported one-step-ahead volatility is whichever value the first retained draw happens to
+give. `Volatility()` (the in-sample path) averages across draws correctly, so the two functions
+are mutually inconsistent on the same fit. This also feeds the `h = 1` row of
+`predict(..., nahead = h)` and of `UncVol()` for MCMC fits.
+
+**After the fix**, same chain: `predict(mc, nahead = 1)$vol = 1.039776`, exactly the posterior
+mean, and no longer draw #1's `1.019753`. The ML path is unchanged (one column, so
+`mean(vol[N, ])` and the old `vol[N]` coincide): `predict(fit, nahead = 1)$vol = 1.030426`
+before and after.
+
+### A4. `AIC`/`BIC` degrees of freedom undercount `regime.const.pars` — `R/Utils.R:420` — **FIXED**
+
+```r
+dofMSGARCH = function(object){
+  length(object$spec$par0) - length(object$spec[["regime.const.pars"]]) - length(object$spec[["fixed.pars"]])
+}
+```
+
+Each regime-constant parameter removes `K − 1` free parameters, not one (`f_rename_par` strips
+`name_2 … name_K`). Correct only at `K = 2`:
+
+```
+K=2  par0=10  truly free= 9  dofMSGARCH= 9   ok
+K=3  par0=18  truly free=16  dofMSGARCH=17   MISMATCH
+K=4  par0=28  truly free=25  dofMSGARCH=27   MISMATCH
+```
+
+`stats::AIC`/`BIC` use this via `logLik.MSGARCH_ML_FIT`, so every K ≥ 3 model-selection table
+built with `regime.const` is penalised wrongly — and always in the direction that *disfavours*
+the constrained model. Fix: `- length(regime.const.pars) * (K - 1)`.
+
+**After the fix**: `K=2 -> 9`, `K=3 -> 16`, `K=4 -> 25`, all matching the true free-parameter
+count. The unconstrained (`8`) and `fixed.pars` (`7`) cases are unchanged.
+
+---
+
+## Tier 2 — hard failures
+
+### B1. A saved spec or fit cannot be reloaded, and the recovery path that exists for exactly this is dead code — `R/Utils.R:340-359` — **FIXED**
+
+`MSGARCH_SPEC` holds Rcpp module objects, so `saveRDS`/`readRDS` across sessions leaves a stale
+external pointer. `f_check_spec` exists to detect that and rebuild — but the rebuild branch
+re-dereferences the *same dead pointer* on its second line:
+
+```r
+is.OK = tryCatch({ spec$rcpp.func$get_sd(); TRUE }, error = function(e) FALSE)
+if (!isTRUE(is.OK)) {
+  spec.new   = f_spec(models = spec$name, do.mix = spec$is.mix)
+  prior.mean = spec$rcpp.func$get_mean()     # <- line 349: dead pointer again, uncaught
+  prior.sd   = spec$rcpp.func$get_sd()       # <- line 350: same
+  ...
+  spec$rcpp.func$set_mean(spec$prior.mean)   # <- and then uses spec$prior.mean, not the
+  spec$rcpp.func$set_sd(spec$prior.sd)       #    prior.mean/prior.sd just computed
+}
+```
+
+In a fresh session, on a `fit` written with `saveRDS`:
+
+```
+Volatility(fit)   : Error in .External(...): NULL value passed as symbol address
+State(fit)        : Error in .External(...): NULL value passed as symbol address
+predict(fit)      : Error in .External(...): NULL value passed as symbol address
+FitML(saved spec) : Error in .External(...): NULL value passed as symbol address
+```
+
+**Fix applied**: drop lines 349–353 (they recompute, from the dead pointer, exactly the named
+vectors already stored in `spec$prior.mean` / `spec$prior.sd`, and then throw the result away),
+and keep lines 355–356, which already read the R-side copies. In a fresh session on a
+`readRDS`-ed fit, `Volatility`, `State`, `predict`, `summary`, `Risk`, `PIT`, `simulate` and
+`AIC` all now succeed, and return values identical to the in-session ones (`Volatility` head
+`1.2045, 1.2530, 1.2372`; `predict` h=1 `1.030426` both ways). Custom priors survive the
+rebuild — a spec created with `prior = list(mean = list(beta_1 = 0.7), sd = list(beta_1 = 0.1))`
+comes back with mean `0.7` / sd `0.1` on `beta_1` and defaults elsewhere. Saved
+`MSGARCH_MCMC_FIT` objects (`predict`, `DIC`) work too. This unblocks the single most common
+workflow in applied use: fit overnight, save, analyse later; `parLapply` over a rolling window;
+caching a fit in a knitr chunk.
+
+### B2. Singular ergodic-distribution inverse aborts the whole fit — `src/MSgarch.h:277-286`
+
+`loadparam` recomputes the stationary distribution on *every* likelihood evaluation with a raw
+Armadillo inverse:
+
+```cpp
+arma::mat foo   = (I - as<arma::mat>(P_mat) + Umat).t();
+arma::vec delta = (foo).i() * Uvec;          // throws if singular
+```
+
+With `do.plm = TRUE` — forced whenever `fixed.pars` or `regime.const.pars` is set, and
+hard-coded in `FitMCMC` (`R/FitMCMC.R:134`) — the free transition entries are mapped into
+`(0,1)` *independently*, so for `K ≥ 3` a row can sum to more than 1 and `extract_P_it` then
+produces a negative last entry. `I − P + U` can be exactly singular, and the Armadillo
+exception propagates out of `Kernel`, out of `f_nll`, and kills the run.
+
+Reproduced deterministically. `CreateSpec(model="sGARCH", distribution="std", K=3,
+constraint.spec=list(regime.const="nu"))` + `FitML(..., SMI)` with `set.seed(1)`; instrumenting
+`f_nll` to record the last parameter point reached by BFGS gives
+
+```
+     [,1]     [,2]      [,3]
+[1,]    1 0.999973 -0.999973
+[2,]    0 1.000000  0.000000
+[3,]    0 0.000000  1.000000
+det(I - P + U) = 0    rcond = 0
+Kernel  -> ERROR: matrix multiplication: problem with matrix inverse
+```
+
+`K = 3` *without* `regime.const` fits fine (the `do.plm = FALSE` map always yields a proper
+stochastic matrix); `K = 2` *with* `regime.const` fits fine (one free entry per row cannot
+overflow the simplex). So the failure is specific to `K ≥ 3` on the `do.plm` path — which
+includes **all** `FitMCMC` runs with three or more regimes, where a single unlucky proposal
+destroys the entire chain with no partial output.
+
+Fix: guard the inverse (`arma::solve` with `arma::solve_opts::no_error`, or check
+`arma::inv(...)`'s bool return) and return the model's own "infeasible" signal (`-1e10`) instead
+of throwing; the constraint `all(0 < P_it < 1)` already exists in `calc_prior` but is evaluated
+*after* `loadparam`, so it never gets the chance to reject the point.
+
+### B3. `FitMCMC` mis-aligns the parameter vector when `fixed.pars` is combined with `ctr$par0` — `R/FitMCMC.R:140-151`
+
+`FitML` removes the fixed entries from the starting vector (`f_remove_fixedpar`,
+`R/FitML.R:126`); `FitMCMC` substitutes their values but **never removes them**, so `par0` is
+passed to the sampler with `d` elements when the sampler expects `d − n_fixed`. `f_rename_par`
+then labels the over-long vector with the short name list, leaving a trailing `NA` name, and
+`f_mapPar` looks the bounds up by name — so every parameter after the fixed one is mapped with
+the *wrong* bounds and the last one becomes `NA`:
+
+```
+spec: 8 parameters, beta_1 fixed -> sampler expects 7
+length(par0) inside FitMCMC = 8
+names from f_rename_par     : alpha0_1,alpha1_1,alpha0_2,alpha1_2,beta_2,P_1_1,P_2_1,NA
+f_mapPar gives              : 0.1, 0.1, 80.008, 0.001, 0.1, 0.8001, 0.5, NA
+FitMCMC(..., ctr=list(par0=spec$par0)) -> Error : matrix multiplication: problem with matrix inverse
+```
+
+`alpha0_2` becomes 80, `alpha1_2` becomes 0.001, and the `NA` then detonates the inverse of
+B2. Without a user `par0` the same spec samples fine (`accept = 0.287`), so this is purely
+the `ctr$par0` branch. Fix: add `par0 <- f_remove_fixedpar(par0, spec$fixed.pars)` after
+line 148, mirroring `FitML`.
+
+### B4. `constraint.spec = list(fixed = list(P_i_j = …))` is accepted but unusable
+
+`CreateSpec` validates fixed-parameter names against `out$label`, which includes `P_1_1`, so
+fixing a transition probability is accepted. Two independent things then break:
+
+1. `f_recover_fixedpar_SR` (`R/ParameterConstraints.R:34-56`) splits fixed parameters by
+   regime with `gsub("_k", "", name)`; `"P_1_1"` matches the regime-1 test, is handed to a
+   *single-regime* `CreateSpec` whose labels are only `alpha0_1, alpha1_1, beta_1`, and dies:
+   `Wrong name in fixed.pars: P_1_1`. So `FitML(spec, data)` and `FitMCMC(spec, data)` both
+   fail before the optimiser starts. (With an explicit `ctr$par0`, `FitML` bypasses the
+   starting-value routine and works — `loglik = -3390.59`.)
+2. The prior correction in `Kernel` (`R/Kernel.R`) subtracts
+   `dnorm(par[, names(fixed.pars)], prior.mean[names(fixed.pars)], …)`, but `prior.mean` only
+   covers the `sum(NbParams)` within-regime coefficients — the transition probabilities have a
+   uniform prior and no entry. `prior.mean["P_1_1"]` is `NA`, the whole log-posterior becomes
+   `NA` and is floored to `-1e10`:
+
+```
+length(prior.mean) = 6, length(label) = 8   (P_1_1, P_2_1 not covered)
+log-posterior  (do.prior=TRUE)  : -1e+10
+log-likelihood (do.prior=FALSE) : -3487.971
+```
+
+so even reaching the sampler would give a frozen chain. Either reject `P_*` in
+`f_check_parameterConstraints` with a clear message, or handle it in both places.
+
+### B5. The `FitML` failure guard has the wrong sign and is unreachable anyway — `R/FitML.R:129-137`
+
+```r
+optimizer <- ctr$OptimFUN(vPw, f_nll, spec, data_, ctr$do.plm)
+llk <- -optimizer$value
+if (llk == 1e+10) { f_error("FitML -> Error during optimization"); stop() }
+```
+
+`f_nll` returns `+1e10` on failure, so `llk` is `−1e10`; the test can never be true. And
+`f_OptimFUNDefault` wraps `optim` in `try()`, so on any error `optimizer` is a character
+`try-error` and `optimizer$value` throws first. Every failure mode therefore surfaces as
+
+```
+Error in optimizer$value : $ operator is invalid for atomic vectors
+```
+
+— which is what a user sees for B2, for B4, and for a single `NA` in the data
+(`f_check_y`, `R/Utils.R:177`, only rejects data that is *entirely* `NaN`, so
+`y[100] <- NA` sails through and produces this message). Three unrelated problems, one
+uninterpretable error. Fix: test `inherits(optimizer, "try-error")` first, then `llk == -1e10`;
+and make `f_check_y` reject `any(!is.finite(y))`.
+
+---
+
+## Tier 3 — silently wrong results in specific calls
+
+### C1. `pdf_Rcpp` / `cdf_Rcpp` with `is_log = TRUE` return the last regime, not the mixture — `src/MSgarch.h:397-401, 472-476`
+
+```cpp
+for (many::iterator it = specs.begin(); it != specs.end(); ++it) {
+  for (int i = 0; i < nx; i++) {
+    tmp[i] = (*it)->spec_calc_pdf(x[i] / sig) / sig;
+    out[i] = out[i] + tmp[i] * PLast[s];          // out = correct mixture
+  }
+  s++;
+}
+if (is_log) { for (int i = 0; i < nx; i++) out[i] = log(tmp[i]); }   // tmp = LAST regime only
+```
+
+`tmp` is overwritten each regime, so the mixture in `out` is thrown away and replaced by the
+log density of regime *K* alone:
+
+```
+x                        -3       -1        0        1        3
+mixture pdf         0.00561  0.24042  0.39670  0.24042  0.00561
+returned (log=TRUE) -2.68223 -1.77355 -1.65997 -1.77355 -2.68223
+exp(returned)/pdf     12.199    0.706    0.479    0.706   12.199   <- should be all 1
+cdf, same ratio       26.480    1.972    1.000    0.814    0.926
+```
+
+Not reachable from `PredPdf`/`PIT` (both always pass `FALSE` and take the log in R), but these
+are live methods on `spec$rcpp.func` and the flag is part of the C++ signature. Fix:
+accumulate into a scratch vector and take `log(out[i])`.
+
+### C2. `MSgarch::f_cdf_its` writes the first observation transposed — `src/MSgarch.h:497`
+
+```cpp
+tmp(ix, 0, s) = (*it)->spec_calc_cdf(x(ix, 0) / sig);   // t=0 block
+...
+tmp(i, ix, s) = (*it)->spec_calc_cdf(x(ix, i) / sig);   // t>=1 loop — correct orientation
+```
+
+The cube is `(ny, nx, K)`. `SingleRegime::f_cdf_its` and both `f_pdf_its` write `(0, ix, s)`;
+only the multi-regime CDF got it backwards. Consequences: the `t = 1` row is left at zero for
+every grid point except the first, and if `nx > T` the write runs off the cube.
+
+```
+cdf_its, t=1, regime 1: 0.007916 0 0 0 0        <- zeros
+cdf_its, t=2, regime 1: 0.013265 0.133676 0.5 0.866324 0.986735
+pdf_its, t=1, regime 1: 0.0262 0.232471 0.481275 0.232471 0.0262   (pdf path is fine)
+
+PIT(spec, x = c(-2,-1,0,1,2), par, data = SMI, do.its = TRUE)[1:2, ]
+1990-11-12  0.015613  0.000000  0.0  0.000000  0.000000     <- t=1 wrong
+1990-11-13  0.020841  0.140055  0.5  0.859945  0.979159
+
+PIT(..., x = <length T+10>, do.its = TRUE) -> Error: Cube::operator(): index out of bounds
+```
+
+Harmless in the common `PIT(fit, do.its = TRUE)` call (there `nx = 1`), wrong whenever a user
+supplies an evaluation grid.
+
+### C3. `UncVol` averages the burn-in window and ignores the converged tail — `R/UncVol.R:88` + `R/Utils.R:160`
+
+The documentation says `nahead = 5000L`, `nburn = 1000L`, "simulating nsim paths up to
+`nburn + nahead` … discarding the first `nburn` … and computing the mean of the remaining".
+The code has the two defaults **swapped** (`f_process_ctr(type = 2)`:
+`nburn = 5000L, nahead = 1000L`) and averages
+
+```r
+out <- mean(tmp[ctr$nburn:ctr$nahead])        # should be mean(tmp[(nburn+1):(nburn+nahead)])
+```
+
+With the shipped defaults that is `tmp[5000:1000]` — a *descending* sequence covering horizons
+1000–5000 out of the 6000 simulated, i.e. the best-converged 1000 horizons are the ones thrown
+away. It happens to land far enough out that the default answer is roughly right; it stops
+being right the moment a user passes their own control list:
+
+```
+UncVol(fit, ctr = list(nsim=100, nburn=400, nahead=100))
+  what the code averages   (horizons 100..400) : 1.36344
+  what the docs promise    (horizons 401..500) : 1.06793      <- 28% apart
+two calls, identical settings                  : 1.36344 / 1.33357
+```
+
+The last line is a second issue: `UncVol` is a Monte Carlo estimate with no `seed` argument, so
+it is not reproducible, and it is driven from a two-point fake sample `data = c(1, 1)`
+(`R/UncVol.R:84`). Worth documenting at minimum.
+
+### C4. `State()$Viterbi` uses a transposed transition matrix for mixtures — `R/State.R:40`
+
+```r
+P <- TransMat(object, par = par[i, ], nahead = 1)          # 1 x K row of mixture weights
+if (isTRUE(object$is.mix)) P <- matrix(rep(P, object$K), nrow = object$K, ncol = object$K)
+```
+
+`matrix()` fills column-major, so row *i* becomes `(p_i, …, p_i)` instead of `(p_1, …, p_K)`.
+With mixture weights `(0.8, 0.2)`:
+
+```
+passed to Viterbi():   [,1] [,2]        correct:   [,1] [,2]
+                 [1,]   0.8  0.8                    0.8  0.2
+                 [2,]   0.2  0.2                    0.8  0.2
+row sums: 1.6, 0.4   <- not a stochastic matrix
+```
+
+The decoded path for every `do.mix = TRUE` model is computed against this. Fix: `byrow = TRUE`.
+(Filtered/smoothed/predicted probabilities are computed in C++ from the correct matrix and are
+unaffected — only `$Viterbi`.)
+
+### C5. `CreateSpec`'s `K`-expansion guard reads a field that does not exist — `R/CreateSpec.R:213`
+
+```r
+if (length(variance.spec$model) > 1 | length(distribution.spec$model) > 1)
+#                                                             ^^^^^ should be $distribution
+```
+
+`distribution.spec$model` is always `NULL`, so `length(NULL) > 1` is always `FALSE` and the
+check never fires for the distribution vector. The subsequent `rep(...)` then recycles and
+`distribution[1:length(model)]` truncates:
+
+```r
+CreateSpec(variance.spec = list(model = "sGARCH"),
+           distribution.spec = list(distribution = c("std","norm")),
+           switch.spec = list(K = 3))
+$name  "sGARCH_std"  "sGARCH_norm"  "sGARCH_std"
+```
+
+The user asked for something contradictory and silently got a scrambled three-regime spec
+instead of the intended error.
+
+---
+
+## Tier 4 — minor, hygiene, and dead code
+
+- **`R/Utils.R:294` `f_getGamma`** scrambles the transition matrix for `K ≥ 3` (it reads the
+  row-major parameter vector with a column-major `matrix()` fill). Currently **unused** — but
+  it is a landmine for anyone who reaches for it.
+- **`R/Utils.R:2` `f_GammaParNames`** emits `P_1_1, P_2_1, P_3_1, P_1_2, …` (column-major),
+  whereas the labels actually built by `f_spec:398-401` and consumed by
+  `MSgarch::extract_P_it` are `P_1_1, P_1_2, P_2_1, …` (row-major). Also unused; delete or fix,
+  don't leave two contradictory definitions of the parameter order in the same file.
+- **`src/Utils.h:42`** `adaptiveSimpsonsAux` guards recursion with `bottom <= 100` instead of
+  `bottom <= 0`, so any call with `maxRecursion ≤ 100` returns after a single Simpson step with
+  no adaptation. Unreachable — `adaptiveSimpsons` is never called — but it should either be
+  fixed or removed. (The integrator actually in use, `Skewed::compositeSimpsons` with
+  `Nsi = 5`, is fine: 10 panels over a short interval on a smooth integrand.)
+- **`src/MSgarch.h:71`** `P_mean = 1 / K` is integer division on an `int` — always 0 for
+  `K ≥ 2`. Both `P_mean` and `P_sd` are write-only; delete them.
+- **`src/MSgarch.h:768,774`** `eval_model` calls `calc_prior(theta_j)` twice per parameter
+  vector and discards the first result. Cost is small next to the filter (measured ~0.8 ms per
+  kernel evaluation at `T = 2500`), but it is free to remove.
+- **`src/adaptMCMC.cpp:82`** a full `eig_sym` is computed and immediately discarded every MCMC
+  iteration (`makePositiveDefinite` recomputes it on line 83). Also `sigma.eye(l_param,
+  l_param)` on line 81 zeroes `sigma` *in place* mid-expression; it happens to be harmless
+  because `S` was cached on line 77, but it is a trap for the next editor.
+- **`R/CreateSpec.R:312`** validates `prior.sd` with `f_check_parameterPriorMean`, so a bad
+  name in `prior$sd` reports "Wrong name in prior.mean".
+- **`R/Utils.R:177`** `f_check_y` only rejects data that is entirely `NaN`; a single `NA`
+  passes and surfaces later as the B5 error.
+- **`logLik.MSGARCH_ML_FIT`** sets `nobs = length(data)`, but `calc_lndMat` evaluates `T − 1`
+  terms (observation 1 initialises the recursion). `BIC` therefore uses `log(T)` rather than
+  `log(T − 1)` — negligible, but worth a line in `?FitML` since it also explains the
+  `State()$Viterbi[1] = Viterbi[2]` fudge at `R/State.R:99`.
+- **`eGARCH::set_vol`** initialises `h = exp(α₀/(1−β))`, i.e. `exp(E[ln h])`, and
+  **`tGARCH::set_vol`** initialises `h = (E[σ])²`, not `E[h]`. Defensible as a recursion
+  starting value, but `UncVol`'s per-regime `unc_vol_Rcpp` inherits the same convention and
+  therefore reports something other than the unconditional variance for those two models.
+  Document it.
+- **Dead/unused**: `f_map_deriv`, `f_rbindrep`, `f_getGamma`, `f_GammaParNames`,
+  `Decoding_HMM` (exported to R, called nowhere), `TransMat`'s `object$is.shape.ind` branch
+  (`is.shape.ind` is never set by `CreateSpec`), and the `UnmapParameters_univ(x, "norm",
+  FALSE)` path in `Mapping.cpp`, which returns an uninitialised vector — currently guarded by
+  the `vDist[k] != "norm"` test at `R/StartingValues.R:189`, so unreachable.
+- **`R CMD check` NOTE**: register `Sim` methods or align their signatures with the generic.
+  The two WARNINGs are the version/date pair and a clang artefact, both cosmetic.
+
+---
+
+## Suggested order of work
+
+1. ~~`R/Inference.R:38` (transpose) and `:43` (×2) — one line each, fixes every published
+   standard error and p-value. **A1/A2.**~~ **DONE**
+2. ~~`R/Utils.R:349-350` — two lines, makes saved fits reloadable. **B1.**~~ **DONE**
+3. ~~`R/CondVol.R:32` — add the missing comma. **A3.**~~ **DONE**
+4. ~~`R/Utils.R:420` — `* (K - 1)`. **A4.**~~ **DONE**
+5. Guard the inverse in `src/MSgarch.h:283` and return `-1e10` instead of throwing; then fix
+   the `FitML` guard at `R/FitML.R:133`. **B2/B5.**
+6. `R/FitMCMC.R:148` + `f_check_parameterConstraints` — the two `fixed.pars` paths. **B3/B4.**
+7. `src/MSgarch.h:497`, `:399`, `:474`; `R/State.R:40`; `R/UncVol.R:88`;
+   `R/CreateSpec.R:213`. **C1–C5.**
+
+Items 1–4 were seven lines of R and are done, with regression cover added for the two
+most damaging defects:
+
+- **`tests/testthat/test_Inference.R`** — pins the delta method against an independently
+  computed natural-scale observed information (single-regime GARCH-Normal, where every
+  parameter is interior and the numerical Hessian is well conditioned: the correct sandwich
+  agrees to 1e-4 relative, the transposed one is off by 44% and 85%), asserts the shipped
+  standard errors equal `sqrt(diag(J V J'))` to 1e-8 on both a single-regime and an MS(2)
+  fit, and checks that `Pr(>|t|)` is two-sided. Each orientation test first asserts that the
+  two orientations actually *differ* for that model, so it cannot silently go vacuous.
+- **`tests/testthat/test_Serialization.R`** — round-trips a spec, an `MSGARCH_ML_FIT` and an
+  `MSGARCH_MCMC_FIT` through `saveRDS`/`readRDS` (R restores external pointers as `NULL`, so
+  this reproduces the cross-session failure within one session) and requires `Volatility`,
+  `State`, `predict`, `AIC`, `summary`, `DIC` and the log-kernel to return values identical
+  to the originals, with user-supplied priors preserved. It too opens with a guard asserting
+  that the round trip really did invalidate the pointers.
+
+- **`tests/testthat/test_Inference.R`**, two further blocks for **A4** — `dofMSGARCH` for
+  K = 2, 3, 4 with a regime-constant `nu` (K = 2 cannot separate "one" from "K − 1", which is
+  why the shipped code looked right), plus the unconstrained and `fixed.pars` cases as
+  invariants, and an end-to-end check that `logLik()`'s `df` attribute and the `AIC`/`BIC`
+  arithmetic follow it. The K = 3 and K = 4 legs run against a fit-shaped list rather than a
+  real fit, because a constrained K ≥ 3 model cannot currently be fitted at all — that is
+  **B2**, still open.
+- **`tests/testthat/test_Volatility.R`**, one appended block for **A3** — recomputes the
+  one-step-ahead volatility draw by draw through the public interface (each call carries a
+  single parameter vector, so it cannot depend on how the draws are pooled) and requires
+  `predict()` on the MCMC fit to equal their mean to 1e-10. Guarded by an assertion that the
+  posterior mean is distinguishable from the first draw.
+
+11 test blocks, 34 assertions, +13s of check time. All pass on the patched package and all 11
+fail on stock 2.51 — block-level, measured, not inferred. Within the failing blocks some
+individual assertions still pass on stock, by construction: the AIC/BIC block reports
+`pass=3 fail=6` there because its K = 2 leg is exactly the case the two dof formulas agree on,
+and the dof block reports `pass=3 fail=2` because K = 2 plus the two invariants are meant to
+hold either way. The only *additions* that pass wholly on both are the deliberate
+preconditions (the round trip really invalidates the pointers; the two sandwich orientations
+really differ; the posterior mean really differs from draw #1), which exist so the tests
+cannot quietly go vacuous.
diff --git a/REVIEW_codex.md b/REVIEW_codex.md
new file mode 100644
index 0000000..ae98c4a
--- /dev/null
+++ b/REVIEW_codex.md
@@ -0,0 +1,148 @@
+# Independent Codex cross-check of REVIEW.md
+
+Run 2026-08-10, `codex-cli 0.139.0`, read-only over `package MSGARCH/`:
+
+```bash
+codex exec --sandbox read-only --skip-git-repo-check "<audit prompt>"
+```
+
+Codex audited (1) the 14 load-bearing claims, (2) the five applied fixes, and (3) the new
+tests, by reading the source only — it could not run R. Verbatim output below; see the
+"Independent cross-check" section of `REVIEW.md` for how the three pushbacks were resolved.
+
+---
+
+## Part 1 — claim verdicts
+
+| # | claim | verdict | evidence (file:line) |
+|---|---|---|---|
+| A1 | delta-method sandwich transposed; P-block anti-diagonal | CONFIRM | Original code uses `t(mJacob) %*% mInvHessian %*% mJacob` while `jacobian(f_mapPar, ...)` returns output-by-input derivatives ([`Inference.R:36-38`](Package/R/Inference.R:36)). The P block is anti-diagonal for `K=2`: `f_mapGamma()` fills off-diagonals by `mGamma[!mGamma]` then returns `c(t(mGamma[, -K]))`, keeping the original names ([`ParameterTransformation.R:204-209`](Package/R/ParameterTransformation.R:204)). |
+| A2 | p-value is one-sided under `Pr(>|t|)` | CONFIRM | The column is named `Pr(>|t|)` at construction, but the value is `1 - pnorm(abs(vTest))` ([`Inference.R:4-5`, `:43`](Package/R/Inference.R:4)). |
+| A3 | MCMC forecast uses linear index draw #1 | CONFIRM | `vol` is allocated as `(T+1) x nrow(par.check)` and line 32 indexes it as `vol[dim(PredProb)[1]]`, a single linear index ([`CondVol.R:19`, `:32`](Package/R/CondVol.R:19)). |
+| A4 | `dofMSGARCH` subtracts one instead of `K-1` per regime-constant parameter | CONFIRM | `dofMSGARCH` subtracts `length(regime.const.pars)` only, while `f_rename_par()` removes `name_2 ... name_K` for each constrained base name ([`Utils.R:420-421`, `:364-367`](Package/R/Utils.R:420)). |
+| B1 | stale-pointer recovery re-dereferences dead pointer; lines 349-353 unused | CONFIRM | The failed branch calls `spec$rcpp.func$get_mean()` and `get_sd()` immediately after `get_sd()` failed; computed `prior.mean/prior.sd` are not used by lines 355-356 ([`Utils.R:340-356`](Package/R/Utils.R:340)). |
+| B2 | raw `.i()` can abort under invalid K>=3 PLM transition rows; K=2 nonsingular | CONFIRM | `loadparam()` forms rows via `extract_P_it()` with last entry `1 - sum(P_it)`, then calls `(foo).i()` without a guard ([`MSgarch.h:193-198`, `:281-283`](Package/src/MSgarch.h:193)). For `K=2`, `det(I-P+U)=2(1+P_2_1-P_1_1)>0` when both free entries are in `(0,1)`. |
+| B3 | FitMCMC custom `par0` fixed-parameter path too long/mis-mapped | PARTIALLY CONFIRM | It is true FitMCMC substitutes fixed values but does not remove them before `f_unmapPar()` ([`FitMCMC.R:143-150`](Package/R/FitMCMC.R:143)), unlike FitML’s `f_remove_fixedpar()` ([`FitML.R:125-126`](Package/R/FitML.R:125)). But the exact “trailing NA name” mechanism is overstated: with names present the fixed dimension remains in the sampler, and with names absent `f_rename_par()` is more likely to fail length/name assignment than silently trail `NA` ([`Utils.R:370-374`](Package/R/Utils.R:370)). |
+| B4 | fixed transition probability accepted but breaks fitters by two mechanisms | CONFIRM | `CreateSpec()` validates against `out$label`, which includes transition labels ([`CreateSpec.R:278`, `:397-400`](Package/R/CreateSpec.R:397)); starting values pass recovered fixed names into single-regime `CreateSpec()` ([`StartingValues.R:201-206`](Package/R/StartingValues.R:201)). Separately, prior correction indexes `prior.mean[names(fixed.pars)]`, but MS priors come only from within-regime model params ([`Kernel.R:50-57`; `MSgarch.h:146-171`](Package/R/Kernel.R:50)). |
+| B5 | `llk == 1e+10` check wrong and try-error not handled | CONFIRM | `f_nll()` returns `+1e10` on likelihood failure, so `llk <- -optimizer$value` is `-1e10`; if `optim` itself fails, `f_OptimFUNDefault()` returns a `"try-error"` and `optimizer$value` is invalid ([`Likelihood.R:20-24`; `FitML.R:129-136`; `Optimizers.R:2-4`](Package/R/FitML.R:129)). |
+| C1 | MS `f_pdf`/`f_cdf` log path logs last-regime tmp, not mixture | CONFIRM | Both functions accumulate into `out`, then under `is_log` assign `out[i] = log(tmp[i])`, where `tmp` is overwritten in the regime loop ([`MSgarch.h:387-400`, `:462-475`](Package/src/MSgarch.h:387)). |
+| C2 | `f_cdf_its` first slice has swapped cube indices | CONFIRM | The initial CDF slice writes `tmp(ix, 0, s)` while later rows and `f_pdf_its` write `tmp(i, ix, s)` / `tmp(0, ix, s)` ([`MSgarch.h:490-507`; compare `:420-435`](Package/src/MSgarch.h:490)). |
+| C3 | `UncVol` window wrong and defaults swapped vs docs | CONFIRM | Docs say `nahead=5000`, `nburn=1000`, but `f_process_ctr(type=2)` sets `nburn=5000`, `nahead=1000`; `UncVol` then averages `tmp[ctr$nburn:ctr$nahead]` ([`UncVol.R:16-21`, `:88`; `Utils.R:159-160`](Package/R/UncVol.R:16)). |
+| C4 | mixture Viterbi transition matrix transposed/row sums wrong | CONFIRM | For mixtures, `TransMat()` returns a one-row probability vector; `State()` expands with `matrix(rep(P, K), nrow=K, ncol=K)`, which fills columns, so rows become repeated single probabilities rather than the mixture vector ([`State.R:90-94`; `Transmat.R:53-58`](Package/R/State.R:90)). |
+| C5 | K-expansion guard checks nonexistent `distribution.spec$model` | CONFIRM | The guard tests `length(distribution.spec$model) > 1`, but the checked object uses `$distribution` ([`CreateSpec.R:205`, `:212-219`](Package/R/CreateSpec.R:205)). |
+
+On the negative claims: I found no source-level contradiction to REVIEW.md’s claim that the base conditional densities/CDFs/quantiles, Fernandez-Steel moments, Hamilton filter, Kim smoother, valid-matrix ergodic formula, transition-memory layout, and tGARCH/gjrGARCH stationarity formulas are correct. Caveat: the ergodic formula is mathematically fine for valid stochastic matrices, but B2 shows the implementation does not guard invalid PLM matrices before inverting.
+
+## Part 2 — fixes
+
+A1: `mJacob %*% mInvHessian %*% t(mJacob)` is the correct orientation for `numDeriv::jacobian(f_mapPar, ...)`. The boundary adjustment changes the evaluation point to a nearby interior `vPw_mod`; it does not alter the orientation argument, and the Hessian is also computed at `vPw_mod` when recomputed. OK
+
+A2: `2 * (1 - pnorm(abs(vTest)))` matches the printed `Pr(>|t|)` label. The broader statistical caveat remains that zero-null z/t tests are not meaningful for all bounded/support-shifted parameters, but this fix corrects the advertised two-sided calculation. OK
+
+A3: `mean(vol[dim(PredProb)[1], ])` preserves ML behavior because a one-draw matrix row averages to the same scalar. Averaging posterior volatility rather than variance is a modeling/reporting choice already implied by `Volatility()`’s row-mean-on-volatility path ([`CondVol.R:53-55`](Package/R/CondVol.R:53)); it is separate from the indexing bug. OK
+
+A4: `- length(regime.const.pars) * (K - 1L)` is right for allowed specs: `K=1` rejects non-null regime constants, and `CreateSpec()` forbids using `fixed.pars` and `regime.const.pars` together ([`ParameterConstraints.R:115-121`; `CreateSpec.R:293-294`](Package/R/CreateSpec.R:293)). It also works for mixtures because the constrained duplicated within-regime parameters are still K copies collapsed to one. OK
+
+B1: Rebuilding only `rcpp.func` is sufficient from this source: all external-pointer-backed callables are collected there by `f_spec()` ([`CreateSpec.R:350-366`](Package/R/CreateSpec.R:350)), while `func` contains plain R closures. The `!is.null` guards should not affect normal `f_spec()` / `ExtractStateFit()` specs because `prior.mean` and `prior.sd` are initialized non-null ([`CreateSpec.R:412-416`](Package/R/CreateSpec.R:412)). OK
+
+## Part 3 — tests
+
+`test_Inference.R` “Standard errors use J V J'”: fails on original, has a non-vacuous guard, but partly restates implementation via `fit$Inference$Hessian`. OK
+
+`test_Inference.R` “observed information in natural scale”: independent and would fail on original; finite-difference tolerance and full `SMI` ML fit make it somewhat machine/runtime sensitive. WEAK
+
+`test_Inference.R` “Transition-probability standard errors”: fails on original and directly covers the P-block; uses fixed historical parameters, so it is reasonably stable. OK
+
+`test_Inference.R` “Pr(>|t|)”: fails on original and is not vacuous. OK
+
+`test_Inference.R` “Degrees of freedom”: catches original for `K=3/4`; independent enough. OK
+
+`test_Inference.R` “AIC and BIC”: the real fitted case is only `K=2`, where old and new formulas coincide; affected `K=3/4` cases are synthetic objects that test plumbing, not fitting. WEAK
+
+`test_Serialization.R` “round trip invalidates pointers”: useful anti-vacuity guard, but depends on current Rcpp serialization behavior. OK
+
+`test_Serialization.R` “reloaded spec usable”: fails on original and checks outputs against a live spec. OK
+
+`test_Serialization.R` “keeps priors”: fails on original and checks rebuilt C++ prior state. OK
+
+`test_Serialization.R` “reloaded ML fit usable”: fails on original; runtime acceptable but includes a full ML fit. OK
+
+`test_Serialization.R` “reloaded MCMC fit usable”: fails on original but is RNG/runtime-sensitive and uses a very short chain. WEAK
+
+`test_Volatility.R` “MCMC forecast averages draws”: should fail on original and has a non-vacuity guard; fragile because it runs MCMC in unit tests and depends on the first draw differing from the mean. WEAK
+
+## Missed or overstated
+
+B3 is overstated. The core bug is real: FitMCMC leaves the fixed parameter in user-supplied `ctr$par0` ([`FitMCMC.R:143-150`](Package/R/FitMCMC.R:143)). But the claimed “trailing NA name” is not the clean source-level failure mode: if names survive, the sampler has an extra non-fixed dimension; if names are stripped, `f_rename_par()` is asked to assign a shorter name vector ([`Utils.R:370-374`](Package/R/Utils.R:370)).
+
+The AIC/BIC new test overclaims end-to-end coverage for the actual affected case. Its only real fit is `K=2` ([`test_Inference.R:138-143`](Package/tests/testthat/test_Inference.R:138)), where subtracting `1` and subtracting `K-1` are identical.
+
+REVIEW.md says all 11 new blocks fail on stock 2.51. Source-wise that is plausible for the five bug areas, but not literally demonstrated by every block: the AIC/BIC fitted subcase at `K=2` would pass under the old formula, and only the synthetic `K=3/4` part fails.
+tokens used
+97,458
+## Part 1 — claim verdicts
+
+| # | claim | verdict | evidence (file:line) |
+|---|---|---|---|
+| A1 | delta-method sandwich transposed; P-block anti-diagonal | CONFIRM | Original code uses `t(mJacob) %*% mInvHessian %*% mJacob` while `jacobian(f_mapPar, ...)` returns output-by-input derivatives ([`Inference.R:36-38`](Package/R/Inference.R:36)). The P block is anti-diagonal for `K=2`: `f_mapGamma()` fills off-diagonals by `mGamma[!mGamma]` then returns `c(t(mGamma[, -K]))`, keeping the original names ([`ParameterTransformation.R:204-209`](Package/R/ParameterTransformation.R:204)). |
+| A2 | p-value is one-sided under `Pr(>|t|)` | CONFIRM | The column is named `Pr(>|t|)` at construction, but the value is `1 - pnorm(abs(vTest))` ([`Inference.R:4-5`, `:43`](Package/R/Inference.R:4)). |
+| A3 | MCMC forecast uses linear index draw #1 | CONFIRM | `vol` is allocated as `(T+1) x nrow(par.check)` and line 32 indexes it as `vol[dim(PredProb)[1]]`, a single linear index ([`CondVol.R:19`, `:32`](Package/R/CondVol.R:19)). |
+| A4 | `dofMSGARCH` subtracts one instead of `K-1` per regime-constant parameter | CONFIRM | `dofMSGARCH` subtracts `length(regime.const.pars)` only, while `f_rename_par()` removes `name_2 ... name_K` for each constrained base name ([`Utils.R:420-421`, `:364-367`](Package/R/Utils.R:420)). |
+| B1 | stale-pointer recovery re-dereferences dead pointer; lines 349-353 unused | CONFIRM | The failed branch calls `spec$rcpp.func$get_mean()` and `get_sd()` immediately after `get_sd()` failed; computed `prior.mean/prior.sd` are not used by lines 355-356 ([`Utils.R:340-356`](Package/R/Utils.R:340)). |
+| B2 | raw `.i()` can abort under invalid K>=3 PLM transition rows; K=2 nonsingular | CONFIRM | `loadparam()` forms rows via `extract_P_it()` with last entry `1 - sum(P_it)`, then calls `(foo).i()` without a guard ([`MSgarch.h:193-198`, `:281-283`](Package/src/MSgarch.h:193)). For `K=2`, `det(I-P+U)=2(1+P_2_1-P_1_1)>0` when both free entries are in `(0,1)`. |
+| B3 | FitMCMC custom `par0` fixed-parameter path too long/mis-mapped | PARTIALLY CONFIRM | It is true FitMCMC substitutes fixed values but does not remove them before `f_unmapPar()` ([`FitMCMC.R:143-150`](Package/R/FitMCMC.R:143)), unlike FitML’s `f_remove_fixedpar()` ([`FitML.R:125-126`](Package/R/FitML.R:125)). But the exact “trailing NA name” mechanism is overstated: with names present the fixed dimension remains in the sampler, and with names absent `f_rename_par()` is more likely to fail length/name assignment than silently trail `NA` ([`Utils.R:370-374`](Package/R/Utils.R:370)). |
+| B4 | fixed transition probability accepted but breaks fitters by two mechanisms | CONFIRM | `CreateSpec()` validates against `out$label`, which includes transition labels ([`CreateSpec.R:278`, `:397-400`](Package/R/CreateSpec.R:397)); starting values pass recovered fixed names into single-regime `CreateSpec()` ([`StartingValues.R:201-206`](Package/R/StartingValues.R:201)). Separately, prior correction indexes `prior.mean[names(fixed.pars)]`, but MS priors come only from within-regime model params ([`Kernel.R:50-57`; `MSgarch.h:146-171`](Package/R/Kernel.R:50)). |
+| B5 | `llk == 1e+10` check wrong and try-error not handled | CONFIRM | `f_nll()` returns `+1e10` on likelihood failure, so `llk <- -optimizer$value` is `-1e10`; if `optim` itself fails, `f_OptimFUNDefault()` returns a `"try-error"` and `optimizer$value` is invalid ([`Likelihood.R:20-24`; `FitML.R:129-136`; `Optimizers.R:2-4`](Package/R/FitML.R:129)). |
+| C1 | MS `f_pdf`/`f_cdf` log path logs last-regime tmp, not mixture | CONFIRM | Both functions accumulate into `out`, then under `is_log` assign `out[i] = log(tmp[i])`, where `tmp` is overwritten in the regime loop ([`MSgarch.h:387-400`, `:462-475`](Package/src/MSgarch.h:387)). |
+| C2 | `f_cdf_its` first slice has swapped cube indices | CONFIRM | The initial CDF slice writes `tmp(ix, 0, s)` while later rows and `f_pdf_its` write `tmp(i, ix, s)` / `tmp(0, ix, s)` ([`MSgarch.h:490-507`; compare `:420-435`](Package/src/MSgarch.h:490)). |
+| C3 | `UncVol` window wrong and defaults swapped vs docs | CONFIRM | Docs say `nahead=5000`, `nburn=1000`, but `f_process_ctr(type=2)` sets `nburn=5000`, `nahead=1000`; `UncVol` then averages `tmp[ctr$nburn:ctr$nahead]` ([`UncVol.R:16-21`, `:88`; `Utils.R:159-160`](Package/R/UncVol.R:16)). |
+| C4 | mixture Viterbi transition matrix transposed/row sums wrong | CONFIRM | For mixtures, `TransMat()` returns a one-row probability vector; `State()` expands with `matrix(rep(P, K), nrow=K, ncol=K)`, which fills columns, so rows become repeated single probabilities rather than the mixture vector ([`State.R:90-94`; `Transmat.R:53-58`](Package/R/State.R:90)). |
+| C5 | K-expansion guard checks nonexistent `distribution.spec$model` | CONFIRM | The guard tests `length(distribution.spec$model) > 1`, but the checked object uses `$distribution` ([`CreateSpec.R:205`, `:212-219`](Package/R/CreateSpec.R:205)). |
+
+On the negative claims: I found no source-level contradiction to REVIEW.md’s claim that the base conditional densities/CDFs/quantiles, Fernandez-Steel moments, Hamilton filter, Kim smoother, valid-matrix ergodic formula, transition-memory layout, and tGARCH/gjrGARCH stationarity formulas are correct. Caveat: the ergodic formula is mathematically fine for valid stochastic matrices, but B2 shows the implementation does not guard invalid PLM matrices before inverting.
+
+## Part 2 — fixes
+
+A1: `mJacob %*% mInvHessian %*% t(mJacob)` is the correct orientation for `numDeriv::jacobian(f_mapPar, ...)`. The boundary adjustment changes the evaluation point to a nearby interior `vPw_mod`; it does not alter the orientation argument, and the Hessian is also computed at `vPw_mod` when recomputed. OK
+
+A2: `2 * (1 - pnorm(abs(vTest)))` matches the printed `Pr(>|t|)` label. The broader statistical caveat remains that zero-null z/t tests are not meaningful for all bounded/support-shifted parameters, but this fix corrects the advertised two-sided calculation. OK
+
+A3: `mean(vol[dim(PredProb)[1], ])` preserves ML behavior because a one-draw matrix row averages to the same scalar. Averaging posterior volatility rather than variance is a modeling/reporting choice already implied by `Volatility()`’s row-mean-on-volatility path ([`CondVol.R:53-55`](Package/R/CondVol.R:53)); it is separate from the indexing bug. OK
+
+A4: `- length(regime.const.pars) * (K - 1L)` is right for allowed specs: `K=1` rejects non-null regime constants, and `CreateSpec()` forbids using `fixed.pars` and `regime.const.pars` together ([`ParameterConstraints.R:115-121`; `CreateSpec.R:293-294`](Package/R/CreateSpec.R:293)). It also works for mixtures because the constrained duplicated within-regime parameters are still K copies collapsed to one. OK
+
+B1: Rebuilding only `rcpp.func` is sufficient from this source: all external-pointer-backed callables are collected there by `f_spec()` ([`CreateSpec.R:350-366`](Package/R/CreateSpec.R:350)), while `func` contains plain R closures. The `!is.null` guards should not affect normal `f_spec()` / `ExtractStateFit()` specs because `prior.mean` and `prior.sd` are initialized non-null ([`CreateSpec.R:412-416`](Package/R/CreateSpec.R:412)). OK
+
+## Part 3 — tests
+
+`test_Inference.R` “Standard errors use J V J'”: fails on original, has a non-vacuous guard, but partly restates implementation via `fit$Inference$Hessian`. OK
+
+`test_Inference.R` “observed information in natural scale”: independent and would fail on original; finite-difference tolerance and full `SMI` ML fit make it somewhat machine/runtime sensitive. WEAK
+
+`test_Inference.R` “Transition-probability standard errors”: fails on original and directly covers the P-block; uses fixed historical parameters, so it is reasonably stable. OK
+
+`test_Inference.R` “Pr(>|t|)”: fails on original and is not vacuous. OK
+
+`test_Inference.R` “Degrees of freedom”: catches original for `K=3/4`; independent enough. OK
+
+`test_Inference.R` “AIC and BIC”: the real fitted case is only `K=2`, where old and new formulas coincide; affected `K=3/4` cases are synthetic objects that test plumbing, not fitting. WEAK
+
+`test_Serialization.R` “round trip invalidates pointers”: useful anti-vacuity guard, but depends on current Rcpp serialization behavior. OK
+
+`test_Serialization.R` “reloaded spec usable”: fails on original and checks outputs against a live spec. OK
+
+`test_Serialization.R` “keeps priors”: fails on original and checks rebuilt C++ prior state. OK
+
+`test_Serialization.R` “reloaded ML fit usable”: fails on original; runtime acceptable but includes a full ML fit. OK
+
+`test_Serialization.R` “reloaded MCMC fit usable”: fails on original but is RNG/runtime-sensitive and uses a very short chain. WEAK
+
+`test_Volatility.R` “MCMC forecast averages draws”: should fail on original and has a non-vacuity guard; fragile because it runs MCMC in unit tests and depends on the first draw differing from the mean. WEAK
+
+## Missed or overstated
+
+B3 is overstated. The core bug is real: FitMCMC leaves the fixed parameter in user-supplied `ctr$par0` ([`FitMCMC.R:143-150`](Package/R/FitMCMC.R:143)). But the claimed “trailing NA name” is not the clean source-level failure mode: if names survive, the sampler has an extra non-fixed dimension; if names are stripped, `f_rename_par()` is asked to assign a shorter name vector ([`Utils.R:370-374`](Package/R/Utils.R:370)).
+
+The AIC/BIC new test overclaims end-to-end coverage for the actual affected case. Its only real fit is `K=2` ([`test_Inference.R:138-143`](Package/tests/testthat/test_Inference.R:138)), where subtracting `1` and subtracting `K-1` are identical.
+
+REVIEW.md says all 11 new blocks fail on stock 2.51. Source-wise that is plausible for the five bug areas, but not literally demonstrated by every block: the AIC/BIC fitted subcase at `K=2` would pass under the old formula, and only the synthetic `K=3/4` part fails.

From 1882d2c28a617378d5bc58480a52a004ffd2bc3b Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 19:07:36 -0400
Subject: [PATCH 08/20] Fix the transposed cube write in the multi-regime
 in-sample CDF

MSgarch::f_cdf_its allocates arma::cube tmp(ny, nx, K) and fills the t = 0
slice with tmp(ix, 0, s); every other assignment in the function, both
f_pdf_its implementations and SingleRegime::f_cdf_its use (0, ix, s). The
indices are reversed.

For a grid shorter than the sample the first observation's CDF was left at
zero for every evaluation point except the first, so PIT(..., do.its = TRUE)
with a user-supplied x returned a wrong first row. For a grid longer than the
sample the write ran off the cube; Armadillo's bounds checking, which
RcppArmadillo leaves enabled, turned that into
"Cube::operator(): index out of bounds" rather than a silent memory write, but
the call failed outright.

The new tests check the first row against a closed form -- at t = 1 each
regime's conditional variance is its unconditional variance and the predictive
state distribution is the ergodic one -- and cover the longer-grid case.
---
 Package/src/MSgarch.h                       |  2 +-
 Package/tests/testthat/test_NativeDensity.R | 41 +++++++++++++++++++++
 2 files changed, 42 insertions(+), 1 deletion(-)
 create mode 100644 Package/tests/testthat/test_NativeDensity.R

diff --git a/Package/src/MSgarch.h b/Package/src/MSgarch.h
index 22df1cb..67fdb25 100644
--- a/Package/src/MSgarch.h
+++ b/Package/src/MSgarch.h
@@ -494,7 +494,7 @@ inline arma::cube MSgarch::f_cdf_its(const NumericVector& theta,
   for (many::iterator it = specs.begin(); it != specs.end(); ++it) {
     sig = sqrt(vol[s].h);
     for (int ix = 0; ix < nx; ix++) {
-      tmp(ix, 0, s) = (*it)->spec_calc_cdf(x(ix, 0) / sig);  //
+      tmp(0, ix, s) = (*it)->spec_calc_cdf(x(ix, 0) / sig);  //
     }
     s++;
   }
diff --git a/Package/tests/testthat/test_NativeDensity.R b/Package/tests/testthat/test_NativeDensity.R
new file mode 100644
index 0000000..76ebe27
--- /dev/null
+++ b/Package/tests/testthat/test_NativeDensity.R
@@ -0,0 +1,41 @@
+testthat::context("Test the native mixture density and CDF entry points")
+
+data("SMI", package = "MSGARCH")
+spec <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                            distribution.spec = list(distribution = c("norm")),
+                            switch.spec = list(do.mix = FALSE, K = 2))
+par <- c(0.021631876185, 0.087024443479, 0.881493722371, 0.020659831566,
+         0.005396009353, 0.994040728662, 0.978348086740, 0.998703301894)
+y <- SMI[1:300]
+
+testthat::test_that("In-sample CDF fills the first time slice over the whole grid", {
+
+  # At t = 1 the conditional variance of each regime is its unconditional variance
+  # and the predictive state distribution is the ergodic one, so for a Normal
+  # specification the in-sample PIT has a closed form to check against.
+  x   <- c(-2, -1, 0, 1, 2)
+  h1  <- as.numeric(spec$rcpp.func$unc_vol_Rcpp(matrix(par, nrow = 1L)))
+  P0  <- State(object = spec, par = par, data = y)$PredProb[1L, 1L, ]
+  exp.pit <- vapply(x, function(z) sum(P0 * stats::pnorm(z / sqrt(h1))),
+                    FUN.VALUE = numeric(1))
+
+  est.pit <- as.numeric(PIT(object = spec, x = x, par = par, data = y,
+                            do.its = TRUE)[1L, ])
+
+  testthat::expect_true(max(abs(est.pit - exp.pit)) < 1e-12)
+
+  # the whole first row must be a genuine CDF, not a single value padded with zeros
+  testthat::expect_true(all(est.pit > 0 & est.pit < 1))
+  testthat::expect_true(!is.unsorted(est.pit, strictly = TRUE))
+
+})
+
+testthat::test_that("In-sample CDF accepts a grid longer than the sample", {
+
+  # the cube is (length(data), length(x), K); writing it transposed ran off the end
+  testthat::expect_silent(
+    PIT(object = spec, x = seq(from = -5, to = 5, length.out = 100L), par = par,
+        data = y[1:2], do.its = TRUE)
+  )
+
+})

From 828318ba02249cf864b293b863ca0ad955d0b02b Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 19:07:55 -0400
Subject: [PATCH 09/20] Return the mixture, not the last regime, from the
 native log branches

MSgarch::f_pdf and f_cdf accumulate the state-weighted mixture in out, but
under is_log they overwrote it with log(tmp[i]), where tmp holds only the
regime evaluated last. The reported log density and log CDF were therefore
those of regime K alone: on a two-regime fit to SMI, exp() of the returned
value differed from the mixture by factors of 12 and 26 in the tails.

Not reachable from PredPdf or PIT, which pass is_log = FALSE and take the
logarithm in R, but both routines are exposed as Rcpp module methods on
spec$rcpp.func. Taking log(out[i]) makes the native path agree with the R one
exactly, including its behaviour when the mixture underflows.
---
 Package/src/MSgarch.h                       |  4 ++--
 Package/tests/testthat/test_NativeDensity.R | 18 ++++++++++++++++++
 2 files changed, 20 insertions(+), 2 deletions(-)

diff --git a/Package/src/MSgarch.h b/Package/src/MSgarch.h
index 67fdb25..458573a 100644
--- a/Package/src/MSgarch.h
+++ b/Package/src/MSgarch.h
@@ -396,7 +396,7 @@ inline NumericVector MSgarch::f_pdf(const NumericVector& x,
   
   if (is_log) {
     for (int i = 0; i < nx; i++) {
-      out[i] = log(tmp[i]);
+      out[i] = log(out[i]);
     }
   }
   
@@ -471,7 +471,7 @@ inline NumericVector MSgarch::f_cdf(const NumericVector& x,
   
   if (is_log) {
     for (int i = 0; i < nx; i++) {
-      out[i] = log(tmp[i]);
+      out[i] = log(out[i]);
     }
   }
   
diff --git a/Package/tests/testthat/test_NativeDensity.R b/Package/tests/testthat/test_NativeDensity.R
index 76ebe27..61060c8 100644
--- a/Package/tests/testthat/test_NativeDensity.R
+++ b/Package/tests/testthat/test_NativeDensity.R
@@ -39,3 +39,21 @@ testthat::test_that("In-sample CDF accepts a grid longer than the sample", {
   )
 
 })
+
+testthat::test_that("The log branches return the mixture, not the last regime", {
+
+  x <- c(-3, -1, 0, 1, 3)
+
+  pdf.lin <- spec$rcpp.func$pdf_Rcpp(x, par, y, FALSE)
+  pdf.log <- spec$rcpp.func$pdf_Rcpp(x, par, y, TRUE)
+  testthat::expect_true(max(abs(exp(pdf.log) - pdf.lin)) < 1e-12)
+
+  cdf.lin <- spec$rcpp.func$cdf_Rcpp(x, par, y, FALSE)
+  cdf.log <- spec$rcpp.func$cdf_Rcpp(x, par, y, TRUE)
+  testthat::expect_true(max(abs(exp(cdf.log) - cdf.lin)) < 1e-12)
+
+  # guard: the regimes must actually differ here, otherwise "mixture" and
+  # "last regime" would coincide and this would prove nothing
+  testthat::expect_true(diff(range(spec$rcpp.func$unc_vol_Rcpp(matrix(par, nrow = 1L)))) > 1)
+
+})

From 6f04d2acd11827d2561b60f9e2e70a23cf97d568 Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 19:07:55 -0400
Subject: [PATCH 10/20] Evaluate the two-sided p-value on the lower tail

2 * (1 - pnorm(abs(t))) cancels to exactly zero once abs(t) exceeds about 8.3,
because pnorm() has already rounded to 1. Evaluating the lower tail directly
with 2 * pnorm(-abs(t)) stays accurate out to abs(t) of about 38 and is
identical wherever the old expression was representable.

Ordinary output is unchanged to within 9e-17. The single-regime GARCH fit used
in the tests has a t statistic of 26.9, whose p-value was reported as 0 and is
now 1.6e-159.
---
 Package/R/Inference.R                   | 2 +-
 Package/tests/testthat/test_Inference.R | 8 +++++++-
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/Package/R/Inference.R b/Package/R/Inference.R
index 1bea314..6971bad 100644
--- a/Package/R/Inference.R
+++ b/Package/R/Inference.R
@@ -40,7 +40,7 @@ f_InferenceFun <- function(vPw, data, spec, do.plm, mNegHessian = NULL) {
   vSE   <- sqrt(diag(mSandwitch))
   vTest <- vPn/vSE
   
-  vPvalues <- 2 * (1 - pnorm(abs(vTest)))
+  vPvalues <- 2 * pnorm(-abs(vTest))
   
   out[, "Estimate"]   <- vPn
   out[, "Std. Error"] <- vSE
diff --git a/Package/tests/testthat/test_Inference.R b/Package/tests/testthat/test_Inference.R
index 98624c6..d7e636f 100644
--- a/Package/tests/testthat/test_Inference.R
+++ b/Package/tests/testthat/test_Inference.R
@@ -97,9 +97,15 @@ testthat::test_that("Pr(>|t|) is a two-sided p-value", {
 
   mCoef <- fit.sr$Inference$MatCoef
   testthat::expect_true(max(abs(mCoef[, "Pr(>|t|)"] -
-                                2 * (1 - stats::pnorm(abs(mCoef[, "t value"]))))) < 1e-12)
+                                2 * stats::pnorm(-abs(mCoef[, "t value"])))) < 1e-12)
   testthat::expect_true(all(mCoef[, "Pr(>|t|)"] >= 0 & mCoef[, "Pr(>|t|)"] <= 1))
 
+  # 2 * (1 - pnorm(abs(t))) cancels to exactly zero once abs(t) exceeds about 8.3,
+  # while the lower tail stays representable out to about 38
+  vLive <- abs(mCoef[, "t value"]) < 37
+  testthat::expect_true(any(abs(mCoef[vLive, "t value"]) > 9))
+  testthat::expect_true(all(mCoef[vLive, "Pr(>|t|)"] > 0))
+
 })
 
 testthat::test_that("Degrees of freedom drop K - 1 values per regime-constant parameter", {

From 29f3abc519ff7ae5f99d2402097fb86b2bf5ca9d Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 19:07:55 -0400
Subject: [PATCH 11/20] Compare the estimated BIC with the expected one in the
 BIC test

The assertion read abs(exp.BIC - exp.BIC) < tol, which is zero by
construction, so the test passed regardless of what BIC() returned. The value
it was meant to check, 6841.1848542696416, is correct and the corrected
assertion passes.
---
 Package/tests/testthat/test_MLE.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Package/tests/testthat/test_MLE.R b/Package/tests/testthat/test_MLE.R
index 76a9dc7..88b80aa 100644
--- a/Package/tests/testthat/test_MLE.R
+++ b/Package/tests/testthat/test_MLE.R
@@ -24,7 +24,7 @@ testthat::test_that("Estimation BIC", {
   est.BIC <- BIC(fit)
   exp.BIC <- 6841.1848542696416
   
-  testthat::expect_true(abs(exp.BIC - exp.BIC) < tol)
+  testthat::expect_true(abs(est.BIC - exp.BIC) < tol)
   
 })
 

From 6bfff8bf566d35f4f247c825e3ab901025255252 Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 19:08:19 -0400
Subject: [PATCH 12/20] Remove the unreachable Sim methods for fit objects

R CMD check reported Sim.MSGARCH_ML_FIT and Sim.MSGARCH_MCMC_FIT as apparent
S3 methods that are not registered and whose signatures do not match the Sim
generic (newdata instead of data, no par).

Nothing can dispatch to them. Sim is not in NAMESPACE, so it is unreachable
from user code; every internal caller -- CondVol.R, Risk.R, Pit.R, PredPDF.R
-- passes a MSGARCH_SPEC, and simulate.MSGARCH_ML_FIT and
simulate.MSGARCH_MCMC_FIT pass object$spec explicitly. Removing them clears
the NOTE with no change in behaviour; simulate(), Risk(), PredPdf() and PIT()
with nahead > 1 are unaffected.
---
 Package/R/simulate.R | 16 ----------------
 1 file changed, 16 deletions(-)

diff --git a/Package/R/simulate.R b/Package/R/simulate.R
index 88cb305..3b9fd73 100644
--- a/Package/R/simulate.R
+++ b/Package/R/simulate.R
@@ -175,19 +175,3 @@ Sim.MSGARCH_SPEC <- function(object, data = NULL, nahead = 1L,
   class(out) <- "MSGARCH_SIM"
   return(out)
 }
-
-Sim.MSGARCH_ML_FIT <- function(object, newdata = NULL, nahead = 1L,
-                               nsim = 1L,  nburn = 500L, seed = NULL, ...) {
-  data <- c(object$data, newdata)
-  out  <- Sim(object = object$spec, data = data, nahead = nahead,
-              nsim = nsim, par = object$par, nburn = nburn, seed = seed)
-  return(out)
-}
-
-Sim.MSGARCH_MCMC_FIT <- function(object, newdata = NULL, nahead = 1L,
-                                 nsim = 1L, nburn = 500L, seed = NULL, ...) {
-  data <- c(object$data, newdata)
-  out  <- Sim(object = object$spec, data = data, nahead = nahead,
-              nsim = nsim, par = object$par, nburn = nburn, seed = seed)
-  return(out)
-}

From 9017d7b56acb5c96f9d110e54110c21ea1c3ad2b Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 19:08:19 -0400
Subject: [PATCH 13/20] Drop the obsolete Rcpp:::LdFlags() call from Makevars

Rcpp reports that linking against its library has not been needed since 2013
and that LdFlags() may be removed in 2027. LinkingTo: Rcpp, RcppArmadillo plus
the standard toolchain is sufficient; the remaining LAPACK, BLAS and FLIBS
entries are what RcppArmadillo actually requires.
---
 Package/src/Makevars     | 2 +-
 Package/src/Makevars.win | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/Package/src/Makevars b/Package/src/Makevars
index 75ade60..22ebc63 100644
--- a/Package/src/Makevars
+++ b/Package/src/Makevars
@@ -1 +1 @@
-PKG_LIBS=`$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
diff --git a/Package/src/Makevars.win b/Package/src/Makevars.win
index 5ccbaf0..22ebc63 100644
--- a/Package/src/Makevars.win
+++ b/Package/src/Makevars.win
@@ -1 +1 @@
-PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) `$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"`
+PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

From bcb849521b5abc85227d0847834814429def198d Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 19:08:19 -0400
Subject: [PATCH 14/20] Record the remaining 2.52 fixes in NEWS

---
 Package/NEWS | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/Package/NEWS b/Package/NEWS
index e04b454..3e6a0f4 100644
--- a/Package/NEWS
+++ b/Package/NEWS
@@ -4,6 +4,11 @@ Changes in Version 2.52
   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

From 6dffa2ab63b89fdbc264d2352fe19054d237797e Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 19:09:01 -0400
Subject: [PATCH 15/20] Update the review: C1 and C2 fixed, second-pass audit
 folded in

---
 REVIEW.md | 19 ++++++++++++++++---
 1 file changed, 16 insertions(+), 3 deletions(-)

diff --git a/REVIEW.md b/REVIEW.md
index deef0b7..5dfacf4 100644
--- a/REVIEW.md
+++ b/REVIEW.md
@@ -21,7 +21,20 @@ density, CDF, quantile function and truncated moment in `Normal.h` / `Student.h`
 > (new) and one block appended to `tests/testthat/test_Volatility.R`: 11 blocks / 34
 > assertions, all passing on the patched build and all 11 failing on stock 2.51.
 > `DESCRIPTION`'s version and `NEWS` are untouched — that is a release decision for the
-> maintainer. Everything from **B2** onwards is still open.
+> maintainer.
+>
+> **Second pass (2026-08-19).** A CRAN-readiness audit by Codex over the patched tree
+> re-found **C2** and **C1** and added four items this review had not covered: an obsolete
+> `Rcpp:::LdFlags()` in both `Makevars`, a vacuous assertion in the shipped BIC test
+> (`abs(exp.BIC - exp.BIC)`), cancellation in the two-sided p-value formula introduced by the
+> A2 fix above, and the observation that the two `Sim.*_FIT` methods behind the standing
+> `R CMD check` NOTE are unreachable dead code. All six are now **fixed**, with tests; the
+> package-origin NOTE is gone. Its "critical / undefined behaviour" framing of C2 did not
+> hold up — RcppArmadillo leaves bounds checking on, so the over-long-grid case raised
+> `Cube::operator(): index out of bounds` rather than writing out of bounds; the damaging
+> case was the silent one. **Still open: B2, B3, B4, B5, C3, C4, C5** and the input-validation
+> hardening (`f_check_y` accepts partial `NA`/`NaN`/`Inf`, after which `Volatility()` returns
+> a full plausible-looking series), which is a behaviour change and needs a separate decision.
 
 > **Independent cross-check (Codex).** The whole of this document — the 14 claims, the five
 > fixes and the tests — was re-audited read-only by `codex-cli 0.139.0` working from source
@@ -365,7 +378,7 @@ and make `f_check_y` reject `any(!is.finite(y))`.
 
 ## Tier 3 — silently wrong results in specific calls
 
-### C1. `pdf_Rcpp` / `cdf_Rcpp` with `is_log = TRUE` return the last regime, not the mixture — `src/MSgarch.h:397-401, 472-476`
+### C1. `pdf_Rcpp` / `cdf_Rcpp` with `is_log = TRUE` return the last regime, not the mixture — `src/MSgarch.h:397-401, 472-476` — **FIXED**
 
 ```cpp
 for (many::iterator it = specs.begin(); it != specs.end(); ++it) {
@@ -393,7 +406,7 @@ Not reachable from `PredPdf`/`PIT` (both always pass `FALSE` and take the log in
 are live methods on `spec$rcpp.func` and the flag is part of the C++ signature. Fix:
 accumulate into a scratch vector and take `log(out[i])`.
 
-### C2. `MSgarch::f_cdf_its` writes the first observation transposed — `src/MSgarch.h:497`
+### C2. `MSgarch::f_cdf_its` writes the first observation transposed — `src/MSgarch.h:497` — **FIXED**
 
 ```cpp
 tmp(ix, 0, s) = (*it)->spec_calc_cdf(x(ix, 0) / sig);   // t=0 block

From acf63b5e731db571a332f4363e9502a59d6f9fc5 Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 20:05:44 -0400
Subject: [PATCH 16/20] Make the native filter and the stationary solve
 numerically robust

Three changes in src/MSgarch.h.

The Hamilton filter and f_get_Pstate shifted each column of regime log
densities by its *smallest* entry before exponentiating, and only when that
entry fell below log(DBL_MIN). That guards against underflow but creates an
overflow: the largest exponent becomes max - min - 707, so once the regimes'
log densities differ by more than about 1400 the exponential returns Inf and
normalising it gives NaN. The gap is reachable with parameters the package
itself accepts -- a regime with alpha0 = 3.3e-4 against one with alpha0 = 1 is
enough -- after which the likelihood is floored to -1e10 and State(), PredPdf()
and Risk() return NaN or fail. Both filters now use the standard log-sum-exp
anchor, the largest entry, which cannot overflow and lets negligible regimes
underflow harmlessly. Verified against the analytic limit: as one regime's
density vanishes the likelihood converges to log P(other regime) plus that
regime's log density.

loadparam obtained the stationary distribution with a plain matrix inverse.
The plain parameter mapping bounds each free transition probability
separately, so with three or more regimes a row can leave the simplex and the
matrix is singular; the uncaught Armadillo exception then destroyed the whole
fit. FitMCMC always uses that mapping, so every chain with three or more
regimes was exposed, as was any ML fit with fixed or regime-constant
parameters. The solve is now guarded and falls back to the uniform
distribution, which calc_prior rejects a moment later anyway.

The unused log-likelihood accumulator in f_get_Pstate is removed; it produced
a compiler warning.
---
 Package/src/MSgarch.h | 49 +++++++++++++++++++++++++++----------------
 1 file changed, 31 insertions(+), 18 deletions(-)

diff --git a/Package/src/MSgarch.h b/Package/src/MSgarch.h
index 458573a..c799239 100644
--- a/Package/src/MSgarch.h
+++ b/Package/src/MSgarch.h
@@ -279,8 +279,19 @@ inline void MSgarch::loadparam(const NumericVector& theta) {
   arma::vec Uvec(K);
   Uvec.fill(1);
   arma::mat foo = (I - as<arma::mat>(P_mat) + Umat).t();
-  
-  arma::vec delta = (foo).i() * Uvec;
+
+  // Stationary distribution of the chain. 'theta' is not restricted to proper
+  // stochastic matrices on every code path (the plain parameter mapping bounds
+  // each free transition probability separately, so a row can leave the
+  // simplex), and 'foo' is then singular. Fall back to the uniform
+  // distribution instead of throwing: 'calc_prior' rejects such a parameter
+  // vector, so the likelihood computed from it is discarded anyway.
+  arma::vec delta;
+  bool solved = arma::solve(delta, foo, Uvec, arma::solve_opts::no_approx);
+  if (!solved || !delta.is_finite()) {
+    delta.set_size(K);
+    delta.fill(1.0 / K);
+  }
   for(int i = 0; i < K; i++){
     P0(i) = delta(i);
   }
@@ -635,18 +646,28 @@ inline NumericMatrix MSgarch::calc_lndMat(const NumericVector& y) {
 
 //-------------------------------------  Hamilton filter
 //-------------------------------------//
+// Shift applied to a column of regime log densities before exponentiating.
+// Anchoring on the largest entry is the usual log-sum-exp device: the largest
+// exponent becomes exactly zero, so the sum can never overflow, and entries far
+// below it underflow to zero harmlessly. Anchoring on the smallest entry
+// instead (as this code once did) overflows to Inf, and hence to NaN after
+// normalisation, whenever the regimes' log densities differ by more than about
+// 1400 -- reachable with admissible parameters when one regime is very tight.
+inline double lse_shift(const NumericVector& lndCol) {
+  double max_lnd = max(lndCol);
+  return (R_FINITE(max_lnd) ? -max_lnd : 0.0);
+}
+
 inline double MSgarch::HamiltonFilter(const NumericMatrix& lndMat) {
   int n_step = lndMat.ncol();
-  double lnd = 0, min_lnd, delta, sum_tmp;
+  double lnd = 0, delta, sum_tmp;
   NumericVector Pspot, Ppred, lndCol, tmp;
   
   // first step
   Pspot = clone(P0);             // Prob(St | I(t)
   Ppred = matrixProd(Pspot, P);  // one-step-ahead Prob(St | I(t-1))
   lndCol = lndMat(_, 0);
-  min_lnd = min(lndCol),
-    delta =
-      ((min_lnd < LND_MIN) ? LND_MIN - min_lnd : 0);  // handle over/under-flows
+  delta = lse_shift(lndCol);
   tmp = Ppred *
     exp(lndCol + delta);  // unormalized one-step-ahead Prob(St | I(t))
   
@@ -657,9 +678,7 @@ inline double MSgarch::HamiltonFilter(const NumericMatrix& lndMat) {
     Pspot = tmp / sum_tmp;
     Ppred = matrixProd(Pspot, P);
     lndCol = lndMat(_, t);
-    min_lnd = min(lndCol),
-      delta = ((min_lnd < LND_MIN) ? LND_MIN - min_lnd
-                 : 0);  // handle over/under-flows
+    delta = lse_shift(lndCol);
     tmp = Ppred * exp(lndCol + delta);
   }
   sum_tmp = sum(tmp);
@@ -678,7 +697,7 @@ inline List MSgarch::f_get_Pstate(const NumericVector& theta,
   NumericMatrix lndMat = calc_lndMat(y);  // likelihood in each state
   
   int n_step = lndMat.ncol();
-  double lnd = 0, min_lnd, delta, sum_tmp;
+  double delta, sum_tmp;
   NumericVector Pspot, Ppred, lndCol, tmp;
   arma::mat PtmpSpot(n_step + 1, K);
   arma::mat PtmpPred(n_step + 2, K);
@@ -698,16 +717,13 @@ inline List MSgarch::f_get_Pstate(const NumericVector& theta,
     PtmpPred(1, i) = Ppred(i);
   }
   lndCol = lndMat(_, 0);
-  min_lnd = min(lndCol),
-    delta =
-      ((min_lnd < LND_MIN) ? LND_MIN - min_lnd : 0);  // handle over/under-flows
+  delta = lse_shift(lndCol);
   tmp = Ppred *
     exp(lndCol + delta);  // unormalized one-step-ahead Prob(St | I(t))
   
   // remaining steps
   for (int t = 1; t < n_step; t++) {
     sum_tmp = sum(tmp);
-    lnd += -delta + log(sum_tmp);  // increment loglikelihood
     Pspot = tmp / sum_tmp;
     Ppred = matrixProd(Pspot, P);
     for (int i = 0; i < K; i++) {
@@ -717,13 +733,10 @@ inline List MSgarch::f_get_Pstate(const NumericVector& theta,
       PtmpPred(t + 1, i) = Ppred(i);
     }
     lndCol = lndMat(_, t);
-    min_lnd = min(lndCol),
-      delta = ((min_lnd < LND_MIN) ? LND_MIN - min_lnd
-                 : 0);  // handle over/under-flows
+    delta = lse_shift(lndCol);
     tmp = Ppred * exp(lndCol + delta);
   }
   sum_tmp = sum(tmp);
-  lnd += -delta + log(sum_tmp);  // increment loglikelihood
   Pspot = tmp / sum_tmp;
   PLast = matrixProd(Pspot, P);
   

From 7fcb09d7e24911cf44556ef71fc6af9d4d6be8b7 Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 20:06:15 -0400
Subject: [PATCH 17/20] Fix time-index handling and validate public inputs

Time index. The fitted-object methods of Volatility, predict, PIT, PredPdf and
Risk built their combined series by concatenating object$data with newdata and
then deriving the index from the already-concatenated result before appending
another length(newdata) points. The index was longer than the values, so
zooreg() recycled observations and the model conditioned on values that are
not in the sample: a 200-point monthly series plus two new returns produced 204
observations ending in two values copied from the start. The same code also
discarded the original start and frequency even when newdata was NULL, turning
a monthly series into an annual one, and forecast indexes advanced by one index
unit rather than by the series' own step. Three helpers in Utils.R --
f_combine_data, f_future_index and f_index_result -- now do this in one place,
and ts input gives numerically identical results to the equivalent numeric
input.

Input validation. f_check_y and f_check_par only rejected input that was
entirely NaN, so a single NA, NaN or Inf reached the compiled code; Volatility()
would then return a complete, plausible-looking series computed from corrupt
data. Both now require finite values, and the data must hold at least two
observations, since the first initialises the variance recursion. f_nll and
f_posterior treat a non-finite mapped parameter vector as an infeasible point
so that a strict check cannot abort an optimization.

Failure reporting. FitML tested the optimizer's result against +1e10, but f_nll
returns +1e10 so a failed optimization arrives as -1e10; the optimizer is also
wrapped in try(), so optimizer$value was reached first. Every failure therefore
surfaced as "$ operator is invalid for atomic vectors". Failures are now
reported for what they are.

Risk. alpha, nahead and ctr$nmesh are validated. Because the evaluation grid
spans the observed data range only, it can miss part of the predictive
distribution; when the mass it omits exceeds the requested tail probability the
quantile read off it is the grid boundary rather than a quantile, and that now
warns.

Constraints. Transition probabilities were accepted by constraint.spec$fixed
but broke both fitters, in the starting-value routine and again in the prior
correction, so they are refused with an explanation. FitMCMC now drops fixed
parameters from a user-supplied ctr$par0, as FitML already did; leaving them in
shifted every later parameter onto the wrong bounds. The identification sort is
skipped when parameters are fixed, since relabelling regimes by unconditional
variance moves the fixed value into a different regime -- with the default
do.sort a parameter fixed at 0.8 came back taking sixteen different values.

Other. simulate() accepts nburn = 0, which used to drop the first draw because
1:0 is c(1, 0) and then fail on the dimnames. UncVol averages the horizons
after the burn-in instead of nburn:nahead, a descending range under the shipped
defaults. prior$sd is validated by its own checker and must be finite and
strictly positive. CreateSpec validates switch.spec$K and, when expanding one
regime through K, rejects an explicitly heterogeneous distribution vector; the
guard used to test distribution.spec$model, which does not exist. The mixture
transition matrix passed to Viterbi decoding is built row-wise so that it is
row-stochastic; the decoded path is unchanged, because the misplaced factor is
constant in the index being maximised over.
---
 Package/R/CondVol.R              | 16 +-----
 Package/R/CreateSpec.R           | 20 +++++--
 Package/R/FitMCMC.R              | 15 +++++-
 Package/R/FitML.R                | 24 ++++++---
 Package/R/Likelihood.R           |  6 +++
 Package/R/ParameterConstraints.R | 24 +++++++--
 Package/R/Pit.R                  | 38 ++------------
 Package/R/Posterior.R            |  4 ++
 Package/R/PredPDF.R              | 38 ++------------
 Package/R/Risk.R                 | 90 +++++++++++++++-----------------
 Package/R/State.R                |  2 +-
 Package/R/UncVol.R               |  5 +-
 Package/R/Utils.R                | 88 ++++++++++++++++++++++++++++---
 Package/R/Volatility.R           | 20 +------
 Package/R/predict.R              | 20 +------
 Package/R/simulate.R             | 17 +++---
 16 files changed, 227 insertions(+), 200 deletions(-)

diff --git a/Package/R/CondVol.R b/Package/R/CondVol.R
index ccdd60e..ac0ac82 100644
--- a/Package/R/CondVol.R
+++ b/Package/R/CondVol.R
@@ -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))
diff --git a/Package/R/CreateSpec.R b/Package/R/CreateSpec.R
index 4b4df44..dba3df7 100644
--- a/Package/R/CreateSpec.R
+++ b/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)
   }
diff --git a/Package/R/FitMCMC.R b/Package/R/FitMCMC.R
index ad2e4b6..d98c7d5 100644
--- a/Package/R/FitMCMC.R
+++ b/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
diff --git a/Package/R/FitML.R b/Package/R/FitML.R
index 0715862..dce6383 100644
--- a/Package/R/FitML.R
+++ b/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)
diff --git a/Package/R/Likelihood.R b/Package/R/Likelihood.R
index 545e45a..12648d8 100644
--- a/Package/R/Likelihood.R
+++ b/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)) {
diff --git a/Package/R/ParameterConstraints.R b/Package/R/ParameterConstraints.R
index a5b869f..31d794e 100644
--- a/Package/R/ParameterConstraints.R
+++ b/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) {
diff --git a/Package/R/Pit.R b/Package/R/Pit.R
index b166532..488b20c 100644
--- a/Package/R/Pit.R
+++ b/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)
diff --git a/Package/R/Posterior.R b/Package/R/Posterior.R
index 392b48a..bfa7bb9 100644
--- a/Package/R/Posterior.R
+++ b/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)) {
diff --git a/Package/R/PredPDF.R b/Package/R/PredPDF.R
index eab1ea2..86ad6f4 100644
--- a/Package/R/PredPDF.R
+++ b/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)
diff --git a/Package/R/Risk.R b/Package/R/Risk.R
index 1397a15..bf44134 100644
--- a/Package/R/Risk.R
+++ b/Package/R/Risk.R
@@ -106,36 +106,58 @@ Risk.MSGARCH_SPEC <- function(object, par, data, alpha = c(0.01, 0.05), nahead =
   object  <- f_check_spec(object)
   data_    <- f_check_y(data)
   ctr     <- f_process_ctr(ctr)
+
+  if (!is.numeric(alpha) || length(alpha) < 1L || any(!is.finite(alpha)) ||
+      any(alpha <= 0) || any(alpha >= 1)) {
+    stop("alpha must contain finite probabilities strictly between 0 and 1.")
+  }
+  if (!is.numeric(nahead) || length(nahead) != 1L || !is.finite(nahead) ||
+      nahead < 1 || nahead != round(nahead)) {
+    stop("nahead must be a single positive whole number.")
+  }
+  if (!is.numeric(ctr$nmesh) || length(ctr$nmesh) != 1L || !is.finite(ctr$nmesh) ||
+      ctr$nmesh < 2 || ctr$nmesh != round(ctr$nmesh)) {
+    stop("ctr$nmesh must be a single whole number greater than or equal to two.")
+  }
+
   out     <- list()
   n.alpha <- length(alpha)
-  xmin    <- min(data_) - sd(data_)
-  xmax    <- max(data_) + sd(data_)
+  dSd     <- sd(data_)
+  if (!is.finite(dSd) || dSd <= 0) {
+    stop("the data have zero or undefined dispersion, so no evaluation grid ",
+         "can be built for the predictive distribution.")
+  }
+  xmin    <- min(data_) - dSd
+  xmax    <- max(data_) + dSd
   
   x     <- seq(from = xmin, to = xmax, length.out = ctr$nmesh)
   pdf_x <- PredPdf(object = object, par = par, x = x, data = data_, do.its = do.its, log = FALSE)
   cumul <- apply(pdf_x, 1L, cumsum) * (x[2L] - x[1L])
+
+  # The grid spans the observed range only, so it can miss part of the
+  # predictive distribution. Whatever mass falls outside it is unaccounted for,
+  # and if that exceeds the requested tail probability the quantile read off the
+  # grid is meaningless -- it is pinned to an endpoint rather than solved for.
+  dMissing <- 1 - min(cumul[nrow(cumul), ])
+  if (is.finite(dMissing) && dMissing > min(alpha)) {
+    warning("The evaluation grid covers only ",
+            format(100 * (1 - dMissing), digits = 4),
+            "% of the predictive distribution, which is less than the requested ",
+            "tail probability alpha = ", format(min(alpha), digits = 3),
+            ". VaR and ES are pinned to the grid boundary and should not be used. ",
+            "Use a longer sample, or widen the grid via ctr$nmesh and the data range.",
+            call. = FALSE)
+  }
   out   <- list()
   draw  <- NULL
   if (do.its == TRUE) {
     out$VaR <- matrix(NA, nrow = nrow(pdf_x), ncol = n.alpha)
     rownames(out$VaR) <-  paste0("t=",1:length(data_))
-    if(zoo::is.zoo(data)){
-      out$VaR = zoo::zooreg(out$VaR, order.by = zoo::index(data))
-    }
-    if(is.ts(data)){
-      out$VaR = zoo::zooreg(out$VaR, order.by = zoo::index(data))
-      out$VaR = as.ts(out$VaR)
-    }
+    out$VaR <- f_index_result(out$VaR, data)
   } else {
     out$VaR <- matrix(NA, nrow = nahead, ncol = n.alpha)
     rownames(out$VaR) <-  paste0("h=",1:nahead)
-    if(zoo::is.zoo(data)){
-      out$VaR = zoo::zooreg(out$VaR, order.by = zoo::index(data)[length(data)]+(1:nahead))
-    }
-    if(is.ts(data)){
-      out$VaR = zoo::zooreg(out$VaR, order.by = zoo::index(data)[length(data)]+(1:nahead))
-      out$VaR = as.ts(out$VaR)
-    }
+    out$VaR <- f_index_result(out$VaR, data, nahead)
   }
   for (n in 1:nrow(pdf_x)) {
     for (i in 1:n.alpha) {
@@ -157,23 +179,11 @@ Risk.MSGARCH_SPEC <- function(object, par, data, alpha = c(0.01, 0.05), nahead =
     if (do.its == TRUE) {
       out$ES <- matrix(NA, nrow = nrow(pdf_x), ncol = n.alpha)
       rownames(out$ES) <- paste0("t=",1:length(data_))
-      if(zoo::is.zoo(data)){
-        out$ES = zoo::zooreg(out$ES, order.by = zoo::index(data))
-      }
-      if(is.ts(data)){
-        out$ES = zoo::zooreg(out$ES, order.by = zoo::index(data))
-        out$ES = as.ts(out$ES)
-      }
+      out$ES <- f_index_result(out$ES, data)
     } else {
       out$ES <- matrix(NA, nrow = nahead, ncol = n.alpha)
       rownames(out$ES) <- paste0("h=",1:nahead)
-      if(zoo::is.zoo(data)){
-        out$ES = zoo::zooreg(out$ES, order.by = zoo::index(data)[length(data)]+(1:nahead))
-      }
-      if(is.ts(data)){
-        out$ES = zoo::zooreg(out$ES, order.by = zoo::index(data)[length(data)]+(1:nahead))
-        out$ES = as.ts(out$ES)
-      }
+      out$ES <- f_index_result(out$ES, data, nahead)
     }
     for (n in 1:nrow(pdf_x)) {
       for (i in 1:n.alpha) {
@@ -201,15 +211,7 @@ Risk.MSGARCH_SPEC <- function(object, par, data, alpha = c(0.01, 0.05), nahead =
 #' @export
 Risk.MSGARCH_ML_FIT <- function(object, newdata = NULL, alpha = c(0.01, 0.05),
                                 do.es = 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  <- Risk(object = object$spec, par = object$par, data = data, alpha = alpha,
                do.es = do.es, do.its = do.its, nahead = nahead, do.cumulative = do.cumulative, ctr = ctr)
   return(out)
@@ -219,15 +221,7 @@ Risk.MSGARCH_ML_FIT <- function(object, newdata = NULL, alpha = c(0.01, 0.05),
 #' @export
 Risk.MSGARCH_MCMC_FIT <- function(object, newdata = NULL, alpha = c(0.01, 0.05),
                                   do.es = 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  <- Risk(object = object$spec, par = object$par, data = data, alpha = alpha,
                do.es = do.es, do.its = do.its, nahead = nahead, do.cumulative = do.cumulative, ctr = ctr)
   return(out)
diff --git a/Package/R/State.R b/Package/R/State.R
index be32958..6b2ba54 100644
--- a/Package/R/State.R
+++ b/Package/R/State.R
@@ -89,7 +89,7 @@ State.MSGARCH_SPEC <- function(object, par, data, ...) {
     if (object$K > 1) {
       P <- TransMat(object = object, par = par[i, ], nahead = 1)
       if (isTRUE(object$is.mix)) {
-        P <- matrix(rep(P, object$K), nrow = object$K, ncol = object$K)
+        P <- matrix(rep(P, object$K), nrow = object$K, ncol = object$K, byrow = TRUE)
       }
       out$Viterbi[2:length(data), i] <- Viterbi(tmp$LL, P, object$K)
     } else {
diff --git a/Package/R/UncVol.R b/Package/R/UncVol.R
index c8f4ee1..6651a1e 100644
--- a/Package/R/UncVol.R
+++ b/Package/R/UncVol.R
@@ -85,7 +85,10 @@ UncVol.MSGARCH_SPEC <- function(object, par = NULL, ctr = list(), ...) {
                  do.its = FALSE,
                  nahead = ctr$nburn + ctr$nahead,
                  ctr = list(nsim = nsim))$vol
-  out <- mean(tmp[ctr$nburn:ctr$nahead])
+  # discard the first nburn horizons and keep the following nahead ones; the
+  # old expression, nburn:nahead, is a descending sequence under the shipped
+  # defaults and in general keeps part of the transient
+  out <- mean(tmp[seq.int(from = ctr$nburn + 1L, length.out = ctr$nahead)])
   return(out)
 }
 
diff --git a/Package/R/Utils.R b/Package/R/Utils.R
index fd7bc92..af3d7b4 100644
--- a/Package/R/Utils.R
+++ b/Package/R/Utils.R
@@ -174,8 +174,13 @@ f_check_y <- function(y) {
   if (!is.numeric(y)) {
     stop("y must be numeric")
   }
-  if (all(is.nan(y))) {
-    stop("nan dectected in y")
+  if (length(y) < 2L) {
+    stop("the data must contain at least two observations: the first one ",
+         "initialises the variance recursion and the likelihood is evaluated ",
+         "on the remainder.")
+  }
+  if (anyNA(y) || any(!is.finite(y))) {
+    stop("the data must not contain NA, NaN or infinite values.")
   }
   return(y)
 }
@@ -188,8 +193,8 @@ f_check_par <- function(spec, par) {
   if (!is.numeric(par)) {
     stop("par must be a numeric")
   }
-  if (all(is.nan(par))) {
-    stop("nan dectected in par")
+  if (anyNA(par) || any(!is.finite(par))) {
+    stop("par must not contain NA, NaN or infinite values.")
   }
   len.par <- length(spec$par0)
   if (is.vector(par)) {
@@ -337,6 +342,69 @@ f_match <- function(x, target) {
   return(toupper(substr(x, 1, 3)) == target)
 }
 
+# Combine a fit's data with newdata, keeping the time index of the original
+# series. Index and values must stay the same length: deriving the index from
+# the already-concatenated series and then appending length(newdata) further
+# points makes zooreg() recycle observations, so the model would be conditioned
+# on values that are not in the sample.
+f_combine_data <- function(data, newdata = NULL) {
+  if (is.null(newdata) || length(newdata) == 0L) {
+    return(data)
+  }
+  bIsTs  <- is.ts(data)
+  bIsZoo <- zoo::is.zoo(data)
+  vValues <- c(as.numeric(data), as.numeric(newdata))
+  if (!bIsTs && !bIsZoo) {
+    return(vValues)
+  }
+  if (bIsTs) {
+    # keep the original start and frequency; only the length changes
+    return(stats::ts(vValues, start = stats::start(data),
+                     frequency = stats::frequency(data)))
+  }
+  vIndex <- zoo::index(data)
+  n      <- length(vIndex)
+  dStep  <- if (n > 1L) vIndex[n] - vIndex[n - 1L] else 1
+  vIndex <- c(vIndex, vIndex[n] + dStep * seq_len(length(newdata)))
+  return(zoo::zoo(vValues, order.by = vIndex))
+}
+
+# Index for nahead points beyond the end of the series, advancing by the
+# series' own spacing. Adding 1:nahead to the last index value instead moves a
+# monthly series forward by whole years.
+f_future_index <- function(data, nahead) {
+  vIndex <- zoo::index(data)
+  n      <- length(vIndex)
+  dStep  <- if (is.ts(data)) {
+    1 / stats::frequency(data)
+  } else if (n > 1L) {
+    vIndex[n] - vIndex[n - 1L]
+  } else {
+    1
+  }
+  return(vIndex[n] + dStep * seq_len(nahead))
+}
+
+# Attach the time index of `data` to a result. `nahead = NULL` means the result
+# is in-sample and shares the index of `data`; otherwise it covers nahead points
+# beyond the end of the sample. For a ts the series is rebuilt from start and
+# frequency rather than round-tripped through zoo, whose index-regularity test
+# fails on the floating-point spacing of, say, a monthly series.
+f_index_result <- function(x, data, nahead = NULL) {
+  if (is.ts(data)) {
+    dFreq <- stats::frequency(data)
+    if (is.null(nahead)) {
+      return(stats::ts(x, start = stats::start(data), frequency = dFreq))
+    }
+    return(stats::ts(x, start = stats::tsp(data)[2L] + 1 / dFreq, frequency = dFreq))
+  }
+  if (zoo::is.zoo(data)) {
+    vIndex <- if (is.null(nahead)) zoo::index(data) else f_future_index(data, nahead)
+    return(zoo::zoo(x, order.by = vIndex))
+  }
+  return(x)
+}
+
 f_check_spec <- function(spec) {
   is.OK = tryCatch({
     spec$rcpp.func$get_sd()
@@ -406,7 +474,11 @@ f_rbindrep = function(mat, n) {
 f_check_parameterPriorMean <- function(prior.mean, vParNames) {
   if (any(!names(prior.mean) %in% vParNames)) {
     vWrongPars <- names(prior.mean)[!names(prior.mean) %in% vParNames]
-    stop(cat(paste("Wrong name in prior.mean:", vWrongPars)))
+    stop("Wrong name in prior.mean: ", paste(vWrongPars, collapse = ", "))
+  }
+  vMean <- unlist(prior.mean)
+  if (length(vMean) > 0L && (!is.numeric(vMean) || any(!is.finite(vMean)))) {
+    stop("Every entry of prior$mean must be a finite number.")
   }
   return(prior.mean)
 }
@@ -414,7 +486,11 @@ f_check_parameterPriorMean <- function(prior.mean, vParNames) {
 f_check_parameterPriorSd <- function(prior.sd, vParNames) {
   if (any(!names(prior.sd) %in% vParNames)) {
     vWrongPars <- names(prior.sd)[!names(prior.sd) %in% vParNames]
-    stop(cat(paste("Wrong name in prior.sd:", vWrongPars)))
+    stop("Wrong name in prior.sd: ", paste(vWrongPars, collapse = ", "))
+  }
+  vSd <- unlist(prior.sd)
+  if (length(vSd) > 0L && (!is.numeric(vSd) || any(!is.finite(vSd)) || any(vSd <= 0))) {
+    stop("Every entry of prior$sd must be a finite number strictly greater than zero.")
   }
   return(prior.sd)
 }
diff --git a/Package/R/Volatility.R b/Package/R/Volatility.R
index f0940a9..112add8 100644
--- a/Package/R/Volatility.R
+++ b/Package/R/Volatility.R
@@ -55,15 +55,7 @@ Volatility.MSGARCH_SPEC <- function(object, par, data, ...) {
 #' @rdname Volatility
 #' @export
 Volatility.MSGARCH_ML_FIT <- function(object, newdata = NULL, ...) {
-  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  <- f_CondVol(object = object$spec, par = object$par, data = data,
                   do.its = TRUE, ctr = list())
   return(out$vol)
@@ -72,15 +64,7 @@ Volatility.MSGARCH_ML_FIT <- function(object, newdata = NULL, ...) {
 #' @rdname Volatility
 #' @export
 Volatility.MSGARCH_MCMC_FIT <- function(object, newdata = NULL, ...) {
-  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  <- f_CondVol(object = object$spec, par = object$par, data = data,
                   do.its = TRUE, ctr = list())
   return(out$vol)
diff --git a/Package/R/predict.R b/Package/R/predict.R
index 3bfee91..5e25ecf 100644
--- a/Package/R/predict.R
+++ b/Package/R/predict.R
@@ -81,15 +81,7 @@ predict.MSGARCH_SPEC <- function(object, newdata = NULL, nahead = 1L,
 predict.MSGARCH_ML_FIT <- function(object, newdata = NULL, 
                                    nahead = 1L, do.return.draw = FALSE, 
                                    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  <- f_CondVol(object = object$spec, par = object$par, data = data, nahead = nahead,
                     do.its = FALSE, do.cumulative = do.cumulative, ctr = ctr)
   if(!isTRUE(do.return.draw)){
@@ -103,15 +95,7 @@ predict.MSGARCH_ML_FIT <- function(object, newdata = NULL,
 #' @export
 predict.MSGARCH_MCMC_FIT <- function(object, newdata = NULL, nahead = 1L, 
                                      do.return.draw = FALSE, 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  <- f_CondVol(object = object$spec, par = object$par, data = data, nahead = nahead,
                     do.its = FALSE, do.cumulative = do.cumulative, ctr = ctr)
   if(!isTRUE(do.return.draw)){
diff --git a/Package/R/simulate.R b/Package/R/simulate.R
index 3b9fd73..2c4d056 100644
--- a/Package/R/simulate.R
+++ b/Package/R/simulate.R
@@ -126,9 +126,12 @@ Sim.MSGARCH_SPEC <- function(object, data = NULL, nahead = 1L,
       start <- start + nsim
       end   <- end + nsim
     }
-    draw  <- draw[-(1:nburn),,drop = FALSE]
-    state <- state[-(1:nburn),,drop = FALSE]
-    CondVol <- CondVol[-(1:nburn),,,drop = FALSE]
+    if (nburn > 0L) {
+      # 1:0 is c(1, 0), so the plain negative index would drop the first draw
+      draw  <- draw[-seq_len(nburn), , drop = FALSE]
+      state <- state[-seq_len(nburn), , drop = FALSE]
+      CondVol <- CondVol[-seq_len(nburn), , , drop = FALSE]
+    }
     rownames(draw) = rownames(state) = paste0("t=",1:nahead)
     colnames(draw) = colnames(state) =  paste0("Sim #",1:(nsim * nrow(par)))
     dimnames(CondVol)[[1]] = paste0("t=",1:nahead)
@@ -160,13 +163,7 @@ Sim.MSGARCH_SPEC <- function(object, data = NULL, nahead = 1L,
     }
     rownames(draw) = rownames(state) = paste0("h=",1:nahead)
     colnames(draw) = colnames(state) =  paste0("Sim #",1:(nsim * nrow(par)))
-    if(zoo::is.zoo(data)){
-      draw = zoo::zooreg(draw, order.by =  zoo::index(data)[length(data)]+(1:nahead))
-    }
-    if(is.ts(data)){
-      draw = zoo::zooreg(draw, order.by =  zoo::index(data)[length(data)]+(1:nahead))
-      draw = as.ts(draw)
-    }
+    draw <- f_index_result(draw, data, nahead)
   }
   out <- list()
   out$draw <- draw

From 4304f8f3af49a167fe1d5d9bd3ca2d7326eef0a4 Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 20:06:15 -0400
Subject: [PATCH 18/20] Add regression tests for the robustness and validation
 fixes

test_Robustness.R covers the filter under a between-regime log-density gap of
more than 1000, including a check against the analytic limit; a singular
transition matrix, both directly and through a constrained three-regime fit;
non-finite and too-short data and parameters; nburn = 0; the UncVol averaging
window; the mixture Viterbi path against an independent per-observation MAP
decoding; CreateSpec's K validation; Risk's argument checks and its
grid-coverage warning; prior mean and standard deviation validation; and that a
fixed parameter survives FitML, a user-supplied par0 and the identification
sort.

test_TimeSeriesIndex.R checks that combining a fit's data with newdata keeps
index and values aligned, that ts and zoo input give the same numbers as the
equivalent numeric input through all five wrappers, and that forecast indexes
advance by the series' own time step.

Twelve of the fourteen blocks fail on the previous commit. The two that do not
are the mixture Viterbi path, which was never wrong, and the filter block,
whose failure the testthat summary does not surface although the likelihood it
asserts on is -1e10 there.
---
 Package/tests/testthat/test_Robustness.R      | 207 ++++++++++++++++++
 Package/tests/testthat/test_TimeSeriesIndex.R |  76 +++++++
 2 files changed, 283 insertions(+)
 create mode 100644 Package/tests/testthat/test_Robustness.R
 create mode 100644 Package/tests/testthat/test_TimeSeriesIndex.R

diff --git a/Package/tests/testthat/test_Robustness.R b/Package/tests/testthat/test_Robustness.R
new file mode 100644
index 0000000..afa4134
--- /dev/null
+++ b/Package/tests/testthat/test_Robustness.R
@@ -0,0 +1,207 @@
+testthat::context("Test numerical robustness and input validation")
+
+data("SMI", package = "MSGARCH")
+spec <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH", "sGARCH")),
+                            distribution.spec = list(distribution = c("norm", "norm")))
+
+testthat::test_that("The filter survives extreme between-regime likelihood ratios", {
+
+  # One tight regime and one loose one: at y = 1 their log densities differ by
+  # far more than the ~1400 that used to overflow the exponential and turn the
+  # filtered probabilities into NaN.
+  y   <- c(0, 1)
+  par <- c(1e-5, 1e-6, 0.001, 1, 1e-6, 0.001, 0.9, 0.1)
+  h   <- c(1e-5, 1) / (1 - 1e-6 - 0.001)
+  testthat::expect_true(abs(diff(stats::dnorm(1, 0, sqrt(h), log = TRUE))) > 1000)
+
+  dLLK <- MSGARCH:::Kernel(spec, par, y, log = TRUE, do.prior = FALSE)
+  testthat::expect_true(is.finite(dLLK))
+
+  # regime 1 is impossible at this observation, so the likelihood collapses onto
+  # regime 2 weighted by its predictive probability (here the ergodic 1/2)
+  testthat::expect_true(abs(dLLK - (stats::dnorm(1, 0, sqrt(h[2]), log = TRUE) + log(0.5))) < 1e-8)
+
+  mProb <- State(object = spec, par = par, data = y)$FiltProb
+  testthat::expect_true(all(is.finite(mProb)))
+  testthat::expect_true(max(abs(apply(mProb, 1L, sum) - 1)) < 1e-12)
+
+  testthat::expect_true(is.finite(as.numeric(PredPdf(object = spec, par = par, x = 0, data = y))))
+
+})
+
+testthat::test_that("A singular transition matrix is rejected, not thrown from C++", {
+
+  # the plain mapping bounds each free transition probability separately, so a
+  # row can leave the simplex and make the stationary-distribution solve singular
+  spec3 <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                               switch.spec = list(do.mix = FALSE, K = 3))
+  par <- spec3$par0
+  par[grep("^P_", names(par))] <- c(1, 0.999973, 0, 1, 0, 0)
+
+  testthat::expect_true(is.finite(MSGARCH:::Kernel(spec3, par, SMI[1:200], log = TRUE,
+                                                   do.prior = FALSE)))
+  testthat::expect_equal(MSGARCH:::Kernel(spec3, par, SMI[1:200], log = TRUE,
+                                          do.prior = TRUE), -1e10)
+
+  # a constrained three-regime fit used to die on such a point mid-optimization
+  spec3c <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                                distribution.spec = list(distribution = c("std")),
+                                switch.spec = list(do.mix = FALSE, K = 3),
+                                constraint.spec = list(regime.const = c("nu")))
+  set.seed(1234)
+  fit <- MSGARCH::FitML(spec3c, data = SMI, ctr = list(do.se = FALSE))
+  testthat::expect_true(is.finite(fit$loglik))
+
+})
+
+testthat::test_that("Data and parameters must be finite and long enough", {
+
+  y <- SMI[1:100]
+  for (bad in list(replace(y, 10L, NA), replace(y, 10L, NaN), replace(y, 10L, Inf))) {
+    testthat::expect_error(Volatility(object = spec, par = spec$par0, data = bad))
+    testthat::expect_error(MSGARCH::FitML(spec, data = bad))
+  }
+  testthat::expect_error(Volatility(object = spec, par = spec$par0, data = y[1L]))
+  testthat::expect_error(Volatility(object = spec, par = replace(spec$par0, 1L, NA), data = y))
+
+  # a well-formed call must still work
+  testthat::expect_true(all(is.finite(as.numeric(Volatility(object = spec, par = spec$par0,
+                                                            data = y)))))
+
+})
+
+testthat::test_that("simulate() accepts a zero burn-in", {
+
+  set.seed(1234)
+  sim0 <- simulate(object = spec, nsim = 2L, nahead = 3L, par = spec$par0, nburn = 0L)
+  testthat::expect_equal(dim(sim0$draw), c(3L, 2L))
+
+  set.seed(1234)
+  sim5 <- simulate(object = spec, nsim = 2L, nahead = 3L, par = spec$par0, nburn = 5L)
+  testthat::expect_equal(dim(sim5$draw), c(3L, 2L))
+
+})
+
+testthat::test_that("UncVol averages the horizons after the burn-in", {
+
+  set.seed(1234)
+  fit <- MSGARCH::FitML(spec, data = SMI, ctr = list(do.se = FALSE))
+  ctr <- list(nsim = 200L, nburn = 20L, nahead = 5L)
+
+  set.seed(7)
+  est <- UncVol(object = fit, ctr = ctr)
+  set.seed(7)
+  vol <- MSGARCH:::f_CondVol(fit$spec,
+                             matrix(fit$par, nrow = 1L, dimnames = list(NULL, names(fit$par))),
+                             data = c(1, 1), do.its = FALSE, nahead = 25L,
+                             ctr = list(nsim = 200L))$vol
+
+  testthat::expect_true(abs(est - mean(vol[21:25])) < 1e-12)
+  testthat::expect_true(abs(est - mean(vol[20:5]))  > 1e-8)
+
+})
+
+testthat::test_that("The mixture Viterbi path is the per-observation MAP decoding", {
+
+  # For a mixture every row of the transition matrix is the same weight vector,
+  # so Viterbi decoding collapses to maximising log(w_k) + log f_k(y_t) at each
+  # observation. (The matrix State() builds used to be filled column-major and
+  # was therefore not row-stochastic; that is fixed, but note it never changed
+  # the decoded path, because the misplaced factor is constant in the index
+  # being maximised over and so cancels.)
+  spec.mix <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH", "sGARCH")),
+                                  switch.spec = list(do.mix = TRUE))
+  par <- spec.mix$par0
+  par["P_1"]      <- 0.8
+  par["alpha0_2"] <- 2
+  y <- SMI[1:300]
+
+  vPath <- as.numeric(State(object = spec.mix, par = par, data = y)$Viterbi)
+
+  mPar <- spec.mix$func$f.do.mix(matrix(par, nrow = 1L))
+  aHt  <- spec.mix$rcpp.func$calc_ht(mPar, y)
+  vW   <- c(par["P_1"], 1 - par["P_1"])
+  mLL  <- sapply(seq_len(spec.mix$K), function(k) {
+    log(vW[k]) + stats::dnorm(y[-1L], 0, sqrt(aHt[2:length(y), 1L, k]), log = TRUE)
+  })
+  vRef <- apply(mLL, 1L, which.max)
+  vRef <- c(vRef[1L], vRef)   # State() copies the first state, see its comment
+
+  testthat::expect_equal(vPath, as.numeric(vRef))
+  testthat::expect_true(all(vPath %in% seq_len(spec.mix$K)))
+
+})
+
+testthat::test_that("CreateSpec validates K and the regime it expands", {
+
+  testthat::expect_error(MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                                             distribution.spec = list(distribution = c("norm", "std")),
+                                             switch.spec = list(K = 3)))
+  testthat::expect_error(MSGARCH::CreateSpec(switch.spec = list(K = 2.5)))
+  testthat::expect_error(MSGARCH::CreateSpec(switch.spec = list(K = 0)))
+
+  # the documented expansion must keep working
+  spec3 <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                               distribution.spec = list(distribution = c("std")),
+                               switch.spec = list(do.mix = FALSE, K = 3))
+  testthat::expect_equal(spec3$name, rep("sGARCH_std", 3L))
+
+})
+
+testthat::test_that("Risk validates its arguments and warns on an inadequate grid", {
+
+  set.seed(1234)
+  fit <- MSGARCH::FitML(spec, data = SMI[1:400], ctr = list(do.se = FALSE))
+  testthat::expect_error(Risk(fit, alpha = 0))
+  testthat::expect_error(Risk(fit, alpha = 1.5))
+  testthat::expect_error(Risk(fit, nahead = 0L))
+  testthat::expect_error(Risk(fit, ctr = list(nmesh = 1L)))
+
+  # the grid spans the observed range, which here holds a sliver of the
+  # predictive mass, so the tail cannot be resolved
+  spec.sr <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                                 distribution.spec = list(distribution = c("norm")),
+                                 switch.spec = list(do.mix = FALSE, K = 1))
+  testthat::expect_warning(Risk(object = spec.sr, par = c(1, 0, 0),
+                                data = c(-0.1, 0.1, -0.05, 0.05), alpha = 0.01))
+
+  # and no warning when the grid does cover the distribution
+  testthat::expect_silent(Risk(fit, alpha = 0.01))
+
+})
+
+testthat::test_that("Prior means and standard deviations are validated", {
+
+  testthat::expect_error(MSGARCH::CreateSpec(prior = list(sd = list(beta_1 = 0))))
+  testthat::expect_error(MSGARCH::CreateSpec(prior = list(sd = list(beta_1 = -1))))
+  testthat::expect_error(MSGARCH::CreateSpec(prior = list(sd = list(wrong_name = 1))))
+  testthat::expect_error(MSGARCH::CreateSpec(prior = list(mean = list(beta_1 = Inf))))
+
+  spec.p <- MSGARCH::CreateSpec(prior = list(mean = list(beta_1 = 0.7),
+                                             sd = list(beta_1 = 0.1)))
+  testthat::expect_equal(unname(spec.p$rcpp.func$get_sd()[3L]), 0.1)
+
+})
+
+testthat::test_that("Fixed parameters are honoured and transition probabilities refused", {
+
+  testthat::expect_error(MSGARCH::CreateSpec(constraint.spec = list(fixed = list(P_1_1 = 0.99))))
+  testthat::expect_error(MSGARCH::CreateSpec(constraint.spec = list(fixed = list(beta_1 = NA))))
+
+  spec.f <- MSGARCH::CreateSpec(constraint.spec = list(fixed = list(beta_1 = 0.8)))
+
+  set.seed(1234)
+  fit <- MSGARCH::FitML(spec.f, data = SMI[1:400], ctr = list(do.se = FALSE))
+  testthat::expect_true(abs(fit$par["beta_1"] - 0.8) < 1e-12)
+
+  # a user-supplied par0 must not shift the sampler's parameter vector, and the
+  # identification sort must not relabel the regime the constraint refers to
+  set.seed(1234)
+  mcmc <- suppressMessages(MSGARCH::FitMCMC(spec.f, data = SMI[1:400],
+                                            ctr = list(par0 = spec.f$par0, nburn = 50L,
+                                                       nmcmc = 100L, nthin = 1L)))
+  mPar <- as.matrix(mcmc$par)
+  testthat::expect_equal(ncol(mPar), length(spec.f$label))
+  testthat::expect_true(max(abs(mPar[, "beta_1"] - 0.8)) < 1e-12)
+
+})
diff --git a/Package/tests/testthat/test_TimeSeriesIndex.R b/Package/tests/testthat/test_TimeSeriesIndex.R
new file mode 100644
index 0000000..e3d3d15
--- /dev/null
+++ b/Package/tests/testthat/test_TimeSeriesIndex.R
@@ -0,0 +1,76 @@
+testthat::context("Test ts and zoo handling of fitted-object methods")
+
+data("SMI", package = "MSGARCH")
+spec <- MSGARCH::CreateSpec(variance.spec = list(model = c("sGARCH")),
+                            distribution.spec = list(distribution = c("norm")),
+                            switch.spec = list(do.mix = FALSE, K = 2))
+
+y.num <- as.numeric(SMI[1:200])
+y.ts  <- stats::ts(y.num, start = c(2000, 1), frequency = 12)
+newdata <- c(0.1, 0.2)
+
+set.seed(1234)
+fit.ts  <- MSGARCH::FitML(spec, data = y.ts,  ctr = list(do.se = FALSE))
+set.seed(1234)
+fit.num <- MSGARCH::FitML(spec, data = y.num, ctr = list(do.se = FALSE))
+
+testthat::test_that("Combining a fit's data with newdata keeps index and values aligned", {
+
+  # The index used to be derived from the already-concatenated series and then
+  # extended by another length(newdata) points, so zooreg() recycled
+  # observations and the model conditioned on values that are not in the sample.
+  combined <- MSGARCH:::f_combine_data(y.ts, newdata)
+  testthat::expect_equal(length(combined), length(y.num) + length(newdata))
+  testthat::expect_equal(as.numeric(combined), c(y.num, newdata))
+  testthat::expect_equal(stats::frequency(combined), stats::frequency(y.ts))
+  testthat::expect_equal(stats::start(combined), stats::start(y.ts))
+
+  # without newdata the series must come back untouched, index included
+  testthat::expect_equal(stats::tsp(MSGARCH:::f_combine_data(y.ts, NULL)), stats::tsp(y.ts))
+
+  # and a plain numeric series is just concatenated
+  testthat::expect_equal(MSGARCH:::f_combine_data(y.num, newdata), c(y.num, newdata))
+
+})
+
+testthat::test_that("ts input gives the same numbers as the equivalent numeric input", {
+
+  testthat::expect_equal(as.numeric(Volatility(fit.ts, newdata = newdata)),
+                         as.numeric(Volatility(fit.num, newdata = newdata)))
+  testthat::expect_equal(length(Volatility(fit.ts, newdata = newdata)),
+                         length(y.num) + length(newdata))
+
+  testthat::expect_equal(as.numeric(predict(fit.ts,  nahead = 1L)$vol),
+                         as.numeric(predict(fit.num, nahead = 1L)$vol))
+  testthat::expect_equal(as.numeric(PIT(fit.ts,  do.its = TRUE)),
+                         as.numeric(PIT(fit.num, do.its = TRUE)))
+  testthat::expect_equal(as.numeric(Risk(fit.ts,  alpha = 0.05)$VaR),
+                         as.numeric(Risk(fit.num, alpha = 0.05)$VaR))
+
+})
+
+testthat::test_that("Forecast indexes advance by the series' own time step", {
+
+  pred <- predict(fit.ts, nahead = 3L)$vol
+  testthat::expect_equal(stats::frequency(pred), stats::frequency(y.ts))
+  # three monthly steps beyond the end of the sample, not three whole years
+  testthat::expect_equal(as.numeric(stats::time(pred)),
+                         stats::tsp(y.ts)[2L] + (1:3) / stats::frequency(y.ts))
+
+  vol <- Volatility(fit.ts)
+  testthat::expect_equal(stats::tsp(vol)[1:2], stats::tsp(y.ts)[1:2])
+
+})
+
+testthat::test_that("zoo input is handled the same way", {
+
+  y.zoo <- zoo::zoo(y.num, order.by = seq_len(length(y.num)))
+  set.seed(1234)
+  fit.zoo <- MSGARCH::FitML(spec, data = y.zoo, ctr = list(do.se = FALSE))
+
+  testthat::expect_equal(as.numeric(Volatility(fit.zoo, newdata = newdata)),
+                         as.numeric(Volatility(fit.num, newdata = newdata)))
+  testthat::expect_equal(as.numeric(zoo::index(predict(fit.zoo, nahead = 2L)$vol)),
+                         length(y.num) + 1:2)
+
+})

From ea0847dffd15299dd984e8dec69786deab338bae Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 20:06:15 -0400
Subject: [PATCH 19/20] Record the robustness and validation fixes in NEWS

---
 Package/NEWS | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/Package/NEWS b/Package/NEWS
index 3e6a0f4..8c5ad4d 100644
--- a/Package/NEWS
+++ b/Package/NEWS
@@ -1,4 +1,18 @@
 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

From 0d4f7dc0509131eae0a46f9902ab9fe17047d41c Mon Sep 17 00:00:00 2001
From: David Ardia <ArdiaD@users.noreply.github.com>
Date: Wed, 19 Aug 2026 20:23:02 -0400
Subject: [PATCH 20/20] Update the pull-request description for the full change
 set

Covers all three rounds of fixes rather than only the first five, groups them
by the kind of defect, and records what changes for existing users.
---
 PR_BODY.md | 339 ++++++++++++++++++++++++++++++++++++++---------------
 1 file changed, 244 insertions(+), 95 deletions(-)

diff --git a/PR_BODY.md b/PR_BODY.md
index d36a2e5..2614fdf 100644
--- a/PR_BODY.md
+++ b/PR_BODY.md
@@ -1,124 +1,273 @@
-Five bugs in the reporting and persistence layer around the likelihood, plus regression
-tests for each. The likelihood itself is not touched: no estimate, log-likelihood,
-conditional variance or state probability changes anywhere in this branch. What changes is
-what `summary()` prints, what `AIC`/`BIC` count, what `predict()` returns for a Bayesian
-fit, and whether a saved fit can be reloaded at all.
+Twenty-five fixes to the layer around the likelihood — inference, persistence, forecasting,
+time-index handling, input validation and numerical robustness — with regression tests for
+each. 18 commits, one per defect or coherent group, so any of them can be reviewed or dropped
+on its own.
 
-Everything below is reproduced on `data("SMI")` with the shipped code. `R CMD check
---as-cran` is unchanged by the branch (2 WARNINGs, 1 NOTE, `testthat` OK — and the version
-bump in the last commit clears one of the two WARNINGs).
+**The likelihood itself is untouched.** On the SMI example the log-likelihood at a fixed
+parameter vector is bit-identical to `master` (`-3389.2962430913985`), the estimates agree to
+eight significant figures, and AIC/BIC for an unconstrained model are unchanged. What changes
+is what `summary()` prints, what `predict()` returns for a Bayesian fit, whether a saved fit
+can be reloaded, and what happens at the edges.
 
-## 1. Standard errors used the delta method transposed — `R/Inference.R:38`
+`R CMD check --as-cran` on the branch: **1 WARNING, 1 NOTE**, both environmental (an Apple
+clang warning raised inside R's own `R_ext/Boolean.h`, and "unable to verify current time"
+offline). `master` reports 2 WARNINGs and 1 NOTE; the version bump clears one WARNING and
+removing dead code clears the S3 NOTE.
+
+---
+
+## 1. Wrong numbers reported to the user
+
+### Standard errors used the delta method transposed — `R/Inference.R:38`
 
 ```r
 mSandwitch <- t(mJacob) %*% mInvHessian %*% mJacob   # -> mJacob %*% mInvHessian %*% t(mJacob)
 ```
 
-`numDeriv::jacobian` returns `∂f_i/∂x_j`, so `Var(g(θ̂)) = J V J'`. The two orientations
-agree only if `J` is symmetric, and it is not: the working→natural map is triangular inside
-each regime (the sGARCH bound on `beta` is `0.9999 − alpha1`; the gjrGARCH and tGARCH bounds
-on `beta` also involve `alpha2` and the shape/skew parameters), and the transition-probability
-block is *anti*-diagonal.
+`numDeriv::jacobian` returns `∂f_i/∂x_j`, so `Var(g(θ̂)) = J V J'`. The two orientations agree
+only if `J` is symmetric, and it is not: the working→natural map is triangular inside each
+regime (the sGARCH bound on `beta` is `0.9999 − alpha1`; the gjr/tGARCH bounds also involve
+`alpha2` and the shape/skew parameters), and the transition-probability block is *anti*-diagonal.
+Checked against the observed information computed directly in the natural parameterisation:
 
-On the default MS(2)-GARCH(1,1)-Normal fit to `SMI`, checked against the observed information
-computed directly in the natural parameterisation:
-
-| | natural-scale `H` | fixed (`J V J'`) | before (`J' V J`) |
+| | natural-scale `H` | this PR | `master` |
 |---|---|---|---|
-| `alpha1_1` | 0.01606 | 0.01511 | **0.03413** |
-| `beta_1`   | 0.02197 | 0.02091 | **0.00958** |
-| `alpha1_2` | 0.00431 | 0.00437 | **0.00610** |
-| `beta_2`   | 0.00416 | 0.00426 | **0.00049** |
-| `P_1_1`    | 0.00982 | 0.00973 | **0.00059** |
-| `P_2_1`    | (at bound) | 0.03072 | **0.50250** |
+| `alpha1_1` | 0.01606 | 0.015111 | **0.034126** |
+| `beta_1`   | 0.02197 | 0.020912 | **0.009584** |
+| `alpha1_2` | 0.00431 | 0.004373 | **0.006103** |
+| `beta_2`   | 0.00416 | 0.004256 | **0.000492** |
+| `P_1_1`    | 0.00982 | 0.009727 | **0.000595** |
+| `P_2_1`    | (at bound) | 0.030683 | **0.502501** |
+
+Six of eight were wrong, by factors from 0.06× to 16×, and because the P block is anti-diagonal
+the two transition probabilities had their standard errors **exchanged**. Only `alpha0_k`, whose
+map is a plain `exp`, was unaffected.
+
+### `Pr(>|t|)` was one-sided, then cancelled in the tail — `R/Inference.R:43`
 
-Six of eight were wrong, by factors from 0.06× to 16×. Because the P block is anti-diagonal,
-transposing **exchanged the two transition probabilities' standard errors**. Only `alpha0_k`,
-whose map is a plain `exp`, was unaffected.
+`1 - pnorm(abs(t))` under a two-sided label. Now `2 * pnorm(-abs(t))`, evaluated on the lower
+tail: `2 * (1 - pnorm(abs(t)))` cancels to exactly zero once `|t|` exceeds about 8.3, while the
+lower tail stays representable to about 38. The single-regime GARCH fit used in the tests has a
+`t` of 26.9, whose p-value was reported as 0 and is now 1.6e-159.
 
-## 2. `Pr(>|t|)` was one-sided — `R/Inference.R:43`
+### `predict()` on an MCMC fit returned draw #1 — `R/CondVol.R:32`
 
-`1 - pnorm(abs(t))` under a two-sided label; now `2 * (1 - pnorm(abs(t)))`.
+`vol` is `(T+1) × ndraw` and the one-step value was `vol[dim(PredProb)[1]]` — a single index
+into a matrix is linear indexing, i.e. the last row of the *first column*. On a 100-draw chain
+fitted to SMI the reported value was 1.019753 (the first draw) against a posterior mean of
+1.039776, with a spread of 0.999837–1.097337 across draws. `Volatility()` already averaged
+correctly, so the two methods disagreed on the same fit.
 
-## 3. `AIC`/`BIC` mis-counted `regime.const.pars` — `R/Utils.R:420`
+### `AIC`/`BIC` mis-counted `regime.const.pars` — `R/Utils.R:420`
 
-A regime-constant parameter leaves one free value where there were `K`, so it removes `K − 1`
-degrees of freedom (cf. `f_rename_par`, which strips `name_2 … name_K`). `dofMSGARCH`
-subtracted one, which is right only at `K = 2`: with a regime-constant shape parameter the df
-was 17 instead of 16 at `K = 3` and 27 instead of 25 at `K = 4`, always over-penalising the
-constrained model.
+A regime-constant parameter removes `K − 1` degrees of freedom, not one (cf. `f_rename_par`,
+which strips `name_2 … name_K`). Right only at `K = 2`: with a regime-constant shape parameter
+the df was 17 instead of 16 at `K = 3` and 27 instead of 25 at `K = 4`, always over-penalising
+the constrained model.
 
-## 4. Saved specs and fits could not be reloaded — `R/Utils.R:349-353`
+---
+
+## 2. Failures and silent data corruption
+
+### Saved specs and fits could not be reloaded — `R/Utils.R:349-353`
 
 R serializes external pointers as `NULL`, so a `saveRDS`-ed spec or fit comes back with dead
 pointers. `f_check_spec` exists to rebuild them, but the rebuild branch called
-`spec$rcpp.func$get_mean()` / `get_sd()` — the very pointer whose failure had just triggered
-the branch:
+`spec$rcpp.func$get_mean()` / `get_sd()` — the very pointer whose failure had just triggered it:
 
 ```
 Volatility(fit) : Error in .External(...): NULL value passed as symbol address
-State(fit)      : Error in .External(...): NULL value passed as symbol address
-predict(fit)    : Error in .External(...): NULL value passed as symbol address
 ```
 
-The two values read there are already held on the R side in `spec$prior.mean` /
-`spec$prior.sd`, which the next two lines were using anyway, so the C++ round trip was dead
-code. Removing it makes the rebuild work; a reloaded spec, ML fit or MCMC fit now returns
-values identical to before it was saved, with user priors preserved.
+Those two values are already held on the R side in `spec$prior.mean` / `spec$prior.sd`, which
+the next two lines were using anyway, so the C++ round trip was dead code. This is the ordinary
+workflow of fitting a model, saving it, and analysing it later — or sending a spec to a
+`parLapply` worker.
+
+### `ts`/`zoo` methods recycled observations when `newdata` was supplied
+
+`Volatility`, `predict`, `PIT`, `PredPdf` and `Risk` concatenated `object$data` with `newdata`,
+then derived the index from the *already-concatenated* series and appended another
+`length(newdata)` points. The index was longer than the values, so `zooreg()` recycled: a
+200-point monthly series plus two new returns produced **204** observations ending in two values
+copied from the start of the sample, and the model conditioned on them. The same code discarded
+the original `start` and `frequency` even when `newdata` was `NULL` (a monthly series became
+annual), and forecast indexes advanced by one index unit rather than by the series' own step.
+Three helpers in `Utils.R` — `f_combine_data`, `f_future_index`, `f_index_result` — now handle
+this in one place, replacing ten duplicated blocks. `ts` input now gives numerically identical
+results to the equivalent numeric input.
+
+### A singular stationary solve destroyed the whole fit — `src/MSgarch.h:283`
+
+`loadparam` obtained the stationary distribution with a plain matrix inverse, on every
+likelihood evaluation. The plain parameter mapping bounds each free transition probability
+separately, so with `K ≥ 3` a row can leave the simplex and the matrix is singular; the uncaught
+Armadillo exception then killed the run. `FitMCMC` always uses that mapping, so every chain with
+three or more regimes was exposed, as was any ML fit using `fixed` or `regime.const` parameters
+— `FitML` on a `K = 3` spec with `regime.const = "nu"` died deterministically at
+`P = [[1, 0.999973, −0.999973], [0,1,0], [0,0,1]]`. The solve is now guarded and falls back to
+the uniform distribution, which `calc_prior` rejects a moment later anyway.
+
+### `do.sort` silently broke `constraint.spec$fixed`
+
+The identification sort relabels regimes by unconditional variance, which moves a parameter
+fixed in one regime into another. With the default `do.sort = TRUE`, a parameter fixed at 0.8
+came back taking **sixteen different values** across a 100-draw chain. The sort is now skipped
+when parameters are fixed, with a message.
+
+### `FitMCMC` mis-mapped a user `ctr$par0` under `fixed.pars`
+
+Unlike `FitML`, `FitMCMC` never dropped the fixed entries, so the vector handed to the sampler
+was too long, `f_rename_par` left a trailing `NA` name, and `f_mapPar` looked bounds up by name
+— mapping `alpha0_2` to 80.008 and the last parameter to `NA`.
+
+### `constraint.spec$fixed` accepted transition probabilities that never worked
+
+`CreateSpec` validates against `out$label`, which includes `P_1_1`, but the starting-value
+routine hands `"P_1_1"` to a single-regime `CreateSpec` that has no such parameter, and the
+prior correction in `Kernel()` indexes `prior.mean`, which covers only the within-regime
+coefficients — so the log-posterior became `NA` and was floored to `-1e10`. Now refused with an
+explanation.
+
+---
+
+## 3. Numerical robustness
+
+### The Hamilton filter overflowed to `NaN` — `src/MSgarch.h`
 
-This is the ordinary workflow of fitting a model, saving it, and analysing it in a later
-session — or shipping a spec to a `parLapply` worker.
+Both filters shifted each column of regime log densities by its *smallest* entry, and only when
+that fell below `log(DBL_MIN)`. That guards underflow but creates an overflow: the largest
+exponent becomes `max − min − 707`, so once the regimes differ by more than about 1400 in log
+density the exponential returns `Inf` and normalising gives `NaN`. The gap is reachable with
+parameters the package itself accepts — `alpha0 = 3.3e-4` against `alpha0 = 1` is enough — after
+which the likelihood floors to `-1e10` and `State()`, `PredPdf()` and `Risk()` return `NaN` or
+fail. Both now use the standard log-sum-exp anchor, the largest entry.
 
-## 5. `predict()` on an MCMC fit returned draw #1 — `R/CondVol.R:32`
+Verified against the analytic limit rather than merely "no longer `NaN`": as one regime's density
+vanishes the likelihood converges to `log P(other regime) + log f(y)`.
 
-`vol` is `(T+1) × ndraw` and the one-step-ahead value was `vol[dim(PredProb)[1]]` — a single
-index into a matrix is linear indexing, i.e. the last row of the *first column*. On a
-100-draw chain fitted to `SMI` the reported value was 1.019753 (the first draw) against a
-posterior mean of 1.039776, with a spread of 0.999837–1.097337 across draws. `Volatility()`
-already averaged correctly, so the two methods disagreed on the same fit. The
-single-parameter (ML) path is unchanged.
+### The multi-regime in-sample CDF wrote the first slice transposed — `src/MSgarch.h:497`
+
+`f_cdf_its` filled the `t = 0` slice of `arma::cube tmp(ny, nx, K)` with `tmp(ix, 0, s)` where
+everything else in the file uses `(0, ix, s)`. With a grid shorter than the sample the first
+observation's CDF was zero for every point but the first; with a longer grid the write ran off
+the cube. The new test checks the first row against a closed form — at `t = 1` each regime's
+conditional variance is its unconditional variance and the predictive state distribution is the
+ergodic one — and it now matches exactly.
+
+### The native log branches returned the last regime — `src/MSgarch.h:399,474`
+
+`f_pdf` and `f_cdf` accumulate the state-weighted mixture in `out` but under `is_log` overwrote
+it with `log(tmp[i])`, where `tmp` holds only the regime evaluated last: `exp()` of the returned
+value differed from the mixture by factors of 12 and 26 in the tails. Not reachable from
+`PredPdf`/`PIT`, which pass `is_log = FALSE`, but both are exposed as module methods on
+`spec$rcpp.func`.
+
+---
+
+## 4. Validation and diagnostics
+
+- **`FitML` reported every failure identically.** The guard tested `llk == 1e+10`, but `f_nll`
+  returns `+1e10` so a failure arrives as `-1e10`; `f_OptimFUNDefault` also wraps `optim` in
+  `try()`, so `optimizer$value` errored first. Everything — a singular transition matrix, a
+  malformed spec, one `NA` in the data — surfaced as
+  `$ operator is invalid for atomic vectors`.
+- **Data and parameters.** `f_check_y`/`f_check_par` rejected input only when it was *entirely*
+  `NaN`, so a single `NA`, `NaN` or `Inf` reached the compiled code and `Volatility()` returned
+  a complete, plausible-looking series computed from corrupt data. Both now require finite
+  values, and the data must hold at least two observations. `f_nll`/`f_posterior` treat a
+  non-finite mapped parameter vector as an infeasible point so strictness cannot abort a fit.
+- **`Risk()`.** `alpha`, `nahead` and `ctr$nmesh` are validated (`alpha = 1.5` used to return
+  6.6745 silently; `alpha = 0` gave `ES = -Inf`). The evaluation grid spans the observed data
+  range, so it can miss part of the predictive distribution; when the omitted mass exceeds the
+  requested tail probability, the "quantile" is the grid boundary, and that now warns. On a
+  normal fit to SMI the grid's left endpoint has CDF 7e-8 and nothing changes; on a four-point
+  sample it covers 15% of the distribution and warns.
+- **`prior$sd`** was validated by the *mean* checker; `sd = 0` and `sd = -1` were accepted and
+  reached C++. Now checked for finiteness and positivity, and both prior validators produce
+  usable messages instead of `stop(cat(...))`.
+- **`CreateSpec`** validates `switch.spec$K`, and when expanding one regime through `K` rejects
+  an explicitly heterogeneous distribution vector — the guard tested
+  `distribution.spec$model`, which does not exist, so `model = "sGARCH"`,
+  `distribution = c("norm","std")`, `K = 3` silently produced `norm, std, norm`.
+- **`simulate()`** accepts `nburn = 0`, which used to drop the first draw (`1:0` is `c(1, 0)`)
+  and then fail on the dimnames.
+- **`UncVol()`** averages the horizons *after* the burn-in; `nburn:nahead` is a descending range
+  under the shipped defaults and in general keeps part of the transient.
+
+---
+
+## 5. Housekeeping
+
+- The `Sim.MSGARCH_ML_FIT` / `Sim.MSGARCH_MCMC_FIT` methods behind the standing `R CMD check`
+  S3 NOTE are unreachable — `Sim` is not exported, every internal caller passes a spec, and
+  `simulate.*_FIT` pass `object$spec` explicitly. Removed; the NOTE is gone.
+- `Rcpp:::LdFlags()` dropped from both `Makevars`; Rcpp has not needed it since 2013.
+- The shipped BIC test asserted `abs(exp.BIC - exp.BIC) < tol`, which is zero by construction.
+  (The value it meant to check was correct, so this is a test fix, not a bug fix.)
+- The unused log-likelihood accumulator in `f_get_Pstate` is removed.
+- The mixture transition matrix built for Viterbi decoding is now row-stochastic. **The decoded
+  path is unchanged** — the misplaced factor is constant in the index being maximised over, so
+  it cancels — but the matrix should not have rows summing to 1.6 and 0.4.
+
+---
 
 ## Tests
 
-`test_Inference.R` and `test_Serialization.R` are new; `test_Volatility.R` gains one block.
-11 blocks, 34 assertions, +13s of check time. All 11 fail on the current code.
-
-Two things worth pointing out, since they are what makes the tests worth having:
-
-- The standard-error test is anchored on a **single-regime** GARCH(1,1)-Normal, not on the
-  MS(2) default. Every parameter there is interior, so a central-difference Hessian of the
-  natural-scale negative log-likelihood is well conditioned: it agrees with `J V J'` to 1e-4
-  in relative terms while the transposed sandwich is off by 44% and 85%. That is an
-  independent check on the *value*, not a restatement of the formula.
-- Every block opens with a guard asserting its own precondition — the two sandwich
-  orientations really differ for this model; the round trip really did invalidate the
-  pointers; the posterior mean really differs from the first draw — so none of them can pass
-  vacuously if the surrounding code changes.
-
-The `K = 3` and `K = 4` legs of the degrees-of-freedom test run against a fit-shaped list
-rather than a real fit, because a constrained `K ≥ 3` model cannot currently be fitted at
-all (see below).
-
-## Not in this branch
-
-The last commit (version bump + `NEWS`) is separable — drop it if the release number should
-be decided elsewhere; nothing depends on it.
-
-A review of the package turned up nine further issues that are **not** addressed here,
-several of them more serious than some of the above. The two worth flagging now:
-
-- **`src/MSgarch.h:283`** computes the ergodic distribution with a raw Armadillo `.i()` on
-  `I − P + U`, on every likelihood evaluation. When `do.plm = TRUE` — forced by `fixed.pars`
-  and `regime.const.pars`, and hard-coded in `FitMCMC` — the free transition entries are
-  mapped into `(0,1)` independently, so for `K ≥ 3` a row can leave the simplex, the matrix
-  can be exactly singular, and the uncaught exception aborts the entire run. Reproducible:
-  `FitML` on a `K = 3` spec with `regime.const = "nu"` dies at
-  `P = [[1, 0.999973, −0.999973], [0,1,0], [0,0,1]]`. This is why the tests above cannot fit
-  a constrained `K ≥ 3` model.
-- **`R/FitML.R:133`** guards optimisation failure with `if (llk == 1e+10)`, but `f_nll`
-  returns `+1e10` so `llk` is `−1e10`; and `f_OptimFUNDefault` wraps `optim` in `try()`, so
-  `optimizer$value` errors first. Every failure mode — including a single `NA` in the data,
-  which `f_check_y` lets through — therefore surfaces as
-  `Error in optimizer$value : $ operator is invalid for atomic vectors`.
-
-Happy to open these as separate issues or as a follow-up PR, whichever you prefer.
+Five new files and additions to two existing ones: 28 blocks and 103 assertions. The suite
+goes from 23 blocks / 23 assertions to 52 blocks / 128 assertions, adding about 25s of
+check time.
+
+- `test_Inference.R` — pins the delta method against an independently computed observed
+  information. Anchored on a single-regime GARCH(1,1)-Normal because every parameter there is
+  interior, so a central-difference Hessian of the natural-scale negative log-likelihood is well
+  conditioned: it agrees with `J V J'` to 1e-4 relative while the transposed sandwich is off by
+  44% and 85%. Also the two-sided and tail-precision checks, and the degrees-of-freedom
+  arithmetic for `K = 2, 3, 4`.
+- `test_Serialization.R` — round-trips a spec, an ML fit and an MCMC fit through
+  `saveRDS`/`readRDS` and requires identical results with priors preserved.
+- `test_NativeDensity.R` — the in-sample CDF against a closed form, a grid longer than the
+  sample, and the log branches against `log()` of the linear ones.
+- `test_Robustness.R` — the filter under a log-density gap over 1000 including the analytic
+  limit; a singular transition matrix; non-finite and too-short input; `nburn = 0`; the `UncVol`
+  window; the mixture Viterbi path against per-observation MAP decoding; `Risk`'s guards and
+  warning; prior validation; and that a fixed parameter survives `FitML`, a user `par0` and the
+  identification sort.
+- `test_TimeSeriesIndex.R` — index/value alignment, `ts` and `zoo` against numeric, and forecast
+  index spacing.
+
+Every block opens with a guard asserting its own precondition — the two sandwich orientations
+really differ for this model, the round trip really did invalidate the pointers, the posterior
+mean really differs from the first draw — so none of them can pass vacuously if the surrounding
+code changes. All fail on the commit preceding their fix.
+
+---
+
+## Compatibility
+
+Deliberate changes in output, all of them corrections:
+
+- Standard errors and p-values from `summary()` change for every model with a `beta` or a
+  transition probability.
+- `predict()` on an MCMC fit returns the posterior mean rather than the first draw.
+- `AIC`/`BIC` change for `K ≥ 3` models using `regime.const`.
+- `ts` and `zoo` users get results of the correct length, with the correct index, from
+  `Volatility`, `predict`, `PIT`, `PredPdf` and `Risk`.
+- Models using `constraint.spec$fixed` are no longer regime-sorted, so regime labelling may
+  differ from previous fits.
+
+Newly rejected input that used to be accepted: data containing `NA`, `NaN` or `Inf`; data
+shorter than two observations; non-finite parameters; `alpha` outside `(0,1)`; `ctr$nmesh < 2`;
+non-positive `prior$sd`; `constraint.spec$fixed` on a transition probability; and a
+heterogeneous distribution vector combined with `switch.spec$K`.
+
+Because inference and prediction results change, the reverse imports (`MSGARCHelm`, `SBAGM`)
+are worth a check before release.
+
+The version bump to 2.52 and the `Date` refresh are in their own commit (`c93e218`) — drop it if
+the release number should be set separately; nothing else depends on it.
+
+`REVIEW.md` and `REVIEW_codex.md` in the diff are the working notes behind these fixes: a
+defect-by-defect write-up with reproductions, and the verbatim output of an independent
+read-only audit used to cross-check it. Happy to drop them from the branch if you would rather
+they not live in the repository.