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/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..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 @@ -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, 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 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/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 c0f405a..ec3e9e7 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 @@ -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; @@ -1819,11 +1810,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; @@ -3412,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(); @@ -3435,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/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/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); diff --git a/src/musecore.h b/src/musecore.h index 005e56d..f137d34 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; @@ -319,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; @@ -365,14 +352,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 +371,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; } @@ -587,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){ 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; } 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); }