From d5b657d99f9c04694cb17b45a8989ea49622abad Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Wed, 22 Jul 2026 08:14:20 +0000 Subject: [PATCH 1/5] fix(engine): correct BoxCox lambda selection (diffuse Jacobian + variance ridge) The joint/profile BoxCox lambda MLE was biased, overshooting toward a bound (iid N(1000,50) -> lambda=2 not 1; Nile -> lambda=2 + spurious global trend). Two coupled defects in the concentrated likelihood (llik, SSpace.h): 1. Diffuse observations received their full-scale BoxCox change-of-variable Jacobian (lambda-1)*sum log(y) but NOT the concentrated-variance scale term -0.5*nDiff*log(s2) (dropped earlier to dodge a +Inf logLik). The KF runs in concentrated units, so log|Finf,t| = log(s2) + log(Finf*,t); dropping log(s2) left the Jacobian dangling -- a monotone-in-lambda term growing with the number of diffuse states, biasing lambda. Restored it: this is the exact-diffuse concentrated BCnorm marginal likelihood. 2. Restoring log(s2) reintroduced the s2->0 singularity on degenerate fits (short over-parameterised series interpolating the data -> +Inf logLik). Replaced the former hard floor with a weakly-informative inverse-gamma prior on the concentrated variance: s2 = (SSR + nu0*s0sq)/(n + nu0), nu0=1, s0sq = 1e-4*Var(g(y)). Bounds s2 smoothly; healthy fits (SSR >> nu0*s0sq) untouched; AICc's k-penalty + the existing n-k-1>0 admissibility do the rest. Verified: lambda recovers truth on clean oracles (iid->1, mult->0, additive->1) and real data (Nile-> local level lambda=1, AirPassengers-> log lambda=0, UKgas-> ~log); short-series singularity resolved (logLik negative); selections stable across prior strength 1e-6..1e-4. Shared engine, so identical in R/Python. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh --- src/SSpace.h | 49 ++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 38 insertions(+), 11 deletions(-) diff --git a/src/SSpace.h b/src/SSpace.h index 4862a11..ca87dda 100644 --- a/src/SSpace.h +++ b/src/SSpace.h @@ -707,7 +707,31 @@ double llik(vec& p, void* opt_data){ data->logJac = logJac; } if (data->cLlik){ // Concentrated Likelihood (MLE) - data->innVariance = v2F(0, 0) / nFinite; + data->innVariance = v2F(0, 0) / nFinite; // plain MLE SSR / n + // Weakly-informative inverse-gamma prior on the concentrated variance: + // s2 = (SSR + nu0*s0sq) / (n + nu0) (penalised / MAP estimate). + // On a degenerate fit (a flexible model interpolating a short series) + // SSR -> 0, so the plain MLE s2 -> 0 and log(s2) -> -Inf, which the diffuse + // term -0.5*nDiff*log(s2) turns into a spurious +Inf log-likelihood (the + // "+986 positive logLik" singularity). The prior bounds s2 away from 0 + // smoothly (no hard cliff): as SSR -> 0, s2 -> nu0*s0sq/(n+nu0) > 0, while + // healthy fits (SSR >> nu0*s0sq) are left essentially untouched. The prior + // scale s0sq is data-relative (a small fraction of the transformed-data + // variance) so it tracks lambda / data magnitude -- s2 spans many orders of + // magnitude across lambda, so a fixed absolute prior would be meaningless. + // NOTE: the density below is evaluated explicitly at this s2 (not via the + // -n/2 log s2 MLE shortcut), so the reported logLik is the honest penalised + // likelihood at the ridged s2. The analytic gradient's concentrated-MLE + // first-order condition (dLL/ds2 = 0) is only mildly violated in the + // degenerate corner (ridge ~ 0 for healthy fits), where the model is being + // rejected anyway. + { + uvec finMask = find_finite(data->y); + double s0sq = 1e-4 * arma::var(data->y.elem(finMask)); // prior scale + double nu0 = 1.0; // prior strength: 1 pseudo-obs + if (std::isfinite(s0sq) && s0sq > 0.0) + data->innVariance = (v2F(0, 0) + nu0 * s0sq) / (nFinite + nu0); + } const double s2 = data->innVariance; // Objective = the BCnorm log-likelihood (the C++ analogue of adam()'s // loss = "likelihood": the concentrated scale s2 is plugged straight into @@ -727,21 +751,24 @@ double llik(vec& p, void* opt_data){ yr.elem(ndIdx), data->y.elem(ndIdx) - data->v.elem(ndIdx), sqrt(s2 * data->F.elem(ndIdx)), lam)); - // Diffuse observations are consumed estimating the diffuse initial state: - // they carry no data-fit term, only the log|Finf| determinant (== - // diffuseTerm) plus their BoxCox Jacobian. The determinant is NOT scaled - // by the concentrated variance s2 -- the exact-diffuse likelihood - // convention (Durbin-Koopman): only the n - d_t non-diffuse observations - // carry the s2 scale. The previous form added -0.5*nDiff*log(s2), which - // is harmless near s2 ~ 1 but, when the optimiser drives s2 to an extreme - // corner (a degenerate trend whose observation variance has collapsed), - // injected a large spurious +likelihood -- the +986 "positive logLik" on - // strongly seasonal series. + // Diffuse observations are consumed estimating the diffuse initial state. + // They carry the log|Finf| determinant (== diffuseTerm) AND their BoxCox + // Jacobian AND the concentrated-variance scale term -0.5*nDiff*log(s2). + // The s2 scale is REQUIRED: the KF runs in concentrated (ratio) units so + // Finf is scale-free, and each diffuse obs' determinant is + // log|Finf,t| = log(s2) + log(Finf*,t); dropping the log(s2) piece leaves + // the full-scale BoxCox Jacobian (lambda-1)*sum_{t<=d} log(y_t) dangling + // with nothing to balance it -- a monotone-in-lambda term that biases the + // BoxCox lambda MLE (worse the more diffuse states there are). This is the + // exact-diffuse concentrated BCnorm marginal likelihood. NOTE: as s2 -> 0 + // (a collapsed/degenerate variance) log(s2) -> -Inf; guarded by the s2 + // floor at its formation above if needed. uvec dIdx = finiteIdx.elem(find(maskFin < 0.5)); if (!dIdx.is_empty()){ double nDiff = (double)dIdx.n_elem; LL += sum(bcnormLogJac(yr.elem(dIdx), lam)) - 0.5 * nDiff * std::log(2.0 * datum::pi) + - 0.5 * nDiff * std::log(s2) - 0.5 * diffuseTerm; } llikValue(0, 0) = -2.0 * LL / nFinite - std::log(2.0 * datum::pi); From c4b6b1a8158d0e0f8e27e18f3b538b601546e774 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Wed, 22 Jul 2026 08:28:03 +0000 Subject: [PATCH 2/5] feat(components): expose component state variances (compV) for confidence bands The engine computed component state variances internally but never returned them ("dead compV"). Surface them through the R interface so component uncertainty bands (and the prediction-variance widening over missing blocks) can be plotted. - components() gains a smoothedVar flag. Default (false) keeps the fast state smoother -> filtered variances (cheap; enough for the widening-over-a-gap picture). When true, the full backward Nt recursion runs so inputs.P holds the two-sided smoothed variances P_{t|T} (O(m^3), gated behind the flag). - compV threaded MuseOutputs -> musecpp2R -> pts object as $compV (BC scale, columns Level/Slope/Seasonal/Irregular), reshaped like $comp. - New pts() argument componentVariance = FALSE; TRUE requests smoothed bands. - One appended .UCompC arg (compVarSmootheds); RcppExports regenerated. Verified: filtered/smoothed compV finite & non-negative; smoothed < filtered (more information); filtered Level variance widens over an artificial gap and re-narrows. Lambda fix from the previous commit unaffected. Note: R front-end only -- musecpp2py.cpp / pts.py not yet mirrored (engine change is backward-compatible: compVarSmoothed defaults false, so Python builds and runs unchanged, just without the compV output). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh --- R/RcppExports.R | 4 ++-- R/pts-internals.R | 25 ++++++++++++++++++++----- R/pts.R | 3 +++ src/PTSmodel.h | 18 ++++++++++++------ src/RcppExports.cpp | 9 +++++---- src/musecore.h | 10 ++++++++-- src/musecpp2R.cpp | 5 ++++- 7 files changed, 54 insertions(+), 20 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 761bdad..77189c9 100644 --- a/R/RcppExports.R +++ b/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_) { diff --git a/R/pts-internals.R b/R/pts-internals.R index 42b17cf..246c797 100644 --- a/R/pts-internals.R +++ b/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 diff --git a/R/pts.R b/R/pts.R index 2bea526..da87765 100644 --- a/R/pts.R +++ b/R/pts.R @@ -164,6 +164,7 @@ pts <- function(data, h = 0, holdout = FALSE, verbose = FALSE, + componentVariance = FALSE, ...){ cl <- match.call() tic <- proc.time() @@ -404,6 +405,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 +521,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, diff --git a/src/PTSmodel.h b/src/PTSmodel.h index c0f405a..cc533f2 100644 --- a/src/PTSmodel.h +++ b/src/PTSmodel.h @@ -139,7 +139,7 @@ class BSMclass : public SSmodel{ // Check whether re-estimation is necessary void checkModel(uvec); // Components - void components(); + void components(bool smoothedVar = false); // Covariance of parameters (inverse of hessian) mat parCov(vec&); // Finding true parameter values out of transformed parameters @@ -1819,11 +1819,17 @@ void BSMclass::estimOutlier(vec p0, bool VERBOSE){ // SSmodel::inputs.u = bestU; } // Components -void BSMclass::components(){ - // components consumes only smoothed states (inputs.a); inputs.P feeds - // the dead compV. Request the fast state smoother (skips the O(m^3) - // backward Nt recursion + smoothed-variance work). - SSmodel::inputs.stateOnly = true; +void BSMclass::components(bool smoothedVar){ + // components consumes the smoothed states (inputs.a) and, into compV, + // the state variances from inputs.P. + // smoothedVar = false (default): fast state smoother (stateOnly) -- + // skips the O(m^3) backward Nt recursion, so inputs.P holds the + // one-sided FILTERED state variances. Cheap; enough to show e.g. + // the prediction variance widening over a missing block. + // smoothedVar = true: run the FULL smoother so the Nt recursion fills + // inputs.P with the two-sided SMOOTHED variances P_{t|T} (used for + // component confidence bands). O(m^3), hence gated behind the flag. + SSmodel::inputs.stateOnly = !smoothedVar; SSmodel::smooth(true); SSmodel::inputs.stateOnly = false; int nCycles = sum(inputs.rhos < 0), k = SSmodel::inputs.u.n_rows; diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 44fe86e..2b1ad4a 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -12,8 +12,8 @@ Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif // UCompC -SEXP UCompC(SEXP commands, SEXP ys, SEXP us, SEXP models, SEXP hs, SEXP lambdas, SEXP outliers, SEXP tTests, SEXP criterions, SEXP periodss, SEXP rhoss, SEXP verboses, SEXP stepwises, SEXP p0s, SEXP armas, SEXP TVPs, SEXP seass, SEXP trendOptionss, SEXP seasonalOptionss, SEXP irregularOptionss, SEXP nsims, SEXP seeds, SEXP lambdaLowers, SEXP aEndIns, SEXP PEndIns, SEXP innVarIns, SEXP betaAugIns); -RcppExport SEXP _muse_UCompC(SEXP commandsSEXP, SEXP ysSEXP, SEXP usSEXP, SEXP modelsSEXP, SEXP hsSEXP, SEXP lambdasSEXP, SEXP outliersSEXP, SEXP tTestsSEXP, SEXP criterionsSEXP, SEXP periodssSEXP, SEXP rhossSEXP, SEXP verbosesSEXP, SEXP stepwisesSEXP, SEXP p0sSEXP, SEXP armasSEXP, SEXP TVPsSEXP, SEXP seassSEXP, SEXP trendOptionssSEXP, SEXP seasonalOptionssSEXP, SEXP irregularOptionssSEXP, SEXP nsimsSEXP, SEXP seedsSEXP, SEXP lambdaLowersSEXP, SEXP aEndInsSEXP, SEXP PEndInsSEXP, SEXP innVarInsSEXP, SEXP betaAugInsSEXP) { +SEXP UCompC(SEXP commands, SEXP ys, SEXP us, SEXP models, SEXP hs, SEXP lambdas, SEXP outliers, SEXP tTests, SEXP criterions, SEXP periodss, SEXP rhoss, SEXP verboses, SEXP stepwises, SEXP p0s, SEXP armas, SEXP TVPs, SEXP seass, SEXP trendOptionss, SEXP seasonalOptionss, SEXP irregularOptionss, SEXP nsims, SEXP seeds, SEXP lambdaLowers, SEXP aEndIns, SEXP PEndIns, SEXP innVarIns, SEXP betaAugIns, SEXP compVarSmootheds); +RcppExport SEXP _muse_UCompC(SEXP commandsSEXP, SEXP ysSEXP, SEXP usSEXP, SEXP modelsSEXP, SEXP hsSEXP, SEXP lambdasSEXP, SEXP outliersSEXP, SEXP tTestsSEXP, SEXP criterionsSEXP, SEXP periodssSEXP, SEXP rhossSEXP, SEXP verbosesSEXP, SEXP stepwisesSEXP, SEXP p0sSEXP, SEXP armasSEXP, SEXP TVPsSEXP, SEXP seassSEXP, SEXP trendOptionssSEXP, SEXP seasonalOptionssSEXP, SEXP irregularOptionssSEXP, SEXP nsimsSEXP, SEXP seedsSEXP, SEXP lambdaLowersSEXP, SEXP aEndInsSEXP, SEXP PEndInsSEXP, SEXP innVarInsSEXP, SEXP betaAugInsSEXP, SEXP compVarSmoothedsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -44,7 +44,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< SEXP >::type PEndIns(PEndInsSEXP); Rcpp::traits::input_parameter< SEXP >::type innVarIns(innVarInsSEXP); Rcpp::traits::input_parameter< SEXP >::type betaAugIns(betaAugInsSEXP); - rcpp_result_gen = Rcpp::wrap(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)); + Rcpp::traits::input_parameter< SEXP >::type compVarSmootheds(compVarSmoothedsSEXP); + rcpp_result_gen = Rcpp::wrap(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)); return rcpp_result_gen; END_RCPP } @@ -65,7 +66,7 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { - {"_muse_UCompC", (DL_FUNC) &_muse_UCompC, 27}, + {"_muse_UCompC", (DL_FUNC) &_muse_UCompC, 28}, {"_muse_UCompARMAC", (DL_FUNC) &_muse_UCompARMAC, 5}, {NULL, NULL, 0} }; diff --git a/src/musecore.h b/src/musecore.h index 005e56d..08ad8c3 100644 --- a/src/musecore.h +++ b/src/musecore.h @@ -54,6 +54,10 @@ struct MuseInputs { // R-side guard (set to 1e-10 when y has zeros). Plumbed straight // into SSinputs::lambdaLower. double lambdaLower = -arma::datum::inf; + // When true, components() runs the full backward smoother so compV holds + // the two-sided SMOOTHED state variances P_{t|T} (for component confidence + // bands). Default false = fast state smoother, filtered variances only. + bool compVarSmoothed = false; // Terminal-state cache (forecastOnly): when aEndIn is non-empty, the // engine reuses it (and PEndIn / innVarIn / betaAugIn) instead of // re-filtering the whole series -- so forecast.pts / predict() are O(h) @@ -117,6 +121,7 @@ struct MuseOutputs { // components / all bool hasComponents = false; arma::mat comp; + arma::mat compV; // state variances (filtered, or smoothed if requested) int m = 0; std::string compNames; @@ -365,14 +370,14 @@ inline void runMuseCommand(MuseInputs in, MuseOutputs& out){ // -- components (part of "all") -- if (command == "all"){ - sysBSM.components(); + sysBSM.components(in.compVarSmoothed); inputs2 = sysBSM.getInputs(); string compNames = inputs2.compNames; if (iniObs > 0){ inputs = sysBSM.SSmodel::getInputs(); uvec missing = find_nonfinite(inputs.y.rows(0, iniObs)); sysBSM.interpolate(iniObs); - sysBSM.components(); + sysBSM.components(in.compVarSmoothed); inputs2 = sysBSM.getInputs(); uvec rowI(1); rowI(0) = 0; if (compNames.find("Level") != string::npos) rowI++; @@ -384,6 +389,7 @@ inline void runMuseCommand(MuseInputs in, MuseOutputs& out){ } out.hasComponents = true; out.comp = inputs2.comp; + out.compV = inputs2.compV; out.m = static_cast(inputs2.comp.n_rows); out.compNames = compNames; } diff --git a/src/musecpp2R.cpp b/src/musecpp2R.cpp index 6a8a6d6..4779170 100644 --- a/src/musecpp2R.cpp +++ b/src/musecpp2R.cpp @@ -26,7 +26,8 @@ SEXP UCompC(SEXP commands, SEXP ys, SEXP us, SEXP models, SEXP hs, SEXP trendOptionss, SEXP seasonalOptionss, SEXP irregularOptionss, SEXP nsims, SEXP seeds, SEXP lambdaLowers, - SEXP aEndIns, SEXP PEndIns, SEXP innVarIns, SEXP betaAugIns){ + SEXP aEndIns, SEXP PEndIns, SEXP innVarIns, SEXP betaAugIns, + SEXP compVarSmootheds){ // --- Marshall SEXP -> MuseInputs (no engine logic in this file) --- MuseInputs in; @@ -58,6 +59,7 @@ SEXP UCompC(SEXP commands, SEXP ys, SEXP us, SEXP models, SEXP hs, in.nsim = as(nsims); in.seed = static_cast(as(seeds)); in.lambdaLower = as(lambdaLowers); // -Inf sentinel = no bound + in.compVarSmoothed = as(compVarSmootheds); // TRUE => smoothed compV bands // Terminal-state cache (empty aEnd => no cache => full re-filter). in.aEndIn = vec(aEndr.begin(), aEndr.size(), false); in.PEndIn = mat(PEndr.begin(), PEndr.nrow(), PEndr.ncol(), false); @@ -102,6 +104,7 @@ SEXP UCompC(SEXP commands, SEXP ys, SEXP us, SEXP models, SEXP hs, } if (out.hasComponents){ output("comp") = out.comp; + output("compV") = out.compV; output("m") = out.m; output("compNames") = out.compNames; } From 18f83bfbe01e48e29ff98213be361bd3e5888c92 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Wed, 22 Jul 2026 12:12:18 +0000 Subject: [PATCH 3/5] refactor(arma): standard (1-phi B)/(1+theta B) conventions internally polyStationary/invPolyStationary now return/consume the STANDARD Box-Jenkins AR coefficients (pacfToAr convention, no internal negation), so the engine's internal ARMA parameter equals the reported coefficient -- AR is (1 - phi B), MA is (1 + theta B), both matching stats::arima. Removes the scattered compensating flips that this convention mismatch required: - armaMatrices/armaMatricesTrue: T = +col (AR), R = +col (MA). - musecore.h reporting + parameterValues: no AR sign flip. - setEstimatedParams: no AR *= -1 on the true-param path. - AR seed: sample PACF maps straight into estimation space (no negate). This fixes the AR-on-irregular basin bug: the old seed fed a standard (1-phi B) PACF into the internal (1+phi B) invPolyStationary, stranding a correct +phi seed in the -phi basin. For pure-structural models the optimiser escaped it, but under the (correct) exact-diffuse log(s2) term that basin is a local minimum it cannot leave -- the RW level stays active and the AR flips sign. With the seed in the right basin the level collapses on stationary data and AR recovers +phi. MA seed: clamp the arToPacf input to the invertible region (|PACF| < 1) as the AR seed already does; otherwise a 2-block SARMA MA seed divides by (1 - phi^2) at the unit boundary and crashes invPolyStationary. Verified: AR(1)=0.7 -> +0.67, MA(1)=0.5 -> +0.4x, ARMA(2,2)/SARMA signs correct (match arima); all 96 R testthat pass; lambda oracles correct (iid->1, Nile-> local level, Air->log). CLAUDE.md updated with the full rationale. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh --- CLAUDE.md | 39 +++++++++++++++++++++++++++++++++++++++ src/ARMAmodel.h | 13 +++++-------- src/PTSmodel.h | 43 ++++++++++++++++++++++++++----------------- src/musecore.h | 26 ++++---------------------- 4 files changed, 74 insertions(+), 47 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 98b8136..f1ae89a 100644 --- a/CLAUDE.md +++ b/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. diff --git a/src/ARMAmodel.h b/src/ARMAmodel.h index 29362e6..ae45dd8 100644 --- a/src/ARMAmodel.h +++ b/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; diff --git a/src/PTSmodel.h b/src/PTSmodel.h index cc533f2..ec3e9e7 100644 --- a/src/PTSmodel.h +++ b/src/PTSmodel.h @@ -1072,19 +1072,10 @@ void BSMclass::estim(vec p, bool VERBOSE){ // variances or on irregular variances < 1e-6). The actual user-supplied // parameters are then handed to this method. void BSMclass::setEstimatedParams(vec userParams){ - // Incoming AR/SAR coefficients are in the standard (1 - phi B) - // convention (the reported convention -- see the output flip in - // runMuseCommand). Convert them back to the engine's internal - // (1 + phi B) form (polyStationary negates, armaMatrices builds - // T = -col) before building the system matrices. MA is already - // standard (armaMatrices uses R = +col), so it is left untouched. - if (inputs.ar > 0){ - arma::vec ncum = arma::cumsum(inputs.nPar); - arma::uword a0 = (arma::uword)ncum(2) + 1; - arma::uword a1 = a0 + (arma::uword)inputs.ar - 1; - if (a1 < userParams.n_elem) - userParams.rows(a0, a1) *= -1.0; - } + // Incoming AR (1 - phi B) and MA (1 + theta B) coefficients are the + // reported convention, which is now also the engine's internal one: + // armaMatricesTrue builds T = +col and R = +col directly, so no sign + // conversion is needed here. SSmodel::inputs.userModel = bsmMatricesTrue; SSmodel::inputs.cLlik = false; SSmodel::inputs.userInputs = &inputs; @@ -3418,17 +3409,20 @@ void BSMclass::initParBsm(){ } } // Converting to estimation space - // AR pars + // AR pars. beta0ARMA holds the standard (1 - phi B) AR + // coefficients (from the sample PACF); polyStationary/ + // invPolyStationary now use that same convention, so the seed + // maps straight into estimation space with no sign flip. if (inputs.ar > 0){ beta0aux = inputs.beta0ARMA(arma::span(0, inputs.ar - 1)); - beta0aux1 = -beta0aux; // Correction for non-stationary polynomial + beta0aux1 = beta0aux; arToPacf(beta0aux1); ind = find(abs(beta0aux1) >= 1); if (ind.n_elem > 0){ beta0aux1(ind) = sign(beta0aux1(ind)) * 0.96; pacfToAr(beta0aux1); - beta0aux = -beta0aux1; + beta0aux = beta0aux1; } invPolyStationary(beta0aux); beta0aux.elem(find_nonfinite(beta0aux)).zeros(); @@ -3441,7 +3435,22 @@ void BSMclass::initParBsm(){ // Bringing MA polynomial to invertibility maInvert(beta0aux); inputs.beta0ARMA(arma::span(inputs.ar, inputs.ar + inputs.ma - 1)) = beta0aux; - // Parameterising polynomial to be invertible + // Guard invPolyStationary the same way the AR seed does: + // maInvert makes the (1 + theta B) polynomial invertible, + // but arToPacf() (inside invPolyStationary) treats the coefs + // as a (1 - phi B) AR polynomial, whose partial-autocorr + // recursion divides by (1 - phi^2) and blows up at the unit + // boundary. Clamp any |PACF| >= 1 to 0.96 before mapping to + // estimation space (a bad SARMA seed then just starts a bit + // inside the invertible region instead of crashing). + beta0aux1 = beta0aux; + arToPacf(beta0aux1); + ind = find(abs(beta0aux1) >= 1); + if (ind.n_elem > 0){ + beta0aux1(ind) = sign(beta0aux1(ind)) * 0.96; + pacfToAr(beta0aux1); + beta0aux = beta0aux1; + } invPolyStationary(beta0aux); aux = regspace(ini + inputs.ar, 1, ini + inputs.ar + inputs.ma - 1); SSmodel::inputs.p0(aux) = beta0aux; diff --git a/src/musecore.h b/src/musecore.h index 08ad8c3..f137d34 100644 --- a/src/musecore.h +++ b/src/musecore.h @@ -324,27 +324,9 @@ inline void runMuseCommand(MuseInputs in, MuseOutputs& out){ out.v = inputs.v; out.covp = inputs.covp; out.coef = inputs.coef; - // Report AR/SAR coefficients in the standard (1 - phi B) convention. - // The engine stores them internally in the (1 + phi B) form - // (polyStationary negates, armaMatrices builds T = -col), so the - // reported coefficient is the negative of the standard one. Flip the - // AR/SAR slice [nparCum(2)+1 .. +ar] for output (MA is already standard - // -- armaMatrices uses R = +col). setEstimatedParams() flips it back - // when a fitted object's coefficients are fed to forecastOnly, so the - // round-trip is exact and forecasts are unaffected. - if (inputs2.ar > 0){ - vec ncum = cumsum(inputs2.nPar); - uword a0 = (uword)ncum(2) + 1; - uword a1 = a0 + (uword)inputs2.ar - 1; - if (a1 < out.coef.n_elem) - out.coef.rows(a0, a1) *= -1.0; - // Keep covp consistent: Cov(-AR_i, x) = -Cov(AR_i, x); the AR-AR - // block and AR variances are unchanged (negated on both row & col). - if (out.covp.n_rows > a1 && out.covp.n_cols > a1){ - out.covp.rows(a0, a1) *= -1.0; - out.covp.cols(a0, a1) *= -1.0; - } - } + // AR/SAR are stored internally in the standard (1 - phi B) convention + // (polyStationary now returns +phi) and MA in (1 + theta B), so + // inputs.coef is already in the reported convention -- no sign flip. // Outlier detection populated by estimOutlier() — empty matrix // when in.outlier == 0 or nothing was detected. out.typeOutliers = inputs2.typeOutliers; @@ -593,7 +575,7 @@ inline void runArmaScore(arma::vec y, arma::ivec arOrders, arma::ivec maOrders, if (pi == 0) continue; vec block = sysOut.p.rows(pIn, pIn + pi - 1); polyStationary(block); - for (int j = 0; j < pi; ++j) coef(pOut++) = -block(j); + for (int j = 0; j < pi; ++j) coef(pOut++) = block(j); // AR (1 - phi B) pIn += pi; } for (uword b = 0; b < maOrders.n_elem; ++b){ From d45803c3d51e2e79e1270197b0a99259eb5ac2a9 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Wed, 22 Jul 2026 12:12:30 +0000 Subject: [PATCH 4/5] python: mirror compV + outlier lambda-pin workaround; fix parity nParam refs - musecpp2py.cpp + core/pts.py: expose component state variances (compV) and the component_variance flag, mirroring the R $compV feature (filtered by default; smoothed P_{t|T} when requested). Same shared C++ engine, so identical output. - core/pts.py: mirror R's joint-lambda + outlier workaround -- when outliers are requested with an auto-lambda ("Z") power, pin lambda from a quick no-outlier fit first, sidestepping the engine's dummy-injection dimensionality bug. - tests/dump_{pts,outlier}_reference.R: dump nparam(m) (scalar total) instead of the 2x5 $nParam matrix, so the parity harness compares like-for-like. Verified: 45 Python functionality tests pass; all six R<->Python parity suites match to <= 1e-6 (engine, pts, forecast, summary, outliers, diagnostics). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh --- python/src/muse/core/pts.py | 48 +++++++++++++++++++++++++-- python/tests/dump_outlier_reference.R | 2 +- python/tests/dump_pts_reference.R | 2 +- src/python/musecpp2py.cpp | 8 +++-- 4 files changed, 54 insertions(+), 6 deletions(-) diff --git a/python/src/muse/core/pts.py b/python/src/muse/core/pts.py index b6555ea..46be172 100644 --- a/python/src/muse/core/pts.py +++ b/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 diff --git a/python/tests/dump_outlier_reference.R b/python/tests/dump_outlier_reference.R index 3ec14b3..2df2e89 100644 --- a/python/tests/dump_outlier_reference.R +++ b/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)), diff --git a/python/tests/dump_pts_reference.R b/python/tests/dump_pts_reference.R index 27e1161..6522d00 100644 --- a/python/tests/dump_pts_reference.R +++ b/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)), diff --git a/src/python/musecpp2py.cpp b/src/python/musecpp2py.cpp index ead22cf..ce1b127 100644 --- a/src/python/musecpp2py.cpp +++ b/src/python/musecpp2py.cpp @@ -94,7 +94,8 @@ static py::dict ucomp(std::string command, py::array_t aEndIn, py::array_t PEndIn, double innVarIn, - py::array_t betaAugIn) { + py::array_t betaAugIn, + bool compVarSmoothed) { MuseInputs in; in.command = command; in.y = np_to_vec(y); @@ -124,6 +125,7 @@ static py::dict ucomp(std::string command, in.PEndIn = np_to_mat(PEndIn); in.innVarIn = innVarIn; in.betaAugIn = np_to_vec(betaAugIn); + in.compVarSmoothed = compVarSmoothed; MuseOutputs out; runMuseCommand(in, out); @@ -161,6 +163,7 @@ static py::dict ucomp(std::string command, } if (out.hasComponents) { d["comp"] = mat_to_np(out.comp); + d["compV"] = mat_to_np(out.compV); d["m"] = out.m; d["compNames"] = out.compNames; } @@ -228,5 +231,6 @@ PYBIND11_MODULE(_musecore, mod) { py::arg("aEndIn") = py::array_t(), py::arg("PEndIn") = py::array_t(), py::arg("innVarIn") = -1.0, - py::arg("betaAugIn") = py::array_t()); + py::arg("betaAugIn") = py::array_t(), + py::arg("compVarSmoothed") = false); } From ca98a7b43a301f149bb03eaffc0fff899912ffd2 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Wed, 22 Jul 2026 15:48:26 +0000 Subject: [PATCH 5/5] docs: document pts() componentVariance argument (fixes R-CMD-check) The componentVariance argument added with the compV feature had no roxygen @param entry, so R CMD check flagged "Undocumented arguments in documentation object 'pts'" -> WARNING (CI failure). Add the @param and regenerate man/pts.Rd. Verified locally: R CMD check documentation sections all OK; vignette builds; Python ruff + mypy + pytest (45) all pass. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh --- R/pts.R | 5 +++++ man/pts.Rd | 8 +++++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/R/pts.R b/R/pts.R index da87765..77f79f4 100644 --- a/R/pts.R +++ b/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 diff --git a/man/pts.Rd b/man/pts.Rd index e7daa57..10bc8ab 100644 --- a/man/pts.Rd +++ b/man/pts.Rd @@ -9,7 +9,7 @@ pts(data, model = "ZZZ", lags = stats::frequency(data), orders = list(ar outliers = c("ignore", "use", "select"), level = 0.99, ic = c("AICc", "BICc", "BIC", "AIC"), lambda_estim = c("likelihood", "guerrero", "decomp-guerrero"), biasadj = FALSE, h = 0, holdout = FALSE, - verbose = FALSE, ...) + verbose = FALSE, componentVariance = FALSE, ...) } \arguments{ \item{data}{response series. Either a univariate \code{ts} / \code{zoo} / @@ -118,6 +118,12 @@ observations of \code{data} are withheld from estimation and returned in \item{verbose}{logical: print intermediate optimisation output.} +\item{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.} + \item{...}{advanced / undocumented passthroughs. Supported keys: \itemize{ \item \code{B} - numeric vector of starting values for the optimiser