Homogeneity Tests

rucola.SNHTTest(alpha=0.05, min_years_from_end=5, min_relative_signal=1.0)

Bases: HomogenizationTest

Standard Normal Homogeneity Test (Alexandersson 1986).

Scans every candidate split m and maximises a likelihood-ratio statistic over the standardised series \(z_i = (q_i - \bar{q}) / \hat{\sigma}\):

\[T_0 = \max_{1 \le m < n} \left[\, m\, \bar{z}_1^{\,2} + (n - m)\, \bar{z}_2^{\,2} \,\right]\]

where \(\bar{z}_1 = \tfrac{1}{m}\sum_{i=1}^{m} z_i\) and \(\bar{z}_2 = \tfrac{1}{n-m}\sum_{i=m+1}^{n} z_i\) are the segment means. Uses tabulated critical values from Alexandersson & Moberg (1997) at alpha ∈ {0.01, 0.05, 0.10}.

Initialise; raise immediately if alpha is not in the critical-value table.

Source code in src/rucola/_homogeneity.py

def __init__(
    self,
    alpha: float = 0.05,
    min_years_from_end: int = 5,
    min_relative_signal: float = 1.0,
) -> None:
    """Initialise; raise immediately if alpha is not in the critical-value table."""
    super().__init__(alpha=alpha, min_years_from_end=min_years_from_end, min_relative_signal=min_relative_signal)
    if alpha not in _SNHT_CRIT:
        msg = f"SNHT: alpha must be one of {sorted(_SNHT_CRIT)}, got {alpha}"
        raise ValueError(msg)

name property

Return test name.

detect(q_series, years)

Apply SNHT to q_series and return a TestResult.

Source code in src/rucola/_homogeneity.py

def detect(self, q_series: pl.Series, years: pl.Series) -> TestResult:
    """Apply SNHT to q_series and return a TestResult."""
    rows = [(q, int(y)) for q, y in zip(q_series.to_list(), years.to_list(), strict=True) if q is not None]
    n = len(rows)
    seg_start = rows[0][1] if rows else (int(years[0]) if len(years) > 0 else 0)
    seg_end = rows[-1][1] if rows else (int(years[-1]) if len(years) > 0 else 0)
    crit = _interpolate_crit(_SNHT_CRIT[self.alpha], max(n, 10))

    null = TestResult(
        test_name=self.name,
        is_significant=False,
        break_year=seg_start,
        test_statistic=0.0,
        critical_value=crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
    )
    if n < _MIN_YEARS:
        return null

    q_vals = [r[0] for r in rows]
    y_vals = [r[1] for r in rows]

    q_mean = sum(q_vals) / n
    q_std = math.sqrt(sum((q - q_mean) ** 2 for q in q_vals) / n)
    if q_std < _NEAR_ZERO:
        return null

    z = [(q - q_mean) / q_std for q in q_vals]
    cumsum = [0.0] * (n + 1)
    for i, zi in enumerate(z):
        cumsum[i + 1] = cumsum[i] + zi

    t_vals: list[float] = []
    for m in range(1, n):
        z1 = cumsum[m] / m
        z2 = (cumsum[n] - cumsum[m]) / (n - m)
        t_vals.append(m * z1 * z1 + (n - m) * z2 * z2)

    t0 = max(t_vals)
    m_opt = t_vals.index(t0) + 1
    break_year = y_vals[m_opt]

    return TestResult(
        test_name=self.name,
        is_significant=t0 > crit,
        break_year=break_year,
        test_statistic=t0,
        critical_value=crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
        series=t_vals,
        years_tested=y_vals,
    )

rucola.BuishandTest(alpha=0.05, min_years_from_end=5, min_relative_signal=1.0)

Bases: HomogenizationTest

Buishand range test (Buishand 1982).

Built on cumulative deviations \(S_k = \sum_{i=1}^{k}(q_i - \bar{q})\). Under the null these form a bridge anchored at \(S_0 = S_n = 0\); a break pulls the bridge away from zero, inflating its range:

\[R = \frac{\max_k S_k - \min_k S_k}{\sqrt{n}\, \hat{\sigma}}\]

The break is located at \(k^\star = \arg\max_k |S_k|\). Uses tabulated critical values at alpha ∈ {0.01, 0.05, 0.10}.

Initialise; raise immediately if alpha is not in the critical-value table.

