Sparce matrices in the C++ code + fixes for the box-cox transform#12
Merged
Conversation
Add an arma::sp_mat overload of KFprediction and build a sparse view Tsp(system.T) once per llik, routing the per-timestep prediction through the O(m^2) sparse path (T is ~98% zeros at high lags). Storage unchanged; KFinit keeps dense T (its eig_gen/schur is sparse-hostile but only ~2% of runtime). Triple products split into explicit binary products. Measured on pts(taylor[1:672], "1NT", lags=336) (m=672): forward filter ~3.4s/pass vs ~14s/pass dense (~4x), ~2.3x overall under same-load conditions. Microbenchmark confirms sp_mat*dense T*Pt*T' is 3.6x faster at this size. Point estimates bit-identical across the golden battery; only ill-conditioned Hessian-derived vcov/confint drift (accepted tolerance-based reordering). Includes temporary MUSE_PROF-gated profiling scaffolding (ScopeTimers), to be stripped before merge. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Apply the Phase-1 sp_mat(view) trick to: - gradLlik backward T-products (rt=T'*rt, Nt=T'*Nt*T) - auxFilter forward KFprediction (smoother/disturbance forward pass) - KFinnovations forward KFprediction - forecast() h-step KFprediction plus a forward declaration of the sparse KFprediction overload so the earlier-defined forecast() can see it. Point estimates stay within tolerance across the golden battery (only ill-conditioned vcov/confint drift, accepted). Note: gradLlik's sparse-T gives little gain on its own -- its backward loop is dominated by the dense Lt'*Nt*Lt term (Lt = I - Kt*Z), which needs rank-1 algebra (next). Extends the MUSE_PROF profiling with aux(smoother)/hess timers. All TEMP profiling scaffolding is env-gated and stripped before merge. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The analytic-gradient backward loop was dominated by the dense triple product Lt'*Nt*Lt with Lt = I - Kt*Z (identity minus rank-1). Expand it in closed form: Lt'*Nt*Lt = Nt - w z' - z w' + (k'w) z z', w = Nt*Kt Lt'*rt = rt - z (Kt'rt) turning two O(m^3) dense products into O(m^2) rank updates (microbench: 5.6x on this kernel at m=672). Same for the diffuse branch (Lt = I-Kinft*Z). Combined with the already-sparse T'*Nt*T (3.7x), the gradient backward pass is now fully O(m^2). Point estimates bit-identical across the golden battery (coef unchanged -> gradient still converges identically); only ill-conditioned vcov/confint drift. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The smoother backward pass spends most of its time on the Nt recursion and Pt*Nt*Pt (smoothed-variance) computation -- irreducibly O(m^3) and the dominant cost at high lags. But these outputs are DEAD: data.P feeds only compV (whose sole consumer is a commented-out block, never surfaced), and the outlier-mode rNrOut/rOut/NOut are written but never read anywhere. components() -- the dominant smoother call -- consumes only smoothed STATES (data.a), which depend on the rt recursion, not Nt. Add a stateOnly flag (SSinputs) gating the entire backward Nt recursion, Pt*Nt*Pt, and the disturbance/outlier blocks. components() sets it, turning its smoother from O(m^3 n) to O(m^2 n). Other callers (disturbance, outlier detection) keep the full path unchanged. Verified: smoothed states/comp bit-identical across the golden battery (only ill-conditioned vcov/confint drift); full testthat suite 988 PASS incl. outliers="use". Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…URE) Remove the MUSE_PROF ScopeTimer instrumentation used during the sparse refactor. Add a NEWS entry (high-lag performance) and an ARCHITECTURE.md section documenting the sparsity exploitation (sparse-T views, rank-1 gradient, fast state smoother) and its numerical transparency. Verified end to end: R 988 tests pass; Python 39 functionality + 21/21 R<->Python parity (worst 3.5e-9); golden battery point estimates bit-identical (only ill-conditioned vcov/confint drift, accepted). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Bring main's decoupled fit/forecast cache (terminal-state caching + xreg auto-forecast, R + Python) onto the sparse-matrix branch. Only NEWS conflicted (both branches added a top entry); resolved by keeping both -- main's "Decoupled fit and forecast" entry plus this branch's "Performance: high-lag models" sparse entry. All code files (PTSmodel.h, musecore.h, RcppExports, pts.R/pts-internals.R, pts.py, forecaster.py) auto-merged. Merged engine builds; R 988 tests pass; forecast cache works on the sparse engine (pts()$forecast == forecast.pts()).
The parameter covariance was computed from a one-sided forward-difference Hessian with a fixed 1e-5 step. That scheme is biased to O(h) and cancels catastrophically in weakly-identified directions: on ill-conditioned models it returned non-PSD "covariances" (implied correlations outside +/-1) and standard errors that swung by tens of percent under the ~1e-10 dense/sparse filter reordering. A. hessLlik (SSpace.h): central second differences with a per-parameter, magnitude-scaled step eps^(1/4)*max(|p|,1). Unbiased to O(h^2) and no longer amplifies the rounding floor in flat directions. B. parCov (PTSmodel.h): when the observed-information Hessian is indefinite or numerically singular at the optimum (boundary / weakly-identified variance, ARMA phi~theta near-cancellation), fall back to the OPG / BHHH estimator sum_t s_t s_t' from per-observation scores, which is PSD by construction. Genuinely PD Hessians keep pinv(H) unchanged. OPG matches the Hessian to machine precision as n grows (information-matrix equality); it is unavailable for augmented/xreg models (no per-obs innovations), which keep the improved Hessian. Effect: all previously-indefinite cases now return valid PSD covariances, and the dense/sparse vcov divergence collapses from ~100-144% to ~1e-7. Point estimates, logLik and ICs are unchanged; R (986) and Python (39) suites pass. Shared C++ engine, so R and Python both get the fix. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
parCov previously gated the whole covariance computation on `hess.is_finite()`, so when the second-difference Hessian evaluated to NaN/Inf at the optimum (e.g. a near-degenerate fit) the OPG fallback never ran and `vcov` was all-NaN. OPG is built from single, small per-observation perturbations and never inverts a near-singular matrix, so it can be finite where the Hessian is not. Restructure so OPG is attempted whenever the Hessian is indefinite, ill-conditioned, OR non-finite (non-augmented models); the Hessian path is unchanged for finite, well-conditioned cases. This does not rescue a *fully* degenerate fit: when the disturbance variances collapse to ~0 (a near-deterministic boundary) the Kalman innovation variances F_t go non-positive, so every per-observation score is non-finite too -- neither estimator is defined and vcov correctly stays NaN. Behaviour on all well-conditioned and ill-conditioned cases is unchanged (988 R tests, 39 Python tests pass). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The quasi-Newton optimiser was fed a wrong gradient for pure structural
models, so it stalled after ~2 iterations ("Q-Newton: Unable to decrease")
at a non-stationary point with near-zero variances -- the "estimation
degeneracy" behind the NaN vcov reports (e.g. pts(taylor, "0NT")).
Two independent defects, both fixed:
1. gradLlik (analytic disturbance-smoother gradient, SSpace.h): the
baseline Q/H (sysmatQ) was left on the stale RATIO scale from the last
objFun() call while the perturbed Qt was built on the ABSOLUTE scale
(xUncon), so dQt = (Qt - sysmatQ)/inc mixed units and blew up by
~1/innVariance when the concentrated variance was small (log /
near-deterministic data -> gradients ~1e10). Rebuild Q/H at p before
sysmatQ. Also normalise by nFinite = n - nMiss (what llik() averages
over), not n - nMiss - d_t - 1 (a few-percent error on seasonal models).
The analytic gradient now matches a central-difference reference to ~6
digits across structural specs.
2. quasiNewtonBSM (PTSmodel.h): descent-direction safeguard -- when the
BFGS inverse-Hessian loses positive-definiteness and yields d with
d.grad >= 0 (not a descent direction), reset iHess = I (steepest
descent) so the line search can make progress instead of stalling.
Together: every structural spec on AirPassengers improves, none regress
(PTS(0,N,T) +129, PTS(1,N,T) +94, PTS(1,L,T) +15 log-lik); the taylor
"0NT" case converges (logLik -14580 vs -18545) instead of giving up.
ARMA/damped/cycle/regressor models already use the numerical gradient
(unaffected by #1) and also benefit from #2. Shared C++ engine, so R and
Python are identical. R suite green (988); one order-selection assertion
relaxed to allow the now-better fit to drop the seasonal ARMA block.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The Python functional test for outliers="use" was failing because two real bugs left an injected spike undetected on a 1LT fit: 1. Disturbance smoother (SSpace.h, auxFilter smooth==2): the per-step standardisation divided by sqrt(diag(Veta)), which is SINGULAR when an estimated disturbance variance collapses (0/0 -> nan), and the post-loop re-standardisation used one matrix pinv() across all components -- its single tolerance treated any component whose variance was far below the largest (e.g. a near-zero observation variance) as singular and zeroed it, destroying the outlier signal there. Now emit the raw smoothed disturbance (Q R' r_t) and standardise EACH component by its own empirical SD (scale-robust). The spike's smoothed observation disturbance now surfaces (eps t-stat ~7.4 at the spike). 2. Outlier acceptance (PTSmodel.h, estimOutlier): a detected, highly significant outlier (backward-deletion t-stat ~47) was dropped by the final "IC worse than no-outlier baseline -> revert" gate. When a flexible model overfits a short series (the unbounded variance->0 likelihood singularity), the baseline IC is artificially good and no genuine outlier model beats it. Acceptance now relies on the per-outlier significance test already enforced by backward deletion; the outlier model is kept unless it fails to converge. Also records the testing-discipline rule in CLAUDE.md (no "pre-existing issues" -- a red test means it was broken; investigate and fix; always run R tests + pytest + ruff + mypy). R suite 988 pass, Python 39 pass, ruff clean, mypy clean. Shared C++ engine, so R and Python detect identically. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… filter With the corrected optimiser, a flexible model on some series drives a variance log-parameter very negative, so exp(2*p) underflows to EXACTLY 0 in bsmMatrices/bsmMatricesTrue. A zero disturbance variance makes the Kalman innovation variance F_t = Z P Z' + H zero, the gain K = P Z'/F_t becomes 0/0 = NaN, the terminal state is NaN, and the forecast explodes (the reported case: ZZZ on a eurostat-tourism series gave max|mean| ~4e12 with NaN terminal states, data maxing near 1.8e6 -- now ~3.4e4, finite, within range). Floor the variance log-params (typePar == 0) at var >= exp(-23) ~ 1e-10 at the top of both matrix builders -- numerically "zero" for the model but keeping the Kalman recursion well-defined. Structural zeros (companion-state Q rows) are assigned afterwards and stay 0; well-identified models (variances far above the floor) are unaffected. R suite 988 pass, Python 39 pass, ruff + mypy clean. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ive variances) Commit 088dee9 added a variance floor to both matrix builders using the same `clamp(p, -11.5, inf)`. That is correct for bsmMatrices, where p holds log-parameters and Q = exp(2*p), so the floor maps to var >= exp(-23) ~ 1e-10. But bsmMatricesTrue takes p as the ABSOLUTE variances directly (Q(i,i) = p(i)), so clamping to [-11.5, inf) permitted NEGATIVE variances down to -11.5 in the absolute-scale (vcov/Hessian) filter. Floor to a positive variance there (exp(2*-11.5)) instead. This is a partial fix: it stops negative variances entering the matrices, but does not by itself cure the negative displayed variance, which comes from the concentrated scale innVariance going negative (a separate non-PSD-filter root, under investigation). 988 R tests pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The augmented Kalman likelihood concentrates the innovation variance as (v2F - snBeta)/n, the residual sum of squares after projecting out the augmented states (diffuse initial states + regressors). That RSS is a sum of squares, hence >= 0 -- but when the augmented states nearly interpolate the data (a degenerate near-perfect fit, which the corrected optimiser can now reach), v2F and snBeta are both huge and nearly equal, so the subtraction loses every significant digit (catastrophic cancellation) and can flip negative. That made innVariance, the displayed variances (ratio * innVariance) and the parameter covariance non-positive (negative B, NaN standard errors). Locate the bug rather than discard the model: the matrix builders and the cLlik likelihood branch are correct; the fault is purely the cancellation in this one subtraction. main (no sparse refactor) produces sane positive variances on the same data, confirming this is a regression exposed by the corrected optimiser reaching the degenerate corner -- not a modelling error. Clamp the RSS to an eps-scaled non-negative floor; it only bites in the cancellation regime, well-identified fits (RSS >> eps*v2F) are untouched and keep their exact value. Verified: R suite FAIL 0 PASS 988; Python 39 passed, ruff + mypy clean. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Replace the previous clamp with the real fix. The concentrated innovation
variance is RSS/n. Forming RSS as the projection identity
v2F - sn'Sn^{-1}sn subtracts two large nearly-equal numbers at a
near-interpolating fit (the classic unstable way to compute a sum of
squares), loses all precision, and could come out negative -- the previous
commit only clamped the sign of that artifact.
Now form the residual explicitly: cache the innovation v_t, the augmented
sensitivity V_t and the weight 1/F_t during the forward pass, then sum the
squared GLS residuals RSS = Sum (v_t - V_t·beta_hat)^2 / F_t after beta_hat
is known. This is a sum of non-negative terms, structurally >= 0 in
floating point, so the negative-variance artifact cannot occur however
degenerate the fit -- the cause is removed, not masked. A tiny eps-scaled
floor remains solely as a log-domain guard against RSS == 0 at an exactly
interpolating fit.
The residual degeneracy (RSS genuinely -> 0, infinite likelihood, overfit
forecast) is unchanged and remains a modelling matter for the
disturbance-variance control, separate from this numerical fix.
Verified: R suite FAIL 0 PASS 988 (incl. R<->Python parity via the
reference dumps); Python 39 passed, ruff + mypy clean; the previously
negative ZLT fit now reports non-negative variances and a finite,
non-NaN parameter covariance.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Make the BCnorm density the single source for the concentrated likelihood in
the non-augmented path, used by both estimation and reporting. Previously
llik() built objFunValue from a hand-written closed form that was only the
"closed-form equivalent of Sum bcnormLogDensityScalar"; the actual bcnorm
density function was never called. Now the objective IS that sum, and the
reported logLik / information criteria derive from the same objFunValue
(computeLLIK), so bcnorm is the one function behind estimation and reporting.
The data likelihood is evaluated in ONE vectorised bcnormLogDensity() call over
the non-diffuse observations, each with its own predictive sd sqrt(innVar*F_t)
(heteroscedastic), then summed -- the C++ analogue of adam()'s
loss = "likelihood" (concentrate the scale, plug it straight into the density).
The d_t diffuse observations are consumed estimating the diffuse initial state
and contribute only the (scale-scaled) log|Finf| determinant plus their BoxCox
Jacobian, computed from the existing accumulator.
To support this:
- bcnormLogDensity() now takes a per-observation sigma vector (was scalar),
and a vectorised bcnormLogJac() overload is added.
- bcnormLogDensityScalar()'s positivity guard is made lambda-aware: y <= 0 is
only rejected for lambda != 1. At lambda == 1 (identity) the BCnorm
distribution is a plain normal valid for all real y -- this matches
greybox::dbcnorm and is required for raw-scale models whose data can be
negative (a porting bug that otherwise turned such a logLik into NaN).
Value-equivalent: R suite FAIL 0 PASS 988 (all logLik/AIC value tests, incl.
the ARMA AR(1) recovery, unchanged); Python 39 passed; ruff + mypy clean.
The augmented path (llikAug) still uses the closed form -- a follow-up.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Convert the augmented (regressor / diffuse) path to the same bcnorm-based
objective as llik(), so the BCnorm density is the single source for the
likelihood across both estimation paths and, via computeLLIK(), for the reported
logLik / information criteria.
The augmented log-likelihood is exactly
LL = Sum_t bcnormLogDensity(y_t, g(y_t) - r_t, sqrt(innVar * F_t), lambda)
- 0.5 * log|Sn|
where r_t = v_t - V_t*beta_hat is the GLS residual (already formed
cancellation-free) and -0.5*log|Sn| is the augmented-state determinant -- the
diffuse initial states + regressors integrated under a flat prior, the augmented
analogue of llik()'s log|Finf| diffuse correction. Unlike llik() there is no
per-observation diffuse split: the augmentation handles the initial states
globally through Sn, so every non-missing observation is an ordinary BCnorm data
term and the determinant is a single scalar correction.
Reuses the existing vectorised bcnormLogDensity() (per-observation sigma); the
GLS residual and RSS are now also formed vectorised. No new helper functions.
Value-equivalent: a regressor fit reproduces the previous logLik to 8 decimals
(-177.57855252 either way). R suite FAIL 0 PASS 988; Python 39 passed; ruff +
mypy clean.
With this, both llik() and llikAug() compute the likelihood through bcnorm.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ance
The concentrated BCnorm likelihood scaled the diffuse-initialisation
determinant by sigma^2 (the d_t diffuse observations entered the
-0.5*nDiff*log(sigma^2) term). That is the standard concentrated convention
near sigma^2 ~ 1, but it is not scale-invariant: the diffuse innovation
variances Finf are the diffuse PRIOR propagated by T and do NOT rescale with
the variance parameters, so sigma^2 * Finf is not invariant. When the
optimiser drives sigma^2 to an extreme corner (a degenerate local-linear-trend
whose observation variance collapses, sigma^2 ~ 1e-117), the diffuse term
injected ~ +1600 of spurious likelihood: the notorious +986 "positive logLik"
on strongly seasonal series, with the model then selected for an explosive
exp-trend forecast (~3.6e7 vs a 1.8e6 actual peak).
Fix: leave the diffuse determinant unscaled by sigma^2 -- the exact-diffuse
likelihood (Durbin-Koopman), where only the n - d_t non-diffuse observations
carry the concentrated scale. innVariance itself is unchanged (still v2F/n,
so de-concentration / B / vcov are untouched); only the objective's diffuse
term changes. The non-diffuse data terms already used the scale-invariant
product sigma^2 * F_t (via bcnorm), so this completes the scale-correctness.
Result: eurostat-323 ZLT logLik +986 -> -737, forecast 3.6e7 -> 7.9e5, and
ZZZ auto-selection moves from the exploding PTS(0,L,T) to a sensible seasonal
PTS(0,N,D).
KNOWN REGRESSION (under investigation, next): the rw+arma(1,0) AR(1) recovery
test (test_pts.R:350) now lands at AR=-0.67 instead of +0.67. This is NOT a
gradient artifact (numeric gradients give the same), so the corrected diffuse
likelihood genuinely shifts which optimum that model prefers between
{RW-var->0, AR=+0.67} and {RW-var>0, AR=-0.67}. Whether -0.67 is the correct
fit under the proper diffuse likelihood is being investigated against
stats::arima before resolving the test.
R suite FAIL 1 PASS 987 (only the AR(1) test above); Python 39 passed.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
muse stores ARMA coefficients internally in the (1 + phi B) convention:
polyStationary() negates (PAR = -PAR) and armaMatrices() builds the companion
as T = -col, which together produce the correct (1 - phi B) state-space
transition AND make the seasonal convolution in fillSARMAcoefs() work out (the
quadratic cross-terms rely on it). But parameterValues() reported the raw
internal block, so the reported AR/SAR coefficient was the NEGATIVE of the
standard Box-Jenkins / stats::arima value -- a +0.7 AR(1) signal was reported
as -0.7. (MA was already standard: armaMatrices uses R = +col.)
This was masked until the diffuse-determinant fix: with the old (scaled)
diffuse likelihood the optimiser converged to the wrong-sign effective AR
(residual ACF ~ 0.6, i.e. a failed fit) which the reporting negation flipped
back to a right-looking number. With the likelihood corrected the model now
fits (white residuals, correct forecast), exposing the reporting sign.
Fix (reporting only -- the internal machinery is correct and untouched):
- runMuseCommand(): flip the AR/SAR slice [nparCum(2)+1 .. +ar] of the
reported coef, and the corresponding rows/cols of its covariance, to the
standard sign.
- setEstimatedParams(): flip that slice back to internal before building the
forecast matrices, so feeding a fitted object's coefficients to
forecastOnly round-trips exactly -- forecasts are unchanged.
Done in the shared C++ engine, so R and Python stay in parity automatically.
Verified: rw+arma(1,0) on a +0.7 AR(1) signal now reports AR(1) = +0.67
(matches stats::arima), forecast still decays correctly, MA still +; the
ARMA AR(1) recovery test passes again. R suite FAIL 0 PASS 988; Python 39
passed; ruff + mypy clean; eurostat-323 ZLT logLik still -737.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Zero-heavy / intermittent series were forced onto the raw scale (lambda=1) and routinely forecast negative values. Two coupled causes, both fixed: 1. BCnorm density guard too aggressive (regression from the bcnorm conversion). bcnormLogDensityScalar rejected y<=0 at every lambda!=1, so any user-supplied lambda!=1 on data with zeros gave a NaN log-likelihood. But g(y)=(y^lambda-1) /lambda is finite at y==0 for lambda>0 (g(0)=-1/lambda, e.g. sqrt(0)=0); it is undefined only for y<=0 with lambda<=0 (log 0 / 0^negative), and y<0 with fractional lambda gives a complex g already caught by the !isfinite(g) check. Reject only that genuinely-undefined corner now -- so any lambda>0 fits zero-containing data. 2. The msdecompose + Guerrero lambda screen bailed to 1 on any(y<=0). It now disqualifies only on NEGATIVE values (y^lambda is complex for y<0). Because the CV only sees positive block levels it would otherwise drive lambda->0 (mapping every zero to -Inf), so a floor lambda>=log(2)/log(max(y)) is applied when zeros are present -- the transformed zero is then no more extreme than the transformed maximum. Capped at 1 for tiny-count series. A power transform in (0,1) fits intermittent series much better AND guarantees non-negative forecasts (the inverse BoxCox cannot go below zero). Example: eurostat-267 moves from PTS(1,G,D) forecasting down to -101, to PTS(0.13,N,D) forecasting in [0,15] (holdout) with logLik -254 vs -430. Density fix is in the shared C++ engine; the R (.pts_guerrero_decomp_lambda) and Python (lambda_screen.guerrero_decomp_lambda) screens are updated identically and verified to return the same lambda (0.0782 on series 267). R suite FAIL 0 PASS 988; Python 39 passed; ruff + mypy clean. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add doc/guerrero.md describing how pts() chooses the Box-Cox lambda when the power slot is "Z": the variance-stabilisation rationale, the Guerrero CV objective and its delta-method derivation, the msdecompose-based block recipe, the clipped [0,2] range, the fallback conditions, and -- in full -- the zero-series handling. The zero floor lambda >= log(2)/log(max(y)) is derived from scratch (the requirement |g(0)| <= |g(max)| solved for lambda) and is flagged explicitly as a muse-specific heuristic with NO literature reference; the doc cites Guerrero (1993) for the screen itself and points to Yeo-Johnson (2000) / shifted Box-Cox as the textbook ways to admit zeros if a pedigreed treatment is wanted later. doc/ is a non-standard top-level directory, so it is added to .Rbuildignore. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…recast
Back-transforming the Box-Cox-scale point forecast g^{-1}(mu) returns the
conditional MEDIAN on the original scale; for lambda < 1 the inverse transform
is convex, so the median lies below the mean and the reported point forecast
(labelled `mean`) systematically under-forecasts -- ~8% low on eurostat-325
(lambda = 0.43), more as lambda -> 0.
Apply the standard second-order (delta-method) bias correction,
mean ~= g^{-1}(mu) * (1 + 0.5 * var * (1 - lambda) / (1 + lambda*mu)^2),
with var the BC-scale forecast variance -- the same correction as
forecast::InvBoxCox(biasadj = TRUE). A new helper .inv_box_cox_mean()
(R) / inv_box_cox_mean() (Python) does the median back-transform plus this
factor; it is used for the in-fit forecast and forecast()$mean. The
correction is guarded at the support boundary (1 + lambda*mu <= 0) and falls
back to no correction there, so heavy-tailed intermittent forecasts stay
bounded rather than blowing up.
Prediction-interval quantiles are NOT bias-corrected -- they are exact
quantiles of the back-transformed distribution. At lambda = 1 nothing changes.
Result: eurostat-325 forecast sum / actual moves 0.92 -> 0.99 (bias removed),
RMSSE 0.597 -> 0.577. Done identically in R and Python (verified: both return
72.8722 for mean(mu=10, var=50, lambda=0.43)).
R suite FAIL 0 PASS 988; Python 39 passed; ruff + mypy clean.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Information criteria now count the estimated diffuse structural initials (initial level, slope for L/G/D trend, cycle, and seasonal states) as degrees of freedom. These are profiled out / diffuse and so are genuine parameters; leaving them free let BICc/AICc under-penalise flexible seasonal and trend shapes. The count ns(0)+ns(1)+ns(2) is read from the engine's own state dimensions, so it is lags-driven (correct for multi-seasonal lags=c(m1,m2,...)) with no hard-coded period sizes; the stationary ARMA block ns(3) is excluded (its initial is the stationary distribution, not free -- the ARMA coefficients stay counted). The engine adds the same quantity (plus regressors) to its own selection k (kFor in BSMclass::estim), so pts(model="ZZZ") selection and the reported nParam agree. $nParam is now a 2x5 matrix mirroring smooth::adam: rows Estimated / Provided, columns nParamInternal, nParamXreg, nParamOccurrence, nParamScale, nParamAll. nparam() reads [Estimated, nParamAll] (the total DoF, unchanged in value). The initials fold into nParamInternal exactly as adam does -- no separate nInitial slot; the concentrated variance is nParamScale=1; regressors are nParamXreg (this also makes the reported count match the engine, which already counted u.n_rows). coef/vcov/ confint still cover only the covariance-bearing coefficients, so length(coef) is now smaller than nparam(). Shared C++ engine, so identical in R and Python (new n_param_table DataFrame mirrors R's matrix). Updated affected tests, ARCHITECTURE.md, man/, and both NEWS files. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Three related changes to Box-Cox lambda handling and forecasting, in both R and Python (shared C++ engine; results verified equivalent). 1. Zero-aware lambda lower bound (fixes likelihood-based lambda estimation on zero-heavy series). Box-Cox maps y = 0 to g(0) = -1/lambda, which at lambda <= 0 is -Inf and silently drops the observation -- so the likelihood (and AICc) at small lambda was computed on fewer points and was not comparable across lambda, luring a likelihood search toward a degenerate small lambda whose back-transform explodes. The joint-lambda lower bound is now the proper zero floor log(2)/log(max(y)) (was a flat 1e-10), shared by the screens and the engine via .pts_lambda_zero_floor / lambda_zero_floor. On an intermittent tourism series the auto lambda moved ~0.09 -> ~0.25 and the forecast ~2e5 -> ~8e3 (actuals ~7e3). 2. New `lambda_estim` argument: "likelihood" (new default -- joint ML in the engine), "guerrero" (classical Guerrero 1993 on raw season blocks), or "decomp-guerrero" (Guerrero on an msdecompose trend, the former default). The Guerrero screens ignore the likelihood and can pick a poor lambda; the likelihood estimator avoids this. Screens stay available, share the CV optimiser (.pts_guerrero_cv_lambda / _guerrero_cv). 3. New `biasadj` argument to pts() / forecast() / predict() (default FALSE). Point forecasts are now the back-transformed conditional MEDIAN by default; biasadj = TRUE applies the second-order bias correction for the conditional MEAN. The mean is not robust under a small-lambda Box-Cox (heavy right tail), so the median is the safer default. Interval quantiles unaffected. Also: default information criterion switched to AICc (was BICc), matching smooth::adam. Updated tests (R + Python), ARCHITECTURE.md, man/, and both NEWS files. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ollapse) When lambda is jointly estimated, llik() clamped the trial lambda only on the lower side; with no upper clamp the BFGS could push p.back() past LAMBDA_BOUND_UPPER = 2 into a region where the leading variance and Jacobian terms cancel and the BCnorm log-likelihood drifts to ~0, while the *reported* lambda was snapped back to 2. That left the objective inconsistent with the reported lambda (e.g. estimated "ZND" reported logLik ~ 0 at lambda = 2, vs -723 for a genuine fixed lambda = 2), which in turn corrupted structure selection -- a strongly seasonal series (eurostat 162) was fit as a no-seasonal damped model that missed the seasonal peak entirely. Clamp the joint-estimated lambda to [lambdaLower, LAMBDA_BOUND_UPPER] so the objective at the boundary matches the reported lambda. estimated == fixed at lambda = 2 now (logLik -723.39 vs -723.30), and series 162 selects the seasonal PTS(2,G,D) (RMSSE 1.28 -> 0.17). Shared engine, so identical in R and Python; both report lambda = 2.0000, logLik = -723.3934. Note: this fixes the collapse, not the deeper degeneracy -- joint Box-Cox MLE with a flexible model still drives lambda to a bound (the flexible seasonal absorbs the transform, so the Jacobian sets lambda). The default-estimator question is tracked separately. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The joint-lambda BFGS optimised lambda directly against a hard clamp to its [lower, upper] bounds (lower = max(LAMBDA_BOUND_LOWER, lambdaLower), upper = LAMBDA_BOUND_UPPER). The clamp made the objective FLAT beyond a bound, so the finite-difference gradient there was exactly 0 and the optimiser got trapped at the bound: it pinned lambda at the floor or cap and missed genuine interior optima. On a 74%-zero series (eurostat 267) it stuck at the zero floor 0.078 -- where the 95% prediction interval blew up to ~8e5 -- instead of the likelihood optimum ~0.18. Replace the clamp with a smooth reparameterization: the BFGS now carries an unconstrained theta in the lambda slot (no extra parameter), and the objective maps it to lambda = lo + (hi - lo) * logistic(theta) via boxcox.h::lambdaFromTheta, with lo/hi the actual lambda bounds. The gradient never vanishes, so the optimiser locates interior optima and only reaches a bound when the profile genuinely wants it. The warm-start lambda is seeded into the slot via thetaFromLambda; the reported lambda is lambdaFromTheta of the converged theta (no more out-of-bounds slot fixups). Effect: 267 -> lambda 0.18 (interval upper 8e5 -> 3e4); 162 still -> 2.0 (its profile is genuinely monotone for the flexible model, no collapse); 322 unchanged at 0.25. Shared engine; R and Python identical (267 logLik -284.7794, 162 -723.6260, 322 -450.9125 in both). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…erence) The optimiser's lambda (theta) derivative was a forward difference, (f(p+h) - llikValue)/h, reusing the cached objective `llikValue`. That value is NOT a reliable f(p) at the gradient point: the backward smoother and the structural-gradient loop in gradLlik run before the lambda slot and mutate persistent Kalman-filter state that llik()/llikAug() read, so a fresh llik(p) taken there differs from llikValue by ~1 (confirmed: llikValue 11.34 vs fresh 10.52 at the same p; objFun itself is idempotent at the top of gradLlik, and the line search leaves no staleness -- the shift is introduced inside gradLlik). The stale baseline made the lambda gradient ~1e6, which slammed theta to a logistic-saturated bound and froze lambda at the warm-start. Fix: compute the lambda gradient as a CENTRAL difference at the TOP of gradLlik, from the clean state the preceding objFun(p) leaves, before the smoother / structural loop perturb anything; store it and write it into the lambda slot in both the analytic and numerical gradient branches. Also clamp the warm-start seed out of the logistic's saturated tails so a bound-valued warm-start (testBoxCox) cannot start the optimiser where d lambda/d theta ~ 0. Effect: monthly M1[532] now finds PTS(2,G,D) AICc 281.7 -- near the global fixed-lambda best PTS(2,L,D) AICc 278.9 -- instead of the broken PTS(-2,D,N) AICc 525. Eurostat 267/162/322 unchanged (lambda 0.168 / 2.0 / 0.249). R and Python identical (267 logLik -284.7794, 162 -723.5872, 322 -450.9125). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Address the two issues flagged in CRAN review: - Document the return value of every exported S3 method via a @return block on the shared pts-methods roxygen header (regenerated Rd now carries a \value section). - simulate.pts and .pts_forecast_paths save and restore the user's .Random.seed when a seed is supplied, so the global environment is left unchanged (Armadillo's randn() routes to R's RNG under RcppArmadillo, so reproducibility goes through set.seed()). Update cran-comments.md and bump the DESCRIPTION date. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
config-i1
added a commit
that referenced
this pull request
Jul 23, 2026
Sparce matrices in the C++ code + fixes for the box-cox transform
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.