Structural Break Detection in Economic Time Series

def make_breakpoint_data(
    n_steps=N_STEPS,
    break_step=BREAK_STEP,
    mean_before=MEAN_BEFORE,
    mean_after=MEAN_AFTER,
    noise_std=NOISE_STD,
    seed=SEED,
):
    """Return a Polars DataFrame with columns 'step' and 'value'.

    The series has a single structural break at break_step:
      - value[t] ~ N(mean_before, noise_std^2)  for t < break_step
      - value[t] ~ N(mean_after,  noise_std^2)  for t >= break_step

    The break is abrupt: there is no gradual transition.
    """
    rng = np.random.default_rng(seed)
    steps = np.arange(n_steps)
    means = np.where(steps < break_step, mean_before, mean_after)
    values = means + rng.normal(0.0, noise_std, n_steps)
    return pl.DataFrame({"step": steps, "value": values})

def cusum(residuals):
    """Cumulative sum (CUSUM) of residuals.

    C_t = sum_{s=0}^{t} e_s  (0-indexed).

    A structural break shifts the expected sign of residuals, so C_t drifts
    away from zero before the break and then reverses afterwards.  The index
    of maximum |C_t| estimates the location of the break.
    """
    return np.cumsum(residuals)


def detect_break(cusum_values):
    """Return the index where |CUSUM| is largest.

    This index is the last step before the detected break: the series
    properties are estimated to shift at index detect_break(...) + 1.
    """
    return int(np.argmax(np.abs(cusum_values)))

def plot_series_with_break(values, detected_break, filename):
    """Save a time-series plot with a vertical rule at the detected break."""
    n = len(values)
    df = pl.DataFrame({"step": np.arange(n, dtype=float), "value": values})
    series = (
        alt.Chart(df)
        .mark_line(color="steelblue", strokeWidth=1.5)
        .encode(
            x=alt.X("step:Q", title="Step"),
            y=alt.Y("value:Q", title="Value"),
        )
    )
    break_df = pl.DataFrame({"step": [float(detected_break + 1)]})
    break_rule = (
        alt.Chart(break_df)
        .mark_rule(color="firebrick", strokeWidth=2, strokeDash=[6, 3])
        .encode(x="step:Q")
    )
    chart = alt.layer(series, break_rule).properties(
        width=450,
        height=250,
        title=f"Time series with detected break at step {detected_break + 1}",
    )
    chart.save(filename)


def plot_cusum_comparison(cusum_mean, cusum_trend, true_break, filename):
    """Save CUSUM trajectories for both OLS models on the same axes."""
    n = len(cusum_mean)
    steps = np.arange(n, dtype=float)
    df = pl.DataFrame(
        {
            "step": np.concatenate([steps, steps]),
            "cusum": np.concatenate([cusum_mean, cusum_trend]),
            "model": ["mean-only"] * n + ["linear trend"] * n,
        }
    )
    lines = (
        alt.Chart(df)
        .mark_line(strokeWidth=2)
        .encode(
            x=alt.X("step:Q", title="Step"),
            y=alt.Y("cusum:Q", title="CUSUM"),
            color=alt.Color("model:N", legend=alt.Legend(title="OLS model")),
        )
    )
    break_df = pl.DataFrame({"step": [float(true_break)]})
    break_rule = (
        alt.Chart(break_df)
        .mark_rule(color="gray", strokeDash=[4, 4], strokeWidth=1.5)
        .encode(x="step:Q")
    )
    chart = alt.layer(lines, break_rule).properties(
        width=450,
        height=300,
        title="CUSUM comparison: mean-only vs. linear-trend residuals",
    )
    chart.save(filename)

import numpy as np
import pytest
from generate_breakpoint import make_breakpoint_data, BREAK_STEP, N_STEPS
from breakpoint import residuals_mean, residuals_trend, cusum, detect_break


def test_cusum_all_zero_residuals():
    # The CUSUM of a zero-residual sequence is identically zero.
    assert np.all(cusum(np.zeros(50)) == 0.0)


def test_cusum_matches_cumsum():
    # By definition CUSUM must equal numpy's cumsum on the same array.
    residuals = np.array([1.0, -2.0, 3.0, -0.5, 0.5])
    np.testing.assert_array_equal(cusum(residuals), np.cumsum(residuals))


def test_mean_residuals_sum_to_zero():
    # OLS residuals from a mean-only fit sum to zero because the least-squares
    # intercept equals the sample mean, centering the residuals exactly.
    df = make_breakpoint_data()
    values = df["value"].to_numpy()
    assert np.sum(residuals_mean(values)) == pytest.approx(0.0, abs=1e-10)


def test_trend_residuals_sum_to_zero():
    # OLS with an intercept always centers the residuals, so they also sum to
    # zero for the linear-trend fit.
    df = make_breakpoint_data()
    values = df["value"].to_numpy()
    assert np.sum(residuals_trend(values)) == pytest.approx(0.0, abs=1e-10)


def test_detect_break_clean_signal():
    # With no noise the CUSUM of mean-only residuals reaches its maximum
    # absolute value at index BREAK_STEP - 1 (the last step before the break),
    # because the cumulative deficit reaches its deepest point exactly there.
    values = np.concatenate(
        [
            np.zeros(BREAK_STEP),
            3.0 * np.ones(N_STEPS - BREAK_STEP),
        ]
    )
    detected = detect_break(cusum(residuals_mean(values)))
    assert detected == BREAK_STEP - 1


def test_detect_break_noisy_within_ten_steps():
    # With seed 7493418 and a signal-to-noise ratio of 3 (mean shift 3.0,
    # noise std 1.0), the detected break must be within 10 steps of the
    # true break.  detect_break returns the last index before the break,
    # so the estimated break location is detected + 1.
    df = make_breakpoint_data()
    values = df["value"].to_numpy()
    detected = detect_break(cusum(residuals_mean(values)))
    assert abs((detected + 1) - BREAK_STEP) <= 10