fix(engine): correct BoxCox lambda selection (diffuse Jacobian + vari…

CLAUDE.md

@@ -129,6 +129,45 @@ Header-only, all included from `src/musecpp2R.cpp`:
 - **Regressor matrix orientation.** User supplies (n × k); internally and in C++ it is (k × n).  The transpose happens once in `.pts_parse_data()` for training data and once in `.pts_forecast_inputs()` for `newdata`.
 - **Harmonic periods.** `lagsAll = lags / (1 : floor(lags/2))`.  The C++ engine may select a subset, but `lagsAll` always stores the full candidate set that was passed in.
 
+### Diffuse concentrated likelihood + the ARMA sign convention (the AR-on-irregular saga)
+
+Two coupled things, both settled — read before touching `llik()` (`src/SSpace.h`)
+or `polyStationary`/`invPolyStationary` (`src/ARMAmodel.h`).
+
+**1. The diffuse term is the exact-diffuse concentrated BCnorm likelihood.**
+`llik()`'s diffuse block carries `-0.5*nDiff*(log(2pi) + log(s2)) - 0.5*diffuseTerm`
+AND the Box-Cox Jacobian over all `n` observations.  The `-0.5*nDiff*log(s2)`
+piece is REQUIRED: the KF runs in concentrated (ratio) units, so each diffuse
+obs' determinant is `log|Finf,t| = log(s2) + log(Finf*,t)`; `s2 = SSR/n` is the
+MLE for this (all-`n`) form.  Dropping `log(s2)` (or the diffuse Jacobian) leaves
+a dangling monotone-in-lambda term that drove the joint Box-Cox lambda MLE to a
+bound (iid N(1000,50) -> lambda=2 not 1; Nile -> lambda=2 + spurious trend).
+As `s2 -> 0` on a degenerate short-series fit this term -> +Inf ("+986 positive
+logLik"); a weakly-informative inverse-gamma ridge on `s2`
+(`s2 = (SSR + nu0*s0sq)/(n+nu0)`, `nu0=1`, `s0sq = 1e-4*Var(g(y))`) bounds it.
+
+**2. The AR sign flip was an INIT/CONVENTION bug, not the likelihood.**
+On stationary AR(1)=+phi data, `pts(y, model="1NN", orders=list(ar=1))` used to
+return `AR(1) ~ -phi` with a non-collapsed `$B["Level"]`.  Cause: the ARMA seed
+(sample PACF, standard `(1 - phi B)`) was fed to `invPolyStationary`, which then
+used the internal `(1 + phi B)` form, so a correct `+phi` seed landed in the
+`-phi` basin.  For pure-structural models the optimiser escaped it; under the
+(correct) diffuse `log(s2)` term that basin is a local min it can't leave, so the
+RW level stays active and the AR flips (the RW-level-vs-AR identification
+manifold: over-detrending a positively-autocorrelated series induces negative
+lag-1 autocorrelation -- `arima(diff(y),c(1,0,0))` confirms it).
+
+**The fix (do NOT revert):** `polyStationary`/`invPolyStationary` now work in the
+STANDARD conventions -- AR `(1 - phi B)`, MA `(1 + theta B)` -- so the internal
+parameter equals the reported coefficient.  Consequences kept in sync:
+`armaMatrices` builds `T = +col` (AR) and `R = +col` (MA); `musecore.h` reporting
+and `parameterValues` emit the coefs with no sign flip; `setEstimatedParams` does
+no AR flip.  The MA seed clamps its `arToPacf` input to the invertible region
+(`|PACF| < 1`) exactly like the AR seed -- without it a 2-block SARMA MA seed
+crashes in `invPolyStationary` (division by `1 - phi^2` at the unit boundary).
+Verified: AR(1)=0.7 -> +0.67, MA(1)=+0.5 -> +0.4x, ARMA(2,2) signs correct, all
+matching `arima`; R + Python parity intact.
+
 ### Documentation generation
 
 `man/` is generated from roxygen comments — never edit `.Rd` files by hand. Shared roxygen fragments live in `man-roxygen/` (e.g. `authors.R`, `keywords.R`) and are included via `@template authors`. `DESCRIPTION` sets `Roxygen: list(old_usage = TRUE)`, so generated usage sections use the legacy style — don't switch it.

R/RcppExports.R

@@ -1,8 +1,8 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-.UCompC <- function(commands, ys, us, models, hs, lambdas, outliers, tTests, criterions, periodss, rhoss, verboses, stepwises, p0s, armas, TVPs, seass, trendOptionss, seasonalOptionss, irregularOptionss, nsims, seeds, lambdaLowers, aEndIns, PEndIns, innVarIns, betaAugIns) {
-    .Call(`_muse_UCompC`, commands, ys, us, models, hs, lambdas, outliers, tTests, criterions, periodss, rhoss, verboses, stepwises, p0s, armas, TVPs, seass, trendOptionss, seasonalOptionss, irregularOptionss, nsims, seeds, lambdaLowers, aEndIns, PEndIns, innVarIns, betaAugIns)
+.UCompC <- function(commands, ys, us, models, hs, lambdas, outliers, tTests, criterions, periodss, rhoss, verboses, stepwises, p0s, armas, TVPs, seass, trendOptionss, seasonalOptionss, irregularOptionss, nsims, seeds, lambdaLowers, aEndIns, PEndIns, innVarIns, betaAugIns, compVarSmootheds) {
+    .Call(`_muse_UCompC`, commands, ys, us, models, hs, lambdas, outliers, tTests, criterions, periodss, rhoss, verboses, stepwises, p0s, armas, TVPs, seass, trendOptionss, seasonalOptionss, irregularOptionss, nsims, seeds, lambdaLowers, aEndIns, PEndIns, innVarIns, betaAugIns, compVarSmootheds)
 }
 
 .UCompARMAC <- function(ys, arOrders_, maOrders_, armaLags_, criterion_) {

R/pts-internals.R

@@ -159,7 +159,8 @@
                            verbose, armaIdent,
                            irregularOptions = "arma(0,0)",
                            outlier = 0,
-                           lambdaLower = -Inf){
+                           lambdaLower = -Inf,
+                           compVarSmoothed = FALSE){
     # u: NULL -> sentinel matrix; vector -> row matrix; otherwise pass through
     if (is.null(u)){
         u_mat <- matrix(0, 1, 2)
@@ -192,7 +193,8 @@
         trendOptions     = "rw/llt/srw/td",
         seasonalOptions  = "none/linear/equal",
         irregularOptions = irregularOptions,
-        lambdaLower      = as.numeric(lambdaLower)
+        lambdaLower      = as.numeric(lambdaLower),
+        compVarSmoothed  = isTRUE(compVarSmoothed)
     )
 }
 
@@ -212,7 +214,10 @@
            if (is.null(args$aEnd)) numeric(0) else as.numeric(args$aEnd),
            if (is.null(args$PEnd)) matrix(0, 0, 0) else as.matrix(args$PEnd),
            if (is.null(args$innVar)) -1 else as.numeric(args$innVar),
-           if (is.null(args$betaAug)) numeric(0) else as.numeric(args$betaAug))
+           if (is.null(args$betaAug)) numeric(0) else as.numeric(args$betaAug),
+           # compVarSmoothed: TRUE => full backward smoother for two-sided
+           # smoothed component variances (bands); FALSE => fast filtered.
+           isTRUE(args$compVarSmoothed))
 }
 
 # .pts_resolve_class: mirror adam's input-class resolution at