Source code in src/rucola/_homogeneity.py

def __init__(
    self,
    alpha: float = 0.05,
    min_years_from_end: int = 5,
    min_relative_signal: float = 1.0,
) -> None:
    """Initialise; raise immediately if alpha is not in the critical-value table."""
    super().__init__(alpha=alpha, min_years_from_end=min_years_from_end, min_relative_signal=min_relative_signal)
    if alpha not in _BUISHAND_CRIT:
        msg = f"Buishand: alpha must be one of {sorted(_BUISHAND_CRIT)}, got {alpha}"
        raise ValueError(msg)

name property

Return test name.

detect(q_series, years)

Apply Buishand range test to q_series and return a TestResult.

Source code in src/rucola/_homogeneity.py

def detect(self, q_series: pl.Series, years: pl.Series) -> TestResult:
    """Apply Buishand range test to q_series and return a TestResult."""
    rows = [(q, int(y)) for q, y in zip(q_series.to_list(), years.to_list(), strict=True) if q is not None]
    n = len(rows)
    seg_start = rows[0][1] if rows else (int(years[0]) if len(years) > 0 else 0)
    seg_end = rows[-1][1] if rows else (int(years[-1]) if len(years) > 0 else 0)
    crit = _interpolate_crit(_BUISHAND_CRIT[self.alpha], n)

    null = TestResult(
        test_name=self.name,
        is_significant=False,
        break_year=seg_start,
        test_statistic=0.0,
        critical_value=crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
    )
    if n < _MIN_YEARS:
        return null

    q_vals = [r[0] for r in rows]
    y_vals = [r[1] for r in rows]

    q_mean = sum(q_vals) / n
    q_std = math.sqrt(sum((q - q_mean) ** 2 for q in q_vals) / n)
    if q_std < _NEAR_ZERO:
        return null

    # cumulative deviations S_k for k = 1, …, n
    s_vals: list[float] = []
    running = 0.0
    for q in q_vals:
        running += q - q_mean
        s_vals.append(running)

    r_stat = (max(s_vals) - min(s_vals)) / (math.sqrt(n) * q_std)

    # break point: last year of pre-break segment → first year of post-break
    k_star = max(range(n), key=lambda k: abs(s_vals[k]))
    break_year = y_vals[k_star + 1] if k_star + 1 < n else y_vals[-1]

    return TestResult(
        test_name=self.name,
        is_significant=r_stat > crit,
        break_year=break_year,
        test_statistic=r_stat,
        critical_value=crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
        series=s_vals,
        years_tested=y_vals,
    )

rucola.PettittTest(alpha=0.05, min_years_from_end=5, min_relative_signal=1.0)

Bases: HomogenizationTest

Pettitt test (Pettitt 1979).

Non-parametric rank test based on the Mann–Whitney statistic. For each candidate split t it counts whether values after the split tend to be larger or smaller than values before:

\[U_{t,n} = \sum_{i=1}^{t} \sum_{j=t+1}^{n} \operatorname{sgn}(q_j - q_i), \qquad K = \max_{1 \le t < n} |U_{t,n}|\]

Because it uses signs rather than magnitudes, the test is robust to outliers and to non-normal distributions. The p-value is approximated analytically as \(p \approx 2 \exp\!\bigl(-6K^2 / (n^3 + n^2)\bigr)\), so any \(\alpha > 0\) is supported.

Source code in src/rucola/_homogeneity.py

def __init__(
    self,
    alpha: float = 0.05,
    min_years_from_end: int = 5,
    min_relative_signal: float = 1.0,
) -> None:
    """Initialise with significance level, edge-effect guard, and signal threshold.

    Parameters
    ----------
    alpha :
        Significance level for the hypothesis test (default: 0.05).
    min_years_from_end :
        Breaks within this many years of either series end are rejected
        (Hawkins 1977 edge effect, default: 5).
    min_relative_signal :
        Minimum ratio of test statistic to critical value required to accept
        a break (default: 1.0, i.e. any significant result). Raise above 1.0
        to require a stronger signal and reduce false positives.

    """
    self.alpha = alpha
    self.min_years_from_end = min_years_from_end
    self.min_relative_signal = min_relative_signal

name property

Return test name.

detect(q_series, years)

Apply Pettitt test to q_series and return a TestResult.

Source code in src/rucola/_homogeneity.py

