Add rolling tail dependence verb

Author: ttuff

docs/reference/verbs_stats.md

@@ -49,6 +49,26 @@ import xarray as xr
 corr = xr.corr(ndvi_z, climate_anom, dim="time")
 ```
 
-Rolling synchrony functions such as `cubedynamics.rolling_corr_vs_center` and `cubedynamics.rolling_tail_dep_vs_center` live outside the verbs namespace.
+Rolling synchrony functions such as `cubedynamics.rolling_corr_vs_center` remain outside the verbs namespace.
 
 Use these stats alongside transform verbs to build climate–vegetation synchrony analyses.
+
+### `v.rolling_tail_dep_vs_center(window, dim="time", min_periods=5, tail_quantile=0.8)`
+
+Compute a rolling contrast between variability in the upper tail and variability across the full window. The result preserves the input cube shape (e.g., `(time, y, x)`).
+
+```python
+from cubedynamics import pipe, verbs as v
+
+td = (
+    pipe(ndvi_z)
+    | v.rolling_tail_dep_vs_center(
+        window=90,
+        dim="time",
+        min_periods=5,
+        tail_quantile=0.8,
+    )
+)
+```
+
+This verb rolls along ``dim`` and, for each window, compares the variance of values above ``tail_quantile`` to the variance of the entire window. Higher values indicate stronger variability in the tails relative to the center.

src/cubedynamics/verbs/__init__.py

@@ -14,7 +14,7 @@
 from .custom import apply
 from .flatten import flatten_cube, flatten_space
 from .models import fit_model
-from .stats import anomaly, mean, variance, zscore
+from .stats import anomaly, mean, rolling_tail_dep_vs_center, variance, zscore
 
 
 def show_cube_lexcube(**kwargs):
@@ -61,6 +61,7 @@ def _op(obj):
     "month_filter",
     "flatten_space",
     "flatten_cube",
+    "rolling_tail_dep_vs_center",
     "variance",
     "correlation_cube",
     "to_netcdf",

src/cubedynamics/verbs/stats.py

@@ -111,4 +111,56 @@ def _op(obj: xr.Dataset | xr.DataArray) -> xr.Dataset | xr.DataArray:
     return _op
 
 
-__all__ = ["anomaly", "mean", "variance", "zscore"]
+def rolling_tail_dep_vs_center(
+    window: int,
+    *,
+    dim: str = "time",
+    min_periods: int = 5,
+    tail_quantile: float = 0.8,
+):
+    """Return a rolling "tail dependence vs center" contrast along ``dim``.
+
+    For each rolling window this computes the difference between variability in
+    the upper tail (values above ``tail_quantile``) and variability across the
+    full window. The verb preserves the original cube shape.
+
+    Parameters
+    ----------
+    window : int
+        Rolling window size in number of time steps.
+    dim : str, optional
+        Dimension to roll over (default: ``"time"``).
+    min_periods : int, optional
+        Minimum periods in window required to compute the statistic.
+    tail_quantile : float, optional
+        Quantile threshold defining the upper tail (default: ``0.8``).
+    """
+
+    def _op(obj: xr.Dataset | xr.DataArray) -> xr.Dataset | xr.DataArray:
+        _ensure_dim(obj, dim)
+
+        window_dim = f"{dim}_window"
+        rolled = obj.rolling({dim: window}, min_periods=min_periods)
+        constructed = rolled.construct(window_dim)
+
+        counts = constructed.count(dim=window_dim)
+        center_var = constructed.var(dim=window_dim, skipna=True, keep_attrs=True)
+
+        q = constructed.quantile(tail_quantile, dim=window_dim)
+        tail_vals = constructed.where(constructed >= q)
+        tail_counts = tail_vals.count(dim=window_dim)
+        tail_var = tail_vals.var(dim=window_dim, skipna=True, keep_attrs=True)
+
+        result = tail_var - center_var
+        valid = (counts >= min_periods) & (tail_counts > 0)
+        result = result.where(valid)
+
+        if isinstance(result, xr.DataArray) and obj.name:
+            result = result.rename(f"{obj.name}_tail_dep_vs_center")
+
+        return result
+
+    return _op
+
+
+__all__ = ["anomaly", "mean", "rolling_tail_dep_vs_center", "variance", "zscore"]

tests/test_rolling_tail_dep.py

@@ -0,0 +1,30 @@
+import numpy as np
+import xarray as xr
+
+from cubedynamics import pipe, verbs as v
+
+
+def test_rolling_tail_dep_vs_center_shape_and_nans():
+    np.random.seed(0)
+    time = np.arange(20)
+    y = np.arange(3)
+    x = np.arange(4)
+    data = xr.DataArray(
+        np.random.randn(len(time), len(y), len(x)),
+        dims=("time", "y", "x"),
+        coords={"time": time, "y": y, "x": x},
+        name="ndvi",
+    )
+
+    out = (pipe(data) | v.rolling_tail_dep_vs_center(window=5, dim="time", min_periods=5)).unwrap()
+
+    assert out.dims == data.dims
+    assert out.shape == data.shape
+
+    # First few windows should be missing because of min_periods
+    early = out.isel(time=slice(0, 4))
+    assert np.isnan(early).all()
+
+    # Later windows should contain finite values for most pixels
+    later = out.isel(time=slice(6, None))
+    assert np.isfinite(later.values).any()