@@ -462,7 +467,7 @@
 .pts_fit <- function(y, u, model, lags, h, criterion, armaIdent, verbose,
                      ar = 0L, ma = 0L, armaLags = 1L, outlier = 0,
                      lambdaLower = -Inf, B = NULL, uFuture = NULL,
-                     biasadj = FALSE){
+                     biasadj = FALSE, compVarSmoothed = FALSE){
     # Flatten per-lag ar / ma vectors into the format pts_to_uc consumes:
     #   length-1 lags -> c(p, q)         (non-seasonal arma(p,q))
     #   length-2 lags -> c(p, q, P, Q, s) (sarma(p,q)(P,Q)_s)
@@ -479,7 +484,8 @@
                                .pts_arma_candidates(ar, ma, armaLags,
                                                     armaIdent),
                            outlier = outlier,
-                           lambdaLower = lambdaLower)
+                           lambdaLower = lambdaLower,
+                           compVarSmoothed = compVarSmoothed)
     # Override the default sentinel (-9999.9) when the caller supplied an
     # explicit starting vector via `B` (adam-style internal hatch, passed
     # through pts(... , B = ...) and used by the loss-surface experiment
@@ -559,6 +565,14 @@
     rawComp <- .pts_ts_comp(out$comp, out$m, y, h)
     colnames(rawComp) <- strsplit(out$compNames, "/")[[1]]
     comp <- .pts_build_comp(rawComp, v)