def detect(self, q_series: pl.Series, years: pl.Series) -> TestResult:
    """Apply Pettitt test to q_series and return a TestResult."""
    rows = [(q, int(y)) for q, y in zip(q_series.to_list(), years.to_list(), strict=True) if q is not None]
    n = len(rows)
    seg_start = rows[0][1] if rows else (int(years[0]) if len(years) > 0 else 0)
    seg_end = rows[-1][1] if rows else (int(years[-1]) if len(years) > 0 else 0)
    # K_crit from p = alpha: K_crit = sqrt(−ln(α/2) · (n³+n²) / 6)
    k_crit = math.sqrt(-math.log(self.alpha / 2) * (n**3 + n**2) / 6) if n >= _MIN_YEARS else 0.0

    null = TestResult(
        test_name=self.name,
        is_significant=False,
        break_year=seg_start,
        test_statistic=0.0,
        critical_value=k_crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
    )
    if n < _MIN_YEARS:
        return null

    q_vals = [r[0] for r in rows]
    y_vals = [r[1] for r in rows]

    # U_{t,n} = Σ_{i<t} Σ_{j≥t} sgn(x_j − x_i)  for t = 1, …, n−1
    u_vals: list[float] = []
    u = 0.0
    for t in range(1, n):
        # incremental update: U_t = U_{t-1} + Σ_{j=t}^{n-1} sgn(x_j − x_{t-1})
        x_new = q_vals[t - 1]
        for j in range(t, n):
            diff = q_vals[j] - x_new
            u += 1.0 if diff > 0 else (-1.0 if diff < 0 else 0.0)
        u_vals.append(u)

    k_stat = max(abs(u) for u in u_vals)
    t_star = max(range(len(u_vals)), key=lambda t: abs(u_vals[t]))
    break_year = y_vals[t_star + 1] if t_star + 1 < n else y_vals[-1]

    return TestResult(
        test_name=self.name,
        is_significant=k_stat > k_crit,
        break_year=break_year,
        test_statistic=k_stat,
        critical_value=k_crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
        series=u_vals,
        years_tested=y_vals[1:],
    )

rucola.WorsleyTest(alpha=0.05, min_years_from_end=5, min_relative_signal=1.0)

Bases: HomogenizationTest

Worsley likelihood ratio test (Worsley 1979).

Maximum standardised two-sample t-statistic over all candidate change points \(k\):

\[W = \max_{1 \le k < n} \left| \sqrt{\frac{k\,(n-k)}{n}} \cdot \frac{\bar{q}_{1:k} - \bar{q}_{k+1:n}}{s_p} \right|\]

where \(s_p\) is the pooled standard deviation of the two segments. Critical values are computed analytically via the Bonferroni approximation \(c = \Phi^{-1}\!\bigl(1 - \alpha / (2(n-1))\bigr)\) and support any \(\alpha > 0\).

Source code in src/rucola/_homogeneity.py

def __init__(
    self,
    alpha: float = 0.05,
    min_years_from_end: int = 5,
    min_relative_signal: float = 1.0,
) -> None:
    """Initialise with significance level, edge-effect guard, and signal threshold.

    Parameters
    ----------
    alpha :
        Significance level for the hypothesis test (default: 0.05).
    min_years_from_end :
        Breaks within this many years of either series end are rejected
        (Hawkins 1977 edge effect, default: 5).
    min_relative_signal :
        Minimum ratio of test statistic to critical value required to accept
        a break (default: 1.0, i.e. any significant result). Raise above 1.0
        to require a stronger signal and reduce false positives.

    """
    self.alpha = alpha
    self.min_years_from_end = min_years_from_end
    self.min_relative_signal = min_relative_signal

name property

Return test name.

detect(q_series, years)

Apply Worsley test to q_series and return a TestResult.

Source code in src/rucola/_homogeneity.py

