Ets initialisation#405
Merged
Merged
Conversation
Replace the length(B)+nStatesBackcasting patchwork and the fractional dfForBack heuristic with a structural, estimate-gated df that is identical across initialisation methods and correct for non-coprime multi-seasonal models. Optimal df is preserved exactly for single/coprime models (verified across ETS/ARIMA/CES/GUM/SSARIMA/om); only the intended changes land. - helper.R: dfInitialsETSLevelSeasonal() -- identifiable ETS level+seasonal initial df via inclusion-exclusion over gcd of seasonal periods (the one new piece of maths; single seasonal -> m, coprime -> naive, non-coprime -> the shared-frequency drop). Plus dfGCD (Euclid). - adam: initials consume the same df however obtained. Under optimal they sit in B at the naive seasonal count, so only the ETS cross-block redundancy is subtracted; otherwise the identifiable count is added (ETS via the formula, ARIMA via initialArimaNumber, xreg initials only under "complete"). gradient/backcasting/complete unified. - CES/GUM/SSARIMA/SPARMA: backcasting initial df = the structural count (CES own formula; GUM/SSARIMA/SPARMA sum of component lags), no non-coprime correction (verified rank == count). Scale always counted. - om (Bernoulli): ETS initial df added under backcasting/complete/gradient, scale stays 0. omg inherits via omEstimator. - Scale is now an estimated parameter for every loss (concentrated likelihood), not only likelihood -- df_scale = 1 for continuous engines, 0 for occurrence. - Remove calculateBackcastingDF/dfDiscounter/dfDiscounterFit and the dfForBack plumbing (globals/utils-adam/adamGeneral/sma). Zero refs remain. Verified: nparam(backcasting) == nparam(optimal) == nparam(gradient) across all engines; non-coprime multi-seasonal (e.g. lags 4,8) corrected for all init methods incl. optimal (17->14); optimal single/coprime byte-identical. R suite 475 pass; the 8 auto.adam baseline pins now reflect the (intended) shifted model selection -- to be re-pinned. Python mirror + NEWS + man pending. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
SMA has fixed equal 1/order weights and no estimated smoothing, so D = F - g w' = 0: only the first fitted value depends on the initials and the order initial states collapse to a single identifiable level (rank 1). Correct df is 1 (level) + 1 (scale) = 2, independent of the order. The generic adam "use" path SMA is built on hardcoded df = 1 (missing the scale); set the SMA nParam table and logLik df explicitly in sma(), and use df = 2 in the order-selection IC (constant, so selection is unchanged). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Earlier commit wrongly set SMA df to a constant 2. SMA is AR(order) with fixed coefficients 1/order that sum to 1 -- a unit root -- so the model is non-stationary and all `order` initial states persist and are identifiable: rank(X) = order (verified by the probe for order 1/2/4/12). df = order (backcast initials) + 1 (scale) = order + 1. This grows with the order, so order selection now carries a proper complexity penalty (the old df = 1 had none). Cap the searched maxOrder to obs-3 so the AICc small-sample correction stays valid now that df can approach the sample size. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Init-method invariance (nparam backcasting==optimal==gradient), non-coprime multi-seasonal gcd correction, the dfInitialsETSLevelSeasonal helper against known ranks, estimate-gating (provided alpha not counted), regressors='adapt' delta parameters, scale for non-likelihood losses, om initials + no scale, CES/GUM/SSARIMA backcasting==optimal, and sma df=order+1. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Model selection shifts once backcasting initials count towards df: simpler models win (ETS persistence 10->3; ARIMA(0,2,2) instead of (2,1,2)) and two distribution flips (dgamma->dnorm, dnorm->dlnorm) as the IC now penalises complexity correctly. AICc values re-pinned to the new selections; the ARIMA test's seasonal arma[[2]] assertion dropped (non-seasonal model now selected). All intended consequences of the structural df fix. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
…a order+1 Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Mirror the R structural df fix. n_param.py: df_initials_ets_level_seasonal (inclusion-exclusion over gcd). Estimator: initials count the same however obtained (add for backcasting/complete/gradient, subtract the ETS seasonal redundancy for optimal), ARIMA via initial_arima_number, xreg under complete; scale counted for every loss. OM: _om_initials_df (Bernoulli, no scale); OMG: both sides' initials so OM-vs-OMG selection is fair; CES: structural initial count + scale. ADAM AICc matches R bit-for-bit (optimal and backcasting). Tests updated for the intended selection/df shifts (ETSX AIC pins, msarima AICc, gradient df-invariance test, AutoOM fixture given real occurrence dynamics). Python suite 738 pass, ruff/mypy clean. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Additive SSOE models have residuals affine in the full initial profile, so the one-shot least-squares gradient solve applies to ARIMA too, not only ETS. Generalise adam_gradientProbeBasis: for additive models it adds a slot-per-state basis over the ARIMA initial cells (each state's lag slots are free; the rank-revealing pivoted QR drops the redundant directions), alongside the existing ETS grouping. Thread componentsNumberARIMA and lagsModelAll through the dispatcher. Verified: at fixed dynamics, gradient reaches the optimal-initials SSE exactly (ARIMA(2,1,2): 604.5721 vs optimal 604.5723) and beats backcasting (618.12), for pure ARIMA, seasonal SARIMA and mixed ETS+ARIMA; df matches optimal. Multiplicative ARIMA (whose Gauss-Newton Jacobian is ETS-only) and xreg fall back to backcasting. No C++ change. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
These are additive SSOE models, so the affine least-squares gradient solve applies (residuals affine in the initial profile). Route their fit sites (cost function + final fit) through adam_fitOrGradient with the additive slot-per-state probe basis (etsModel=FALSE, arimaModel=TRUE, the engine's Etype), instead of calling adamCpp$fit directly. No C++ change. Verified: at fixed dynamics gradient matches the optimal-initials SSE exactly (SPARMA 368.557, CES-none 367.303) and differs from backcasting (genuinely solving). For GUM and seasonal CES where the joint optimisation of the transition/measurement elements plus initials diverges, gradient is markedly more robust (GUM SSE 9732 -> 368; CES full 543 -> 475). Multiplicative CES/GUM and xreg fall back to backcasting. R suite 516 pass; 40 gradient tests incl. a CES/GUM/SSARIMA/SPARMA case. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Generalise adam_gradient_probe_basis (error_type / components_number_arima / lags_model_all args + ARIMA per-slot probes with the additive check) and the adam_fit_or_gradient dispatcher so the affine least-squares solve covers every additive single-source-of-error engine, mirroring the R side. Route both CES fit sites (ces/cost_function.py, ces_model.py) through the dispatcher. At fixed dynamics the gradient fit matches the optimal-initials SSE and differs from backcasting for MSARIMA and CES; multiplicative-error ARIMA and xreg/custom-loss still fall back to backcasting. Update test_gradient.py: MSARIMA solves additive ARIMA initials, CES solves additive SSOE initials (replacing the stale fall-back assertion). NEWS on both sides documents the additive-SSOE extension and reconciles the launch bullets. Python suite 738 passed; ruff/mypy clean. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Occurrence models now solve their initials by gradient beyond pure ETS: - om() with ARIMA and/or xreg components. The analytic Jacobian companions (adamGradient.h) are ETS-only, so gradientSolve now drops the occurrence Gauss-Newton to its finite-difference Jacobian when nArima/nXreg/constant are present — the FD path drives the same loop through the fully general adamW/F/Gvalue residual pass. The probe basis gains xreg state cells and (on the occurrence path) allows ARIMA/xreg; the demand path is unchanged (xreg and multiplicative ARIMA still fall back there). - omg(), the coupled two-model occurrence. New adamCore::gradientSolveGeneral solves both sub-models' initial profiles jointly over the shared probability residual r = o - aFit/(aFit+bFit) by Gauss-Newton with an FD Jacobian (line search + Levenberg-Marquardt), covering ETS/ARIMA/xreg on either side. Bound for R (Rcpp) and Python (pybind). omgCF_local runs the coupled solve in the objective; omgFinalFit routes each reported side through the single-model dispatcher so its standalone initials are gradient-solved too. Shared C++, so R and Python stay in parity. Verified at fixed dynamics the gradient fit reaches optimal-quality likelihood and differs from backcasting (om ARIMA/xreg; omg ETS/ARIMA/xreg). Non-gradient initials are byte-unchanged. R gradient tests 40 -> 50; test_oes green. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Mirror the R-side occurrence gradient extension (shared C++, so parity is automatic): - Probe basis (utils/gradient.py) gains xreg state cells and o_type-aware gating: the occurrence path allows ARIMA/xreg (finite-difference Jacobian); the demand path is unchanged. adam_fit_or_gradient threads xreg_number. - om.py / om_cost.py pass xreg_number through, so OM with ARIMA/xreg profiles its initials by gradient instead of falling back. - omg_cost.py: _omg_gradient_profiles runs the coupled adamCore.gradientSolveGeneral over both sub-models' shared probability residual, then one forward pass from the solved profiles; omg_cf routes gradient through it (backcasting fallback out of scope / on the direct-link seed saturation). Verified: OM ARIMA logLik -75.85 (beats backcasting), OMG ETS -69.90 (near optimal), OMG ARIMA -63.63 (vs backcasting -71.92) — mirrors R. Three new tests; Python suite 741 passed; ruff/mypy clean. NEWS on both sides. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
omg reported fitted from the standalone per-side refits but logLik from the coupled cost — so the reported fitted and logLik disagreed (Bernoulli of the fitted != logLik), in both languages identically. Fix, applied identically to both so parity is preserved: the top-level fitted probability now comes from the same coupled recursion the optimiser minimised (re-run the joint cost at the final B via a new returnFitted / return_fitted flag), and logLik is the Bernoulli of that coupled probability. So sum(o*log p + (1-o)*log(1-p)) == logLik exactly, for every loss. The standalone modelA/modelB sub-objects are unchanged (kept for diagnostics). - R: omgCF_local gains returnFitted; the final object's fitted/logLik derive from it. returnFitted is carried in nloptrArgs (so nloptr's .checkfunargs passes) and overridden via modifyList for the post-hoc coupled fitted. - Python: omg_cf gains return_fitted; OMG.fit computes the coupled fitted and its Bernoulli logLik. The now-dead _joint_fitted fallback is removed. For loss="likelihood" logLik/IC are unchanged; fitted/residuals change numerically for every omg fit; fit-only losses now report a proper Bernoulli logLik (as om() does). Verified R and Python bit-identical on shared data (backcasting + gradient, likelihood + MSE); 32 live R-parity omg tests pass; R gradient tests +4 consistency assertions (54); Python 741 + test_omg green. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
The occurrence Gauss-Newton solve used a finite-difference Jacobian for
ARIMA/xreg models (nFree forward passes per iteration) because the analytic
companions are ETS-only. But for additive-error occurrence (E="A", the
odds-ratio/inverse-odds default) the whole state->yhat map is linear over the
full component vector, so the exact Jacobian needs one pass, not nFree.
- gradientPassJacobianOccurrenceLinear: single-model additive-error occurrence,
raw-matrix propagation (d(yhat)=w'S, update F*S + g*dError) covering
ARIMA/xreg/constant. Bit-identical to the companion pass for pure additive
ETS (there jacGv=jacGy=0, jacGe=g).
- combinedJacobian in gradientSolveGeneral: coupled omg analytic Jacobian when
both sides are additive-error -- propagates the two sides' sensitivities
jointly through the shared probability p=aFit/(aFit+bFit), with closed-form
dp/dfitA, dp/dfitB and the occurrenceError('g') split derivatives.
Shared C++, so R and Python stay bit-identical. Multiplicative-error occurrence
with extra components keeps FD; pure-ETS occurrence unchanged.
Controlled before/after (same data): om-ARIMA 550->353ms (1.6x), omg-ARIMA
4966->2761ms (1.8x), omg two-sided additive-ETS 1076->129ms (8.3x). The cleaner
Jacobian reaches an equal-or-slightly-better optimum, so these fits' logLik
shift marginally (e.g. omg-ARIMA -44.48 -> -43.93).
Tests: R gradient 54, oes 34; Python gradient/omg 93, full suite 741; live
R-parity omg 32 and om 114 all pass. The sharper convergence widens the flat
ARIMA-ridge coefficient spread between the R(-O0 dev) and Python(-O2) builds
from ~1e-2 to ~5e-2 while every observable still matches tightly, so test_coef
uses a coefficient-only identifiability-floor tolerance for ARIMA/xreg
scenarios (observable checks stay tight). NEWS on both sides.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Update the branch
Dated copy of the exponential-smoothing benchmark (M1 + M3 + Tourism, 5315 series, automatic ETS selection) re-run against the current package versions, with aeon AutoETS added alongside statsforecast / sktime / skforecast and smooth's own ADAM / ES. - New aeon branch: aeon.forecasting.stats.AutoETS with iterative_forecast. - Forced a fork multiprocessing context so the ProcessPool workers inherit the notebook namespace (Python 3.14 no longer defaults to fork). - Records the full environment (python 3.14.4, numpy 2.3.5, pandas 2.3.3, smooth 1.0.7, statsforecast 2.1.1, sktime 1.0.1, skforecast 0.23.0, aeon 1.5.0) so the run is reproducible; numpy/pandas were pulled back by the statsforecast update and now appear in the version cell. Results (mean RMSSE / median time, 0 failures across 14 methods, 17.3 min): statsforecast AutoETS best accuracy (1.924); smooth ES/ADAM optimal ~1.947 close behind; aeon AutoETS the speed leader (0.023s) at 2.035. sktime AutoETS blows up on some series (mean RMSSE ~1e5) -- a known upstream statsmodels ETSModel instability, not a config issue. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
The "fitted equals pA/(pA+pB)" test asserted the old standalone per-side combination, which the fitted/logLik consistency fix (commit 5bf9f40) intentionally replaced: omg's fitted now comes from the coupled recursion the optimiser minimised. The Python mirror of this test was already updated; this is the missing R side. Reworked to the new contract -- fitted in (0,1) and its Bernoulli equals the reported logLik exactly. test_omg.R: 110 pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
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.
What issue does this PR fix?
Fixes in the df calculation for backcasting.
Language
AI assistance
If yes — how was AI used?
Claude. Always Claude!