[lag_transforms] Honor `skipna` in the stateful `update()` path (ExpandingMean/Std/Min/Max, ExponentiallyWeightedMean)

[simonez-tuidi]

Description

skipna (added in 0.0.17/0.0.18 — it is absent throughout 0.0.16) works correctly in transform()
for every transform in coreforecast.lag_transforms. It is silently ignored by update() for five
of them, so a NaN that arrives through the incremental path poisons the accumulator permanently even
when skipna=True was requested.

Verified on coreforecast 0.0.18 (Linux, cp310). skipna=True on every transform; history
[1, 2, 3, 4, 5]; then a NaN observation arrives via update(), then a valid 9.0:

transform          upd(+nan)  upd(+nan,+9)   expected
ExpandingMean            nan           nan     4.0000   <-- poisoned
ExpandingStd             nan           nan     2.8284   <-- poisoned
ExpandingMin             nan           nan     1.0000   <-- poisoned
ExpandingMax             nan           nan     9.0000   <-- poisoned
EWM(alpha=0.5)           nan           nan     6.5312   <-- poisoned
ExpandingQuantile     3.0000        3.5000     3.5000   ok
RollingMean(w=3)      4.5000        7.0000     7.0000   ok
RollingStd(w=3)       0.7071        2.8284     2.8284   ok

expected is transform() over the same full array read at the same position — i.e. the
transform/update pair disagree, and once the accumulator is NaN it never recovers.

For contrast, when the NaN is already in the history at transform() time, transform() and the
following update() are both correct for all of these — so the state built by transform is fine;
only the incremental step is broken:

history = [1., 2., nan, 4., 5., 6.], skipna=True
transform          update(hist)       expected   ok
ExpandingMean          3.600000       3.600000   OK
ExpandingStd           2.073644       2.073644   OK
ExpandingMin           1.000000       1.000000   OK
ExpandingMax           6.000000       6.000000   OK
EWM(0.5)               4.937500       4.937500   OK
RollingMean3           5.000000       5.000000   OK
RollingStd3            1.000000       1.000000   OK

Root cause

There is an exact invariant: update() honors skipna if and only if it delegates to _lib.

_lib.grouped_array._GroupedArrayFloat64 exposes _rolling_{mean,std,min,max,quantile}_update,
_seasonal_rolling_{mean,std,min,max,quantile}_update and _expanding_quantile_update — and every
transform whose update() calls one of those is correct, because skipna is forwarded into C++:

# coreforecast/lag_transforms.py:120-123  (_RollingBase)
def update(self, ga: "GroupedArray") -> np.ndarray:
    return getattr(ga, f"_rolling_{self.stat_name}_update")(
        self.lag - 1, self.window_size, self.min_samples, self.skipna
    )

There is no _expanding_mean_update, _expanding_std_update, _expanding_min_update,
_expanding_max_update or _exponentially_weighted_mean_update in _lib. Those five transforms
therefore carry their own NumPy accumulator update in Python, and none of them reads self.skipna:

# coreforecast/lag_transforms.py:391-394  (ExpandingMean.update)
self.stats_[:, 0] += 1.0                              # count incremented unconditionally
self.stats_[:, 1] += ga._index_from_end(self.lag - 1)  # NaN poisons the sum
return self.stats_[:, 1] / self.stats_[:, 0]

# coreforecast/lag_transforms.py:409-417  (ExpandingStd.update)
x = ga._index_from_end(self.lag - 1)
self.stats_[:, 0] += 1.0                     # Welford update, no NaN guard;
n = self.stats_[:, 0]                        # mean and M2 both become NaN
prev_avg = self.stats_[:, 1].copy()
self.stats_[:, 1] = prev_avg + (x - prev_avg) / n
self.stats_[:, 2] += (x - prev_avg) * (x - self.stats_[:, 1])

# coreforecast/lag_transforms.py:429-431  (_ExpandingComp.update -> ExpandingMin/ExpandingMax)
self.stats_ = self._comp_fn(self.stats_, ga._index_from_end(self.lag - 1))
# _comp_fn is np.minimum / np.maximum (lines 443, 455), which propagate NaN.
# np.fmin / np.fmax ignore it.

