diff --git a/PR_BODY.md b/PR_BODY.md
new file mode 100644
index 0000000..2614fdf
--- /dev/null
+++ b/PR_BODY.md
@@ -0,0 +1,273 @@
+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.
+
+**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.
+
+`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 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:
+
+| | natural-scale `H` | this PR | `master` |
+|---|---|---|---|
+| `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`
+
+`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.
+
+### `predict()` on an MCMC fit returned draw #1 — `R/CondVol.R:32`
+
+`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.
+
+### `AIC`/`BIC` mis-counted `regime.const.pars` — `R/Utils.R:420`
+
+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.
+
+---
+
+## 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 it:
+
+```
+Volatility(fit) : Error in .External(...): NULL value passed as symbol address
+```
+
+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`
+
+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.
+
+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)`.
+
+### 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
+
+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.
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..8c5ad4d 100644
--- a/Package/NEWS
+++ b/Package/NEWS
@@ -1,3 +1,29 @@
+Changes in Version 2.52
+  o The Hamilton filter now stabilises the regime log densities on their largest entry; anchoring on the smallest one overflowed to NaN when the regimes differed by more than about 1400 in log density
+  o The stationary distribution is obtained with a guarded solve: a transition matrix outside the simplex, reachable whenever fixed or regime-constant parameters are used with three or more regimes, no longer aborts the fit with an Armadillo exception
+  o Fitted-object methods no longer recycle observations when newdata is supplied with ts or zoo data, and forecast indexes advance by the series' own time step
+  o simulate() accepts nburn = 0
+  o FitML reports optimizer failures instead of failing on '$ operator is invalid for atomic vectors'
+  o UncVol averages the horizons after the burn-in rather than a descending range inside it
+  o The data must now contain at least two observations and no NA, NaN or infinite values; parameters are checked likewise
+  o Risk() validates alpha, nahead and ctr$nmesh, and warns when the evaluation grid does not cover enough of the predictive distribution to resolve the requested tail
+  o prior$sd is validated with the standard-deviation checker and must be finite and strictly positive
+  o Transition probabilities can no longer be passed to constraint.spec$fixed, which never worked
+  o FitMCMC drops fixed parameters from a user-supplied ctr$par0, as FitML already did
+  o The identification sort is skipped when constraint.spec$fixed is used, since relabelling regimes would move the fixed value into another regime
+  o CreateSpec validates switch.spec$K and rejects an explicitly heterogeneous distribution vector when expanding a single regime
+  o The mixture transition matrix built for Viterbi decoding is now row-stochastic (the decoded path is unchanged)
+  o Fixed the standard errors reported by summary(): the delta-method sandwich was transposed
+  o Pr(>|t|) in summary() is now the two-sided p-value, as its label states
+  o AIC/BIC now drop K-1 degrees of freedom per regime-constant parameter, not one
+  o Saved specifications and fits (saveRDS/readRDS) are usable again: the Rcpp modules are rebuilt
+  o predict() on an MCMC fit now averages over the posterior draws instead of returning the first draw
+  o Fixed a transposed cube write in the multi-regime in-sample CDF, which corrupted the first row of PIT(do.its = TRUE) and could overrun the buffer
+  o The is_log branches of the native mixture density and CDF now return the mixture instead of the last regime
+  o Two-sided p-values are evaluated on the lower tail, so they no longer cancel to zero for large statistics
+  o Removed the unreachable Sim methods for fit objects, which R CMD check flagged as unregistered
+  o Dropped the obsolete Rcpp:::LdFlags() call from src/Makevars
+  o Added regression tests for all of the above
 Changes in Version 2.51
   o Fix warning: use of bitwise '|' with boolean operands
 Changes in Version 2.5
diff --git a/Package/R/CondVol.R b/Package/R/CondVol.R
index fa0318d..ac0ac82 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) {
@@ -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/Inference.R b/Package/R/Inference.R
index 1196a0c..6971bad 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 * pnorm(-abs(vTest))
   
   out[, "Estimate"]   <- vPn
   out[, "Std. Error"] <- vSE
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 8ffd8d8..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()
@@ -345,15 +413,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)
 }
@@ -403,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)
 }
@@ -411,12 +486,19 @@ 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)
 }
 
 
 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"]]))
 }
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 88cb305..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
@@ -175,19 +172,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)
-}
diff --git a/Package/src/MSgarch.h b/Package/src/MSgarch.h
index 22df1cb..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);
   }
@@ -396,7 +407,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 +482,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]);
     }
   }
   
@@ -494,7 +505,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++;
   }
@@ -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);
   
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)
diff --git a/Package/tests/testthat/test_Inference.R b/Package/tests/testthat/test_Inference.R
new file mode 100644
index 0000000..d7e636f
--- /dev/null
+++ b/Package/tests/testthat/test_Inference.R
@@ -0,0 +1,173 @@
+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 * 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", {
+
+  # 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_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)
   
 })
 
diff --git a/Package/tests/testthat/test_NativeDensity.R b/Package/tests/testthat/test_NativeDensity.R
new file mode 100644
index 0000000..61060c8
--- /dev/null
+++ b/Package/tests/testthat/test_NativeDensity.R
@@ -0,0 +1,59 @@
+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)
+  )
+
+})
+
+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)
+
+})
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_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_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)
+
+})
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
diff --git a/REVIEW.md b/REVIEW.md
new file mode 100644
index 0000000..5dfacf4
--- /dev/null
+++ b/REVIEW.md
@@ -0,0 +1,604 @@
+# 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.
+>
+> **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
+> 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` — **FIXED**
+
+```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` — **FIXED**
+
+```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.