arima: fix coef stderr - #1206

Conflicts in python/statsforecast/arima.py and tests/test_arima.py:

- upARIMA: kept main's cached, non-mutating `Z` dict, which supersedes this
  branch's `mod = {**mod}` copies (and drops the unused `T` update).
- ML/CSS optimizations: kept this branch's `partial`-based `obj_fn` +
  `approx_hess3` Hessian, layered onto main's `ml_obj`/`dist_tail_fit`
  bookkeeping (`args=ml_args` dropped since the objective is now bound).
- maInvert re-score: kept main's direct re-evaluation instead of re-running
  `minimize`, carrying `res.hess` through.
- Fell back to NaN when `cho_factor` finds the Hessian isn't positive
  definite, matching the existing singular-matrix handling; main's
  near-degenerate distribution tests hit this.
- tests: kept both sides' new tests.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_0119VbLrxzmn89AWrseL74fn

Three fixes on top of the numerical Hessian:

- Re-estimate the Hessian after maInvert re-parameterises the MA part, as
  R's optim(maxit = 0, hessian = TRUE) does. The Jacobian `A` was already
  evaluated at the re-parameterised coefficients while the Hessian came from
  the pre-inversion fit point. On AirPassengers ARIMA(2,1,1)(0,1,0)[12] this
  moves ma1's standard error from 0.03027 to 0.02919 (R: 0.02920); ar1/ar2
  and the ML fit now match R to 5 decimals too.

- Marginalise the distribution parameters out of var_coef. The t/skew-normal
  /ged branches sliced the arma block out of the Hessian and inverted that,
  giving variances conditional on sigma2 and the shape parameter. Inverting
  first and slicing after gives the marginal ones; the difference is under
  1% for t and ged but understates the intercept by ~99% for skew-normal,
  whose alpha is strongly correlated with the mean.

- Fold the three copies of the invert-or-NaN logic into `_coef_var`. The
  transformed branch used cho_factor, which additionally rejects Hessians
  that are merely indefinite and raises ValueError rather than LinAlgError
  on NaN input, so it needed a separate isnan pre-check and still crashed on
  the near-degenerate distribution fits.

test_coef_stderr now pins R's values (atol 2e-3 -> 1e-4) and covers ML as
well as CSS-ML; every one of these tests fails without the changes above.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_0119VbLrxzmn89AWrseL74fn

4.36 was recorded during the intermediate states of this branch, where
_coef_var used cho_factor: that also rejected merely-indefinite Hessians and
raised ValueError, which myarima swallows into ic = inf, so candidates were
dropped and 17 of the 174 Other series reselected.

Switching to np.linalg.inv dropped that. The search still rejects 4
candidates main keeps, but none of them win their series: the selected model
and AICc are identical to main for all 174, and the forecasts are
bit-identical. Both main and HEAD measure 4.4868, which the test's 1% rtol
does not admit against 4.36.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_019YXn1fursQmqpMazreQccZ

The parallel path runs every worker through GroupedArray's
_single_threaded_* wrappers, which cap BLAS at one thread. The serial
branches called GroupedArray directly, so they never got that cap.

It costs little until something makes LAPACK calls in the fitting loop. The
numerical Hessian adds ~1235 np.linalg.inv calls per M3 Other run on
matrices of size <= 8, and LAPACK's getri hands those to threaded BLAS3,
where waking the pool costs far more than the arithmetic.

Routing the six serial branches through the wrappers that already exist cuts
this branch's overhead on that run from +7.8% to +5.0% against main idle, and
from +22.4% to +12.9% with the other cores busy. It also stops the timings
swinging with machine load: run-to-run spread drops from ~1s to under 0.05s.
Forecasts are bit-identical.

Only numpy's OpenBLAS is capped -- the module-level ThreadpoolController
snapshots at import, before scipy's copy loads -- but capping scipy's too
measured no gain here, so that is left alone.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_019YXn1fursQmqpMazreQccZ

auto_arima_f raised "No suitable ARIMA model found" whenever the xreg was
rank-deficient -- two constant columns, which is what a static feature becomes
on a per-series model, and what generate_series feeds the dask/spark flow
tests. Every candidate carries the same xreg, so every candidate was rejected.

np.linalg.inv does not raise on a numerically singular matrix; it returns a
covariance with negative diagonal entries, whose square root is NaN, which is
exactly what checkarima rejects on. Main never saw this because scipy's BFGS
hess_inv is positive definite by construction, so the degeneracy stayed hidden
until the Hessian became a real one.

Decompose instead and read the spectrum. approx_hess3's central differences
resolve the Hessian to about sqrt(eps), so eigenvalues under that are
indistinguishable from zero: one below -tol is a genuine saddle and still
yields NaN, while the flat ones are directions collinear xreg leaves
unidentified and are dropped, inverting over the identified subspace. A
well-conditioned Hessian keeps every direction and the result is the plain
inverse, so the standard errors pinned against R are untouched.

This is a guard, not a cure -- a single constant regressor is equally
unidentified and stays undetected, since the design matrix is full rank and
only the likelihood is flat. Dropping collinear regressors up front is filed
separately.

eigh also raises LinAlgError on a Hessian holding NaN, which the previous code
only caught for inv.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_019YXn1fursQmqpMazreQccZ