# coreforecast/lag_transforms.py:500-503  (ExponentiallyWeightedMean.update)
x = ga._index_from_end(self.lag - 1)
self.stats_ = self.alpha * x + (1 - self.alpha) * self.stats_

The EWM case also contradicts its own documented contract — ExponentiallyWeightedMean's docstring
at coreforecast/lag_transforms.py:487-488 promises:

skipna (bool): If True, exclude NaN values from calculations using forward-fill behavior.
When False (default), NaN values propagate through the calculation.

and exponentially_weighted_mean's docstring in coreforecast/exponentially_weighted.py:16-20 says
the same. transform() does forward-fill; update() does not.

The ExpandingMean/ExpandingStd/ExpandingMin/ExpandingMax docstrings
(lag_transforms.py:382-383, 402-403, 439-440, 451-452) all state "If True, exclude NaN values
from calculations" without scoping it to transform, so the current behaviour reads as a bug rather
than a documented limitation.

Suggested fix

Guard the five Python accumulators on self.skipna. No C++ change is strictly required:

• ExpandingMean.update / ExpandingStd.update — when skipna and x is NaN, skip the update
  entirely (do not increment the count, do not touch the running mean/M2) and return the current
  statistic.
• _ExpandingComp.update — use np.fmin/np.fmax as _comp_fn when skipna=True.
• ExponentiallyWeightedMean.update — when skipna and x is NaN, return self.stats_ unchanged
  (that is the forward-fill the docstring describes).

All four are per-group vectorized, so a np.where(np.isnan(x), ...) form keeps them branch-free.
Alternatively, add the missing _expanding_*_update / _exponentially_weighted_mean_update kernels
to _lib so all transforms follow the single delegating code path and the invariant becomes
structural rather than something each Python update() has to remember.

A regression test shaped like a transform/update consistency property would catch the whole
class: for each transform and each skipna, update() applied step by step over a series must equal
transform() over the full series read at the corresponding positions.

Use case

This is related to mlforecast's issue 704. We use mlforecast for retail/product-level demand forecasting. Products have a life-cycle: for part of the history a SKU does not exist, or is delisted / out of range. mlforecast requires a gap-free time grid, so those periods have to be materialized as rows, and today the only allowed filler is a number — in practice 0. That makes "did not sell" indistinguishable from "did not exist", and a 12-week rolling mean spanning 8 weeks of non-existence is computed over 12 observations instead of the 4 real ones.

The fix on our side is NaN: keep the row so the grid stays dense and features stay positional, but
exclude it from the statistics. That is exactly what skipna=True gives us, and it already works —
in transform(). We have filed the corresponding request against mlforecast to allow NaN targets
at fit/preprocess and to pass skipna through to these transforms.

But mlforecast calls transform() once at preprocess time and then update() once per step
during recursive prediction, and per new batch of observations in MLForecast.update(). So with the
current update() behaviour, skipna=True would give correct training features and then silently
diverge at prediction time for ExpandingMean, ExpandingStd, ExpandingMin, ExpandingMax and
ExponentiallyWeightedMean — the worst possible failure mode, since nothing raises and the training
matrix looks right. A NaN arriving through MLForecast.update() (a SKU going out of range mid-stream)
would additionally poison the state permanently, for every future horizon.

Closing this gap is what makes skipna usable end-to-end for the "sparse / intermittent life-cycle"
case rather than only for one-shot feature construction.

Repro

import numpy as np
from coreforecast.grouped_array import GroupedArray
import coreforecast.lag_transforms as ct

def mk(y):
    return GroupedArray(np.asarray(y, dtype=np.float64),
                        np.array([0, len(y)], dtype=np.int32))

hist = [1.0, 2.0, 3.0, 4.0, 5.0]