+    # Component state variances (filtered by default; smoothed P_{t|T} when
+    # compVarSmoothed=TRUE was requested).  Same layout as the raw component
+    # matrix, so reshape identically and label with the component names.
+    compV <- if (!is.null(out$compV) && length(out$compV) > 0L){
+        cv <- .pts_ts_comp(out$compV, out$m, y, h)
+        colnames(cv) <- strsplit(out$compNames, "/")[[1]]
+        cv
+    } else NULL
     # fitted on the original scale (back-transformed from comp[, "Fit"]).
     # residuals stay as the engine's BC-scale innovations -- they are the
     # white-noise sequence used by the validation table's diagnostics.
@@ -631,6 +645,7 @@
         yFor         = yFor,           # original scale (length h, possibly 0)
         v            = v,
         comp         = comp,           # BC scale, additive
+        compV        = compV,          # component state variances (BC scale) or NULL
         fitted       = fitted,         # original scale
         residuals    = residuals,      # BC scale (engine innovations)
         scale        = scale,          # MLE sigma on the BC scale

R/pts.R

@@ -99,6 +99,11 @@
 #' observations of \code{data} are withheld from estimation and returned in
 #' \code{$holdout} for later accuracy assessment.
 #' @param verbose logical: print intermediate optimisation output.
+#' @param componentVariance logical (default \code{FALSE}).  When \code{TRUE}
+#'   the fitted object gains a \code{$compV} matrix of two-sided \emph{smoothed}
+#'   state variances \eqn{P_{t|T}} (for component confidence bands), obtained by
+#'   running the full backward smoother.  When \code{FALSE} the cheaper filtered
+#'   variances are returned in \code{$compV} instead.
 #' @param ... advanced / undocumented passthroughs.  Supported keys:
 #' \itemize{
 #'   \item \code{B} - numeric vector of starting values for the optimiser
