Method overview¶

This page explains what RSD does, why it exists, and when to prefer it over classical baselines. For the full mathematical treatment, simulation study, and global validation, see the companion paper (Tsilimigkras 2026, submitted to Water Resources Research).

The problem RSD solves¶

Classical standardized drought indices (SPI, SSI, SDI) standardize each series against its own distribution. That works for monitoring a single station against its own history, but it fails for the question that motivates most climate-attribution and climate-projection work:

How does this scenario / model / period compare to that one?

When each series is mapped to its own reference distribution, the real magnitude differences between scenarios disappear. A 25 % drier future and the historical baseline both look like climatology under self-standardization. The thing you actually want to measure is gone.

RSD fixes this by anchoring all comparisons to a single, fixed reference distribution: target values are evaluated against the reference's CDF rather than against their own.

Fixing the reference creates three interlocking problems, each of which requires its own design step. The three steps are not optional - skip one and the others stop working.

The three-problem chain¶

Problem 1: A fixed reference breaks the ECDF¶

A purely empirical CDF (ECDF) only spans the values present in the reference sample. Target values that fall outside the reference range get mapped to rank 1 or rank n, which collapses them all to the same extreme z-score regardless of how far out they actually are. The information you wanted to preserve (extreme magnitude) is lost.

Solution. Splice GPD tails onto the ECDF core. Fit a Generalized Pareto Distribution to each tail of the reference sample so that out-of-range target values get a smoothly extrapolated CDF probability instead of being capped.

Problem 2: GPD needs enough exceedances to fit¶

Even with a long reference (50 years monthly = 600 values), a per-month fit has only ~50 values; a 10 % tail threshold leaves you ~5 exceedances - far too few for a stable GPD. The fit is dominated by noise.

Solution. Pool all twelve months into one sample so the GPD sees ~600 reference values total. A 10 % tail then provides ~60 exceedances, which is enough for a stable shape and scale estimate.

Problem 3: Pooling needs deseasonalization¶

You cannot pool raw monthly values across the seasonal cycle - the pool would just describe the seasonal envelope, not the residual variability the tails are supposed to model.

Solution. Deseasonalize first: subtract a per-month location and divide by a per-month scale (median and IQR by default; mean and std as an option). Pool the residuals, then fit the GPD on those.

How the three solutions interlock¶

You cannot pool without deseasonalizing first. You cannot fit tails without pooling first. You need tails because you chose a fixed reference. Removing any one step breaks the others.

The full pipeline¶

flowchart LR
    A[Raw monthly series] --> B[Accumulate<br/>scale-month windows]
    B --> C[Deseasonalize<br/>per calendar month]
    C --> D[Pool residuals<br/>across all months]
    D --> E[ECDF core<br/>+ GPD tails]
    E --> F[Standard normal<br/>z-score]

For bounded variables an optional logit pre-transform is applied at the very front of the pipeline, between the raw series and the accumulation step; everything downstream then operates on real-valued data (see Bounded variables below).

The three methods, side by side¶

rsd bundles two classical baselines next to the hybrid method so you can compare them on the same input with one function call. All three share the same standardize(...) / standardize_xr(...) signature.

`method=`	Pipeline	When to use
`"rsd"`	Deseasonalize -> pool -> ECDF core + GPD tails -> invnorm	Cross-series comparison under a fixed reference (the package's reason for being).
`"monthwise_ecdf"`	Per-month ECDF -> invnorm	Single-series monitoring; matches the classical SPI ECDF baseline.
`"monthwise_parametric"`	Per-month parametric fit (gamma by default, `floc=0`) -> invnorm	Single-series monitoring with the canonical Stagge SPI convention. Useful as a benchmark.

The two baselines are not new methodology; they are bundled so the paper and your downstream work can compare RSD against the canonical alternatives without re-implementing them.

Bounded variables¶

Some hydroclimate drought variables live on a bounded interval [L, U] where both bounds are physical. Soil moisture is the canonical case (bounded between 0 and the saturation moisture content θ_sat); relative humidity is another (bounded between 0 and 1). Most other drought variables (precipitation, discharge, P-PET, TWS, ET, SWE) are non-negative but unbounded above and do not need this pathway.

When the data approaches either bound, the residual distribution becomes left- or right-skewed in a way that ordinary location-and-scale deseasonalization cannot remove (IQR collapses against the wall; the distribution shape becomes month-dependent). Exchangeability may break for the compressed months.

The optional logit pre-transform (bounds=(L, U)) maps the bounded interval into unbounded space before deseasonalization:

y = log( (x - L) / (U - x) )

The transform is applied per timestep, before accumulation, so the rest of the pipeline (accumulation, deseasonalization, pooling, tail fitting) operates on real-valued data with no boundary effects. Choose between agg="sum" and agg="mean" based on what you want to standardize (an accumulated water-content quantity vs an average state); both work cleanly in logit space.

Each bound may be a scalar (applied uniformly) or, for standardize_xr, a DataArray broadcastable against the input grid (per-pixel saturation values). bounds=(L, U) is only valid with method="rsd".

Verifying the pooling assumption¶

Pooling across months assumes that, after deseasonalization, the residual distributions are exchangeable across calendar months (same shape, possibly different finite-sample noise). On well-behaved data this holds; on data with month-dependent compression or other quirks it may not.

rsd.diagnose(values, months, name, ...) runs the full exchangeability diagnostic in one call: it prints a configuration overview, renders before/after KDE panels, runs an Anderson-Darling omnibus test and per-month Kolmogorov-Smirnov leave-one-out tests, and reports the results. Pass bounds=(L, U) (and optionally auto_bounds=True) to compare the baseline pathway against the logit-bounded pathway side by side and see whether bounded handling restores exchangeability.

See the Diagnostics showcase for a worked example.

What RSD does not claim¶

RSD does not produce a "true probability under the target distribution" when applied out of sample. It produces a standardized anomaly relative to a fixed reference distribution; that distinction matters and is stated explicitly in the package documentation and manuscript.
RSD does not share the per-month statistical frame of classical monthwise methods. The classical baselines (and the bundled monthwise_ecdf / monthwise_parametric) fit one model per calendar month, so a z = -2 in January is "two SD below January's distribution". RSD pools across months after removing the seasonal location and scale, so its z = -2 is "two SD below the pooled deseasonalized residual distribution". Both are valid standardized indices, but cross-comparison with monthwise SPI values requires recognising the change in reference frame.
RSD is not a generic time-series standardization framework. The scope is hydroclimate drought indices under reference-based evaluation; other applications may need different design choices.