OPG vcov calculation#406
Merged
Merged
Conversation
The observed Fisher Information (numerical Hessian of the log-likelihood) is indefinite when a smoothing parameter is estimated at an active bound (e.g. a trend smoothing beta pinned at 0), so its inverse is not a valid covariance and summary() warns "Observed Fisher Information is not positive semi-definite". New covarOPG(): the outer-product-of-gradients estimator J = sum_t s_t s_t' of the per-observation scores, PSD by construction, so it returns finite standard errors for boundary-but-identified parameters. Scores come from central- differencing pointLik() (the per-observation log-density, already universal over distributions and model types) at parameters perturbed through the standard adam() refit, preserving the selected ETS structure and lags. Exposed as vcov(object, opg=TRUE); the default (opg=FALSE) is unchanged. Verified: at the AirPassengers MAM boundary (beta=0) OPG is PSD and gives a finite beta SE where the Hessian is indefinite; in an interior case OPG matches the Hessian (ANN alpha SE 0.041 vs 0.050) and both track the bootstrap. v1 scope is pure ETS (persistence, phi, level/trend/seasonal initials, backcast or optimal). ARIMA / xreg / constant fall back to the Hessian with a message (their initial-perturbation map is a follow-up). loss must be "likelihood". New tests: test_vcov_opg.R (9). test_adam/test_es unaffected. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Providing a fitted model's own ARIMA initials back to adam() produced a different (wrong) likelihood for identical reported parameters, e.g. ARIMA(1,1,1) on BJsales: -265.88 estimated vs -779.33 when the same arma and initial$arima were provided. Root cause: an asymmetry between how ARIMA initials are placed during estimation vs when provided. The estimated path (adam_filler) spreads the free initial parameters across the ARI state rows as ariPolynomial %*% initials, and the collector reports initial$arima by reading the last ARI row and dividing by tail(ariPolynomial). But the provided path (adam_creator) dumped the reported values raw into the last row only, dropping the ariPolynomial factor -> a different initial state -> a different likelihood. adam_creator now applies the same ariPolynomial spread as the filler, so the collected initials round-trip exactly. The ARIMA "initials" are the minimal free parameters (not the raw, rank- deficient state block), so transforming them on placement is the correct fix rather than storing raw states. Verified round-trip for ARIMA(1,1,1) and (2,1,2); test_adam / test_ssarima / test_gradient unaffected. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Now that provided ARIMA initials round-trip (previous commit), the OPG covariance is no longer restricted to pure ETS: - Initial states are perturbed as free parameters ONLY for initial="optimal". For backcasting / gradient / complete the initials are derived, not estimated, so the score is taken over the dynamics parameters (persistence / phi / arma) only and the initials are re-derived (via the model's own initialType) at each perturbed evaluation. This matches how the Hessian FI treats them and avoids the meaningless "provide the derived initials forward" path. - Extended the parameter -> slot map to arma (ar/ma), ARIMA initial states (ARIMAState*), and xreg coefficients, alongside the ETS persistence/level/ trend/seasonal. - Added a reproduction guard: the model rebuilt at the unperturbed parameters must reproduce logLik(object); otherwise (e.g. regressors="select"/"adapt", which re-selects on refit) it falls back to the Hessian. This replaces the hard pure-ETS gate with a correctness check. Verified PSD with finite SEs for MAM (backcasting 3x3, optimal 16x16 incl initials), ARIMA (optimal 4x4, backcasting 2x2), ETSX (optimal, incl xreg); regressors="select" falls back. Moved covarOPG from adam.R to helper.R. test_vcov_opg.R expanded to 16; test_adam/ssarima/gradient unaffected. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
covarOPG now picks the refit engine from the fitted object instead of always using adam(): the state-space ARIMA engines have their own ARIMA representation, so an adam() refit reconstructs a different model. ssarima is now wired (its data argument is `y`, and it reuses the arma / ARIMAState perturbation map), so vcov(opg=TRUE) computes a PSD OPG covariance for ssarima models (optimal: incl ARIMA initials; backcasting: dynamics only) instead of falling back to the Hessian. CES and GUM (complex-smoothing / free-transition parameterisations) and sparma (distinct coef layout and its own reproduction quirk) are not yet mapped and return NULL so the caller falls back to the Hessian -- safe, and tracked as a follow-up. The reproduction guard still protects every path. test_vcov_opg.R -> 22 (ssarima optimal + backcasting, CES/GUM fallback); adam / ssarima unaffected. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Split the OPG covariance into a reusable core and an engine-specific closure, so the bespoke CES / GUM / sparma variants can share the assembly: - covarOPGCore(object, parametersNames, perturbedPointLik, stepSize) does the reproduction guard, the central-difference score loop, J = S'S and the inversion with the near-zero-information drop. perturbedPointLik(nameJ, delta) returns the per-observation log-likelihood for the model rebuilt with that parameter perturbed (nameJ=NULL -> base fit), or NULL on failure. - covarOPG builds that closure for the ETS/ARIMA family (adam + ssarima) and calls the core. Behaviour unchanged (adam + ssarima compute OPG; CES/GUM/sparma still fall back to the Hessian via the reproduction guard). test_vcov_opg.R 22 pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
CES has its own parameterisation (complex smoothing alpha_0/alpha_1 = Re/Im of a, beta_0/beta_1 for seasonal b; the remaining coefficients are initial states in profileInitial), so covarOPG's adam refit cannot rebuild it. New covarOPGces perturbs a perturbed clone and re-fits via ces(model=clone) -- which reads parameters$a, parameters$b and profileInitial and holds them -- so both the smoothing parameters and the (optimal) initials are perturbed exactly. This is the same mechanism the CES Hessian FI uses (model=object), and it sidesteps the ces(initial=<numeric>) path, which re-optimises rather than holding the initial. vcov(opg=TRUE) dispatches to covarOPGces for CES objects. Works for the non-seasonal model (optimal: complex smoothing + level/potential initials; backcasting: smoothing only). Seasonal CES stores its initials across the seasonal lags with fewer free coefficients than stored cells, so the flat mapping is guarded off and falls back to the Hessian -- a follow-up. GUM/sparma still fall back. test_vcov_opg.R -> 26. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Seasonal CES stores its initial states in a (states x head-lags) profileInitial with fewer free coefficients than stored cells and unnamed seasonal coefficients, so the flat name-based map fell back. Now mapped by position: - covarOPGCore is index-based (perturbedPointLik(j, delta)), so it handles the duplicate/empty coefficient names seasonal CES produces; it takes the named coefficient vector for step sizes and result dimnames. - covarOPGces maps the initial coefficients by position: the non-seasonal states (level, potential) are a single free value repeated across the head columns (perturb the whole row), and the seasonal states fill their rows across the head columns in column-major order (verified against coef()). adam/ssarima refactored onto the index-based core (behaviour unchanged). Seasonal CES now computes a PSD OPG covariance (finite SEs for the identified parameters, Inf for the weakly-identified initials). test_vcov_opg.R 28 pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
GUM's coefficients split across three slots that gum(model=<clone>) reads: the persistence g1..gK (persistence vector), the free transition matrix F<row><col> (transition, row-major), and the initial states vt1..vtK (held via the coefficient vector B). New covarOPGgum perturbs the matching slot on a clone and re-fits -- the same mechanism the GUM Hessian FI uses. vcov(opg=TRUE) dispatches to it via gumChecker. Works for optimal initialisation (non-seasonal and seasonal, incl the vt initials). Backcasting falls back: the backcasted state is not forward- reproducible from a clone, and gum()'s fully-provided backcasting refit hits an empty-B nlopt error (a separate GUM bug), so the reproduction guard correctly routes it to the Hessian. sparma still falls back. test_vcov_opg.R 31 pass; test_gum unaffected. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
sparma does not reproduce exactly even via model=object (its own reconstruction bug), so the OPG reproduction guard routes it to the Hessian. Add a test pinning that behaviour until sparma's reconstruction is fixed. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Two issues surfaced by the GUM OPG backcasting path: 1. gum() errored "x0 must have length > 0" when persistence, transition and the initials were all provided (empty parameter vector): nloptr was called with a zero-length x0. Added an empty-B guard that evaluates the cost once instead (mirrors the om() guard). This is a general fix -- providing a fully-specified GUM now works, not just for OPG. 2. covarOPGgum now handles backcasting like the other engines: initials are perturbed as free parameters only for initial="optimal" (clone perturbation); for backcasting/gradient/complete the score spans g / F only and the initials are re-derived by re-fitting with the perturbed persistence / transition and the model's own initialType. This needs (1) to refit. 3. covarOPGCore gained a Moore-Penrose pseudo-inverse fallback (symmetric eigen-decomposition, dropping near-zero-eigenvalue directions) for when the score matrix is ill-conditioned/collinear -- as GUM backcasting's is (the F22 direction dominates). Still PSD; benefits every engine. GUM OPG now computes for optimal, seasonal AND backcasting. test_vcov_opg.R 36; test_gum unaffected. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
sparma() inverted the arma-provided flag (providing coefficients set arEstimate=TRUE, re-estimating them instead of holding); the creator already fills matF/vecG from armaParameters when !arEstimate, so the fix is to hold on a non-NULL side and estimate only a NULL one. A numeric initial is now front-padded to the dense companion length and placed via a new initialType=="provided" branch in the filler, so a fitted model's free ARIMA state initials round-trip exactly. An empty-B case (all held) evaluates the cost once (res<-NULL) instead of erroring. With held arma + provided initials reproducing the fit, covarOPGsparma computes a genuine PSD OPG covariance (initials perturbed under optimal; arma-only under backcasting/gradient), wired into vcov.adam. Follow-up: the Hessian vcov path for sparse-order sparma still errors in the adam initial collector (tail(ariPolynomial) multi-valued at utils-adam.R:1305/1707) -- orthogonal to OPG, which is becoming default. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
The omg re-entry (vcov's re-fit path) mis-reconstructed the fit: it re-supplied each sub-model's persistence/phi/arma as *fixed* values, which flips their Estimate flags off in the checker so adam_filler stops consuming those slots from the injected joint B and shifts the whole vector -- a boundary alpha=0 pushed the level initial to 0, so the re-fit logLik drifted from the original (e.g. -83.91 vs -83.03). The joint-B split also used the initialiser's reduced nParamsA (which drops provided params) instead of the fitted-object A-count, leaking A's params into B. Now the re-entry keeps every coefficient estimated and fed back 1:1 from the joint B, and the split uses nParamsA_use, so omg(model=object) reproduces the original fit exactly (optimal and backcasting). This also corrects the observed-FI vcov.omg, which relied on the same re-fit. With reproduction restored, covarOPGom / covarOPGomg compute the OPG/BHHH covariance for the occurrence models (Bernoulli / coupled Bernoulli score), wired into vcov.om / vcov.omg via opg=TRUE. A shared perturbAdamClone helper mutates a fitted clone's structured slots for the single-model om; omg perturbs the joint B directly. covarOPGCore now falls back to a one-sided difference when a central step leaves the feasible region (a negative persistence gives NaN under the occurrence link at a boundary alpha=0), keeping the estimator usable at bounds. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Ports R's pointLik.om to the Python occurrence models: point_lik(log=True) returns log(p_t) when the event occurred and log(1-p_t) otherwise, over the fitted (coupled, for OMG) probability. sum(point_lik()) equals loglik. Verified bit-for-bit against R on shared binary data (max abs diff ~2e-16 for both OM and OMG). This is the per-observation score the OPG/BHHH covariance path consumes on the R side. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Introduces type=c("opg","hessian","bootstrap") across vcov / confint /
summary / reapply / reforecast for the adam / om / omg methods, resolved by
a shared covarTypeResolver(). Default is now "opg" (the OPG/BHHH covariance
wired in this branch), so vcov() and the standard errors in summary() are
PSD-by-construction and finite at boundary estimates. type="hessian" keeps
the observed Fisher Information; type="bootstrap" the bootstrap. The old
bootstrap=TRUE switch still works but warns (deprecated).
reapply()/reforecast() hold any parameter with a non-finite OPG variance
(one the data does not identify -- an initial that washes out at a boundary
smoothing parameter) at its estimate, so the coefficient resampling stays
feasible instead of erroring in mvrnorm.
Docs: @param type + OPG/BHHH reference (Berndt-Hall-Hall-Hausman 1974) on the
reapply/reforecast page; NEWS updated. Tests migrated to the type= API.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
vcov(type="hessian") on a sparse-order sparma (e.g. ar=2,ma=1 with phi1 absent) errored: the observed-FI path re-fitted via adam(model=object), whose dense ARIMA-polynomial reconstruction yields an NA for the sparse structure (tail(ariPolynomial) -> 'missing value where TRUE/FALSE needed', utils-adam.R:1305/1707). sparma now computes its FI natively via a central- difference numerical Hessian of its own log-likelihood, reusing the same coefficient->refit map (sparmaModelAtParams) as the OPG path, so type="hessian" returns a finite covariance for every order combination. type="opg" (the default) was already fine. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Ports R's pointLik.adam: point_lik(log=True) returns the per-observation log-likelihood of the fitted model via calculate_likelihood (distribution/ scale aware) plus, for occurrence models, the zero-observation differential entropy and the occurrence Bernoulli term. sum(point_lik()) equals loglik. Verified bit-for-bit against R on shared data (max abs diff ~1e-13 across dnorm/dgamma/dlaplace/dlnorm). Completes the point_lik family (ADAM/OM/OMG), the per-observation score foundation for the OPG covariance. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Ports the OPG covariance and the type=('opg','hessian','bootstrap')
selector to Python for ADAM, OM and OMG, mirroring the R side; default is
now 'opg'. New covar_opg core (var_covar.py) builds J = sum_t s_t s_t'
from per-observation scores with a reproduction guard, one-sided-at-bound
differences, pseudo-inverse fallback and infinite variance for
unidentified directions -- a direct translation of R's covarOPGCore. The
score at a perturbed B is obtained via a new return_fitted flag on CF
(ADAM) and om_cf (OM); omg_cf already had one. Each engine forms its
density: ADAM the concentrated distribution likelihood (scaler +
calculate_likelihood), OM/OMG the (coupled) Bernoulli. Non-likelihood
losses and intermittent ADAM fall back to the Hessian, as in R.
type= is threaded through vcov/confint/summary (resolve_covar_type
handles the deprecated bootstrap=True -> DeprecationWarning). Verified
against R vcov(type='opg') on shared data: OM ~7e-17, ADAM/OMG ~1e-5
(finite-difference noise); PSD and genuinely distinct from the Hessian.
The r_parity summary/confint comparisons are pinned to type='hessian' on
both sides so their Hessian-tuned tolerances stay valid.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
reapply()/reforecast() gain the type= selector and, under the new opg default, hold any parameter with a non-finite OPG variance (one the data does not identify) at its estimate before the multivariate-normal draw -- otherwise rng.multivariate_normal chokes on the infinite entry. Mirrors the R-side reapply.adam fix. Full Python suite green (755 passed). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
Clarify (R covarOPGCore + Python covar_opg) that dropping a parameter to infinite variance is intentional, not a numerical artifact: at a boundary MLE such a parameter is weakly identified (its OPG score is tiny at the optimum even though it still enters the likelihood), so its OPG variance is genuinely enormous. Inf is a clearer 'effectively unidentified' signal than a huge finite number and makes reapply()/reforecast() hold the parameter at its estimate instead of injecting nonsense-wide draws. Behaviour unchanged; type='hessian' gives a finite (still large) SE there. Investigated a scale-invariant (correlation-normalised) rank check that turns the Inf into the finite ~1/diagJ variance; reverted because it made reapply sample the weakly-identified parameter wildly with no gain in meaning. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRdQhv3c5fgxm7GaequJtv
# Conflicts: # NEWS # python/NEWS.md
Python: reapply()/reforecast() crashed with "assignment destination is read-only" when the OPG covariance had a non-finite diagonal entry (a weakly-identified parameter, e.g. AAdN optimal initials at a smoothing bound). np.asarray(vcov_df) returns a read-only view into the DataFrame block, and the new non-finite row/column zeroing wrote into it in place. Take a writable copy (np.array) before zeroing. R: the reapply generic gained type/heuristics args, but reapply.adamCombined kept the old (object, nsim, bootstrap, ...) signature, tripping the S3 generic/method consistency WARNING (which the r-lib check treats as an error). Align the method signature and thread type/heuristics/nsim into the inner reapply() call (nsim was hardcoded to 1000). Import utils::combn/modifyList to clear the undefined-globals NOTE from the new helper.R/omg.R code. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014uFNLn5iC7hvWyBG8LtbcF
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?
Implements the new way to calculate vcov().
Language
AI assistance
If yes — how was AI used?
Claude Code