@@ -164,6 +169,7 @@ pts <- function(data,
                 h          = 0,
                 holdout    = FALSE,
                 verbose    = FALSE,
+                componentVariance = FALSE,
                 ...){
     cl  <- match.call()
     tic <- proc.time()
@@ -404,6 +410,7 @@ pts <- function(data,
                     B         = B,
                     uFuture   = u_held,   # future xreg for the auto-forecast
                     biasadj   = biasadj,
+                    compVarSmoothed = isTRUE(componentVariance),
                     verbose   = verbose)
     # When h > 0 we cache the engine's forecast (length h, original scale).
     # When h == 0 we still populate $forecast with a 1-period NA placeholder
@@ -519,6 +526,7 @@ pts <- function(data,
         fitted     = res$fitted,        # original scale (back-transformed)
         residuals  = res$residuals,     # BC scale (engine innovations)
         comp       = res$comp,          # pts-specific BC-scale additive decomposition
+        compV      = res$compV,         # component state variances (filtered, or smoothed if componentVariance=TRUE)
         states     = statesMat,         # adam-aligned structural state evolution
         ## --- forecast convenience cache (NULL when pts is called with h = 0) ---
         forecast     = cachedFor,

python/src/muse/core/pts.py

@@ -62,6 +62,7 @@ def __init__(
         h: int = 0,
         holdout: bool = False,
         verbose: bool = False,
+        component_variance: bool = False,
     ):
         self.model = model
         self.lags = lags
@@ -82,6 +83,9 @@ def __init__(
         self.h = int(h)
         self.holdout = bool(holdout)
         self.verbose = bool(verbose)
+        # When True, components carry SMOOTHED state variances P_{t|T} (for
+        # confidence bands); default False = fast filtered variances.
+        self.component_variance = bool(component_variance)
         self._fit = None  # populated by fit()
 
     # ---- fit ------------------------------------------------------------
@@ -187,13 +191,30 @@ def fit(self, y, X=None):
             arma_lags = [int(v) for v in np.atleast_1d(sel["lags"])]
 
         # outlier detection threshold (adam-style level -> two-sided z).
-        # The lambda screen above already pins lambda to a number before the
-        # engine runs, so R's joint-lambda + outlier workaround is moot here.
         from scipy.stats import norm
 
         outlier_z = 0.0 if self.outliers == "ignore" else float(norm.ppf((1 + self.level) / 2))
 
+        # Mirror R's joint-lambda + outlier workaround (pts.R): the engine's
+        # outlier-injection refit has a parameter-vector dimensionality bug when
+        # it is ALSO jointly estimating lambda -- the extra lambda slot doesn't
+        # survive the dummy injection and the BFGS aborts (4x1 vs 3x1).  A
+        # numeric lambda screen already pins lambda, but the default
+        # lambda_estim="likelihood" leaves the "Z" power for the engine, so when
+        # outliers are requested with an auto-lambda power we pin lambda first
+        # from a quick no-outlier fit, then rewrite the spec with that number.
         arma_tuple, irregular = translate.arma_spec(ar_v, ma_v, arma_lags)
+        nm_s = len(model_str)
+        if outlier_z > 0 and model_str[: nm_s - 2].lower() == "z":
+            model_uc0, lam0 = translate.pts_to_uc(model_str, arma_orders=arma_tuple)
+            pre = self._engine(
+                y, self._u, model_uc0, lam0, lags, criterion, irregular,
+                arma_ident=False, outlier=0.0, lambda_lower=lambda_lower,
+            )
+            lam_chosen = round(float(pre["lambda"]), 4)
+            model_str = _fmt_lambda_str(lam_chosen) + model_str[nm_s - 2 :]
+            arma_tuple, irregular = translate.arma_spec(ar_v, ma_v, arma_lags)
+
         model_uc, lam = translate.pts_to_uc(model_str, arma_orders=arma_tuple)
         out = self._engine(
             y,
@@ -206,6 +227,7 @@ def fit(self, y, X=None):
             arma_ident=False,
             outlier=outlier_z,
             lambda_lower=lambda_lower,
+            comp_var_smoothed=self.component_variance,
         )
         if out.get("model") == "error":
             raise RuntimeError("muse engine returned an error for this spec.")
@@ -267,6 +289,7 @@ def _engine(
         arma_ident,
         outlier=0.0,
         lambda_lower=-math.inf,
+        comp_var_smoothed=False,
     ):
         periods = lags / np.arange(1, max(1, lags // 2) + 1)
         rhos = np.ones_like(periods)
@@ -294,6 +317,7 @@ def _engine(
             1,
             0,
             float(lambda_lower),
+            compVarSmoothed=bool(comp_var_smoothed),
         )
 
     def _arma_candidate(self):
@@ -698,6 +722,17 @@ def _post_process(self, out, y_train, y_full, held, lam_in, lambda_screened):
         v = np.asarray(out["v"], dtype=float)
         self._comp, self._comp_names = _build_comp(comp, names, v)
 
+        # Component state variances (filtered, or smoothed if
+        # component_variance=True was requested).  Same (m x T) -> (T x m)
+        # reshape as comp; labelled by the component names.  None if absent.
+        cv_raw = out.get("compV", None)
+        if cv_raw is not None and np.asarray(cv_raw, dtype=float).size > 0:
+            self._comp_v = np.asarray(cv_raw, dtype=float).T  # (T, m)
+            self._comp_v_names = list(names)
+        else:
+            self._comp_v = None
+            self._comp_v_names = None
+
         ns = len(y_train)
         error = self._comp[:, 0]
         fit_bc = self._comp[:, 1]
@@ -828,6 +863,15 @@ def sigma(self):
     def comp(self):
         return self._comp, self._comp_names
 
+    @property
+    def comp_variance(self):
+        """Component state variances (T x m) and their names, or (None, None).
+
+        Filtered variances by default; smoothed P_{t|T} when the model was
+        constructed with component_variance=True.  Mirrors R's `$compV`.
+        """
+        return self._comp_v, self._comp_v_names
+
     @property
     def orders(self):
         return self._orders

python/tests/dump_outlier_reference.R

@@ -14,7 +14,7 @@ for (cs in cases) {
              outliers = "use", level = cs$level)
     out[[cs$name]] <- list(
         model     = m$model,
-        nParam    = m$nParam,
+        nParam    = nparam(m),
         logLik    = as.numeric(logLik(m)),
         coefNames = names(coef(m)),
         coef      = as.numeric(coef(m)),

python/tests/dump_pts_reference.R

@@ -53,7 +53,7 @@ for (sp in specs) {
         AICc     = as.numeric(greybox::AICc(m)),
         BICc     = as.numeric(greybox::BICc(m)),
         nobs     = nobs(m),
-        nParam   = m$nParam,
+        nParam   = nparam(m),
         sigma    = as.numeric(sigma(m)),
         fitted   = as.numeric(fitted(m)),
         residuals = as.numeric(residuals(m)),

src/ARMAmodel.h

@@ -208,23 +208,20 @@ void arToPacf(vec& PAR){
       / (1 - PAR(i) * PAR(i));
   }
 }
-// Returns stationary polynomial from an arbitrary one
+// Map unconstrained params -> STANDARD AR coefficients phi of a stationary
+// polynomial (1 - phi_1 B - ... ) y(t) = a(t)  (Box-Jenkins / pacfToAr convention).
 void polyStationary(vec& PAR){
-  // (1 + PAR(1) * B + PAR(2) *B^2 + ...) y(t) = a(t)
   vec limits(2);
   limits(0) = -0.98;
   limits(1) = 0.98;
   constrain(PAR, limits);
   pacfToAr(PAR);
-  PAR = -PAR;
 }
-// Inverse of polyStationary
+// Inverse of polyStationary: STANDARD AR coefficients -> unconstrained params.
 void invPolyStationary(vec& PAR){
-  // (1 + PAR(1) * B + PAR(2) *B^2 + ...) y(t) = a(t)
   mat limits(PAR.n_elem, 2);
   limits.col(0).fill(-0.98);
   limits.col(1).fill(0.98);
-  PAR = -PAR;
   arToPacf(PAR);
   unconstrain(PAR, limits);
 }
@@ -367,7 +364,7 @@ void armaMatrices(vec p, SSmatrix* model, void* userInputs){
     int arDeg = inp->arDeg > 0 ? inp->arDeg : (int)arma::sum(orders % lags);
     arma::vec col(arDeg, arma::fill::zeros);
     fillSARMAcoefs(p.rows(1, inp->ar), orders, lags, true, col);
-    model->T.submat(0, 0, arDeg - 1, 0) = -col;
+    model->T.submat(0, 0, arDeg - 1, 0) = col;   // AR (1 - phi B): T top = +phi
   }
   // MA.
   if (inp->ma > 0){
@@ -404,7 +401,7 @@ void armaMatricesTrue(vec p, SSmatrix* model, void* userInputs){
     int arDeg = inp->arDeg > 0 ? inp->arDeg : (int)arma::sum(orders % lags);
     arma::vec col(arDeg, arma::fill::zeros);
     fillSARMAcoefs(p.rows(1, inp->ar), orders, lags, false, col);
-    model->T.submat(0, 0, arDeg - 1, 0) = -col;
+    model->T.submat(0, 0, arDeg - 1, 0) = col;   // AR (1 - phi B): T top = +phi
   }
   if (inp->ma > 0){
     arma::ivec orders, lags;