for make in [lambda: ct.ExpandingMean(1, skipna=True),
             lambda: ct.ExpandingStd(1, skipna=True),
             lambda: ct.ExpandingMin(1, skipna=True),
             lambda: ct.ExpandingMax(1, skipna=True),
             lambda: ct.ExponentiallyWeightedMean(1, 0.5, skipna=True),
             lambda: ct.ExpandingQuantile(1, 0.5, skipna=True),
             lambda: ct.RollingMean(1, 3, 1, skipna=True)]:
    tfm = make()
    tfm.transform(mk(hist))
    tfm.update(mk(hist))                        # incorporate the last real obs
    tfm.update(mk(hist + [np.nan]))             # a NaN observation arrives
    got = float(tfm.update(mk(hist + [np.nan, 9.0]))[0])
    exp = float(make().transform(mk(hist + [np.nan, 9.0, 0.0]))[-1])
    print(f"{type(tfm).__name__:28s} update={got!r:>10}  transform={exp:.4f}")

ExpandingMean                update=       nan  transform=4.0000
ExpandingStd                 update=       nan  transform=2.8284
ExpandingMin                 update=       nan  transform=1.0000
ExpandingMax                 update=       nan  transform=9.0000
ExponentiallyWeightedMean    update=       nan  transform=6.5312
ExpandingQuantile            update=       3.5  transform=3.5000
RollingMean                  update=       7.0  transform=7.0000

Versions: coreforecast 0.0.18, numpy 2.x, Python 3.10, Linux x86_64.

Two smaller, related inconsistencies

Separable from the above; happy to split into their own issues if preferred.

1. LocalBoxCoxScaler is the only _BaseLocalScaler subclass without a skipna parameter.
   LocalMinMaxScaler (coreforecast/scalers.py:187-188), LocalStandardScaler (:200-201) and
   LocalRobustScaler (:214-219) all accept it and forward it via
   _BaseLocalScaler.fit (:131: getattr(ga, f"_{self._scaler_type}_stats")(self.skipna)).
   LocalBoxCoxScaler.__init__ (:234-245) does not, so it inherits the class-level default
   skipna: bool = False (:121) and its overridden fit (:247-268) calls _boxcox_loglik /
   _boxcox_guerrero with no skipna argument — a NaN anywhere in a group makes that group's lambda
   NaN with no way to opt out. Related: the positivity guard at :255
   (if self.method == "loglik" and any(ga.data < 0)) does not fire for NaN, since nan < 0 is
   False, so NaN input reaches loglik silently.

   Confirmed working today for the scalers that do have it:

   x = [1., 2., nan, 4., 5.]
   LocalStandardScaler(skipna=True)   stats=[3.0, 1.58113883]   -> [-1.265, -0.632, nan, 0.632, 1.265]
   LocalStandardScaler(skipna=False)  stats=[nan, nan]          -> [nan, nan, nan, nan, nan]
   LocalMinMaxScaler(skipna=True)     stats=[1.0, 4.0]          -> [0.0, 0.25, nan, 0.75, 1.0]
   

2. scalers.Difference / differences.diff have no skipna and no documented NaN contract.
   Difference(1).fit_transform on [1, 2, nan, 4, 5] gives [nan, 1, nan, nan, 1], which is the
   mathematically right answer — a difference against a missing value is undefined. We are not asking
   for that to change. But Difference.tails_ (scalers.py:311, ga._tail(self.d)) can hold NaN,
   and Difference.update / inverse_transform (:314-340) then propagate it into the inverted
   series. A note in the docstring on what NaN tails mean for inverse_transform, or a skipna-style
   option that carries the last valid value in the tail, would remove the ambiguity for callers doing
   differencing on top of life-cycle-sparse series.

Note for downstream pinning

skipna does not exist in 0.0.16 (zero occurrences in rolling.py, expanding.py,
lag_transforms.py, scalers.py); 0.0.18 has it throughout. mlforecast currently declares
coreforecast>=0.0.15, so any mlforecast-side use of skipna needs that floor raised to whichever
release introduced it — worth stating explicitly in the release notes if it is not already.