def detect(self, q_series: pl.Series, years: pl.Series) -> TestResult:
    """Apply Worsley test to q_series and return a TestResult."""
    rows = [(q, int(y)) for q, y in zip(q_series.to_list(), years.to_list(), strict=True) if q is not None]
    n = len(rows)
    seg_start = rows[0][1] if rows else (int(years[0]) if len(years) > 0 else 0)
    seg_end = rows[-1][1] if rows else (int(years[-1]) if len(years) > 0 else 0)
    crit = _worsley_crit(n, self.alpha) if n >= _MIN_YEARS else 0.0

    null = TestResult(
        test_name=self.name,
        is_significant=False,
        break_year=seg_start,
        test_statistic=0.0,
        critical_value=crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
    )
    if n < _MIN_YEARS:
        return null

    q_vals = [r[0] for r in rows]
    y_vals = [r[1] for r in rows]

    # O(n) running sums for means and sum-of-squares
    cum_x = [0.0] * (n + 1)
    cum_x2 = [0.0] * (n + 1)
    for i, x in enumerate(q_vals):
        cum_x[i + 1] = cum_x[i] + x
        cum_x2[i + 1] = cum_x2[i] + x * x

    total_x = cum_x[n]
    total_x2 = cum_x2[n]

    t_vals: list[float] = []
    for k in range(1, n):
        n1, n2 = k, n - k
        s1 = cum_x[k]
        s2 = total_x - s1
        ss1 = cum_x2[k] - s1 * s1 / n1
        ss2 = (total_x2 - cum_x2[k]) - s2 * s2 / n2
        pooled_var = (ss1 + ss2) / (n - 2)
        if pooled_var < _NEAR_ZERO:
            t_vals.append(0.0)
            continue
        t_k = math.sqrt(n1 * n2 / n) * (s1 / n1 - s2 / n2) / math.sqrt(pooled_var)
        t_vals.append(abs(t_k))

    w_stat = max(t_vals)
    k_star = t_vals.index(w_stat)
    break_year = y_vals[k_star + 1] if k_star + 1 < n else y_vals[-1]

    return TestResult(
        test_name=self.name,
        is_significant=w_stat > crit,
        break_year=break_year,
        test_statistic=w_stat,
        critical_value=crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
        series=t_vals,
        years_tested=y_vals[1:],
    )

rucola.EasterlingPetersonTest(alpha=0.05, min_years_from_end=5, min_relative_signal=1.0)

Bases: HomogenizationTest

Easterling–Peterson two-phase regression test (Easterling & Peterson 1995).

For each candidate break \(k\), fits a linear trend plus a step:

\[q_i = a + b\, t_i + c \cdot \mathbf{1}_{i > k} + \varepsilon_i\]

The step change \(c\) is tested via its OLS t-statistic \(t_k(c)\), and the test statistic is the maximum absolute value over all candidate breaks:

\[T = \max_{1 \le k < n} |t_k(c)|\]

Unlike mean-shift tests, the linear trend term \(b\) prevents a real trend from being mistaken for a break — useful for temperature series. Critical values use the same Bonferroni normal approximation as WorsleyTest.

Source code in src/rucola/_homogeneity.py

def __init__(
    self,
    alpha: float = 0.05,
    min_years_from_end: int = 5,
    min_relative_signal: float = 1.0,
) -> None:
    """Initialise with significance level, edge-effect guard, and signal threshold.

    Parameters
    ----------
    alpha :
        Significance level for the hypothesis test (default: 0.05).
    min_years_from_end :
        Breaks within this many years of either series end are rejected
        (Hawkins 1977 edge effect, default: 5).
    min_relative_signal :
        Minimum ratio of test statistic to critical value required to accept
        a break (default: 1.0, i.e. any significant result). Raise above 1.0
        to require a stronger signal and reduce false positives.

    """
    self.alpha = alpha
    self.min_years_from_end = min_years_from_end
    self.min_relative_signal = min_relative_signal

name property

Return test name.

detect(q_series, years)

Apply Easterling–Peterson test to q_series and return a TestResult.

Source code in src/rucola/_homogeneity.py

def detect(self, q_series: pl.Series, years: pl.Series) -> TestResult:
    """Apply Easterling–Peterson test to q_series and return a TestResult."""
    rows = [(q, int(y)) for q, y in zip(q_series.to_list(), years.to_list(), strict=True) if q is not None]
    n = len(rows)
    seg_start = rows[0][1] if rows else (int(years[0]) if len(years) > 0 else 0)
    seg_end = rows[-1][1] if rows else (int(years[-1]) if len(years) > 0 else 0)
    crit = _worsley_crit(n, self.alpha) if n >= _MIN_YEARS else 0.0

    null = TestResult(
        test_name=self.name,
        is_significant=False,
        break_year=seg_start,
        test_statistic=0.0,
        critical_value=crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
    )
    if n < _MIN_YEARS:
        return null

    q_vals = [r[0] for r in rows]
    y_vals = [r[1] for r in rows]

    # Center time axis so that X'X has zero off-diagonal in the (1,2) block,
    # which simplifies the 3×3 normal equations considerably.
    t_mean = sum(y_vals) / n
    t_c = [y - t_mean for y in y_vals]

    s_tt = sum(tc * tc for tc in t_c)
    if s_tt < _NEAR_ZERO:
        return null

    s_y = sum(q_vals)
    s_ty = sum(tc * q for tc, q in zip(t_c, q_vals, strict=True))
    s_yy = sum(q * q for q in q_vals)
    # RSS under H0 (single linear trend, no break)
    r0 = s_yy - s_y * s_y / n - s_ty * s_ty / s_tt

    # Initialise running sums for post-break segment {1, …, n-1} (k = 0)
    s_z: float = n - 1
    s_tz: float = sum(t_c[1:])
    s_zy: float = sum(q_vals[1:])

    t_stats: list[float] = []
    for k in range(n - 1):
        if k > 0:
            s_z -= 1
            s_tz -= t_c[k]
            s_zy -= q_vals[k]

        denom_c = s_z - s_z * s_z / n - s_tz * s_tz / s_tt
        if denom_c < _NEAR_ZERO:
            t_stats.append(0.0)
            continue

        numer_c = s_zy - s_z * s_y / n - s_tz * s_ty / s_tt
        rss = r0 - numer_c * numer_c / denom_c
        if rss < 0:
            rss = 0.0

        df = n - 3
        if df <= 0 or rss < _NEAR_ZERO:
            t_stats.append(0.0)
            continue

        t_k = numer_c / math.sqrt(denom_c * rss / df)
        t_stats.append(abs(t_k))

    if not t_stats:
        return null

    t_max = max(t_stats)
    k_opt = t_stats.index(t_max)
    break_year = y_vals[k_opt + 1] if k_opt + 1 < n else y_vals[-1]

    return TestResult(
        test_name=self.name,
        is_significant=t_max > crit,
        break_year=break_year,
        test_statistic=t_max,
        critical_value=crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
        series=t_stats,
        years_tested=y_vals[1:],
    )

rucola.StarsTest(l=10, alpha=0.05, min_years_from_end=5, min_relative_signal=1.0)

Bases: HomogenizationTest

STARS: Sequential T-test Analysis of Regime Shifts (Rodionov 2004).

Scans the Q-series sequentially and declares a regime shift at year j when the Regime Shift Index (RSI) stays consistently positive (or negative) over the following l years. Unlike the other tests that maximise a global statistic, STARS propagates a running regime mean, so it is sensitive to shifts that accumulate gradually across the series rather than appearing as a single dominant peak.

Variance is estimated from the mean-square successive differences (MSSD), which is robust to the very breaks being detected:

\[\hat{\sigma}^2 = \frac{1}{2(n-1)} \sum_{t=1}^{n-1}(x_{t+1} - x_t)^2\]

The minimum detectable mean shift is \(\delta = t_{\text{crit}} \cdot \hat{\sigma} \sqrt{2/l}\), where \(t_{\text{crit}}\) is the two-tailed t-distribution critical value for \(\nu = 2(l-1)\) degrees of freedom.

A shift at year j is confirmed when the RSI

\[\text{RSI}(k) = \sum_{t=j}^{k} \frac{x_t - (\bar{x}_{\text{prev}} \pm \delta)}{\hat{\sigma}\, l}\]

does not change sign for all \(k \in [j,\, j+l-1]\). If multiple shifts are confirmed, the one with the largest t-statistic is returned as the primary break year.

Parameters:

Name	Type	Description	Default
l	int	Cut-off length in years — the minimum regime duration and the RSI confirmation window. Sets \(\nu = 2(l-1)\) degrees of freedom for the t-test. Typical values for annual climate data: 10–15.	10
alpha	float	Significance level for the t-test. Supported: {0.01, 0.05, 0.10}.	0.05
min_years_from_end	int	Inherited edge-effect guard (default: 5).	5
min_relative_signal	float	Minimum ratio of the confirmed t-statistic to the critical value (default: 1.0).	1.0

References

Rodionov, S. N. (2004): A sequential algorithm for testing climate regime shifts. Geophys. Res. Lett., 31, L09204. https://doi.org/10.1029/2004GL019448

Initialise with cut-off length, significance level, and edge-guard settings.

Source code in src/rucola/_homogeneity.py

def __init__(
    self,
    l: int = 10,  # noqa: E741
    alpha: float = 0.05,
    min_years_from_end: int = 5,
    min_relative_signal: float = 1.0,
) -> None:
    """Initialise with cut-off length, significance level, and edge-guard settings."""
    super().__init__(alpha=alpha, min_years_from_end=min_years_from_end, min_relative_signal=min_relative_signal)
    self.l = l
    if alpha not in _T_CRIT:
        msg = f"StarsTest: alpha must be one of {sorted(_T_CRIT)}, got {alpha}"
        raise ValueError(msg)

name property

Return test name.

__repr__()

Return short summary string.

Source code in src/rucola/_homogeneity.py

def __repr__(self) -> str:
    """Return short summary string."""
    return (
        f"{self.__class__.__name__}("
        f"l={self.l}, "
        f"alpha={self.alpha}, "
        f"min_years_from_end={self.min_years_from_end}, "
        f"min_relative_signal={self.min_relative_signal})"
    )

detect(q_series, years)

Apply STARS to q_series and return a TestResult.

Source code in src/rucola/_homogeneity.py

def detect(self, q_series: pl.Series, years: pl.Series) -> TestResult:
    """Apply STARS to q_series and return a TestResult."""
    rows = [(q, int(y)) for q, y in zip(q_series.to_list(), years.to_list(), strict=True) if q is not None]
    n = len(rows)
    seg_start = rows[0][1] if rows else (int(years[0]) if len(years) > 0 else 0)
    seg_end = rows[-1][1] if rows else (int(years[-1]) if len(years) > 0 else 0)

    df = 2 * (self.l - 1)
    crit = _interpolate_crit(_T_CRIT[self.alpha], df) if df > 0 else float("inf")

    null = TestResult(
        test_name=self.name,
        is_significant=False,
        break_year=seg_start,
        test_statistic=0.0,
        critical_value=crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
    )

    if n < max(_MIN_YEARS, 2 * self.l):
        return null

    q_vals = [r[0] for r in rows]
    y_vals = [r[1] for r in rows]

    # σ² via mean-square successive differences — robust to the breaks we detect
    mssd = sum((q_vals[i + 1] - q_vals[i]) ** 2 for i in range(n - 1)) / (2 * (n - 1))
    if mssd < _NEAR_ZERO:
        return null
    sigma = math.sqrt(mssd)

    # Minimum detectable mean difference at this significance level
    delta = crit * sigma * math.sqrt(2.0 / self.l)

    # Sequential detection ─────────────────────────────────────────────────
    # Track the best (highest t-statistic) confirmed break.
    # After a confirmed shift at index j, the next eligible check is j+l
    # (minimum regime-length constraint).
    best_break_idx: int | None = None
    best_t: float = 0.0

    mu = sum(q_vals[: self.l]) / self.l  # regime mean, initialised over first l values
    regime_start = 0
    next_check = self.l  # first index eligible for a shift test

    for j in range(self.l, n):
        if j >= next_check and j + self.l <= n and abs(q_vals[j] - mu) >= delta:
            direction = 1 if q_vals[j] > mu else -1
            boundary = mu + direction * delta

            # Confirm via RSI: must stay on the same side for l consecutive years
            rsi = 0.0
            confirmed = True
            for k in range(j, j + self.l):
                rsi += (q_vals[k] - boundary) / (sigma * self.l)
                if direction * rsi < 0:
                    confirmed = False
                    break

            if confirmed:
                t_stat = abs(q_vals[j] - mu) / (sigma * math.sqrt(2.0 / self.l))
                if t_stat > best_t:
                    best_t = t_stat
                    best_break_idx = j
                # Start new regime: reset mean to single-point estimate and grow it
                mu = q_vals[j]
                regime_start = j
                next_check = j + self.l
                continue

        # Welford incremental update of the current regime mean
        count = j - regime_start + 1
        mu += (q_vals[j] - mu) / count

    if best_break_idx is None:
        return null

    return TestResult(
        test_name=self.name,
        is_significant=best_t > crit,
        break_year=y_vals[best_break_idx],
        test_statistic=best_t,
        critical_value=crit,
        n_years=n,
        segment_start=seg_start,
        segment_end=seg_end,
        years_tested=y_vals,
    )