xtbreak Examples

Runnable Scripts

The examples/ directory contains self-contained, runnable scripts that demonstrate typical workflows:

File	Topic	Key concepts
`01_nile_river.py`	Nile River flow (classic dataset)	Single mean-shift test, estimation, CI, regime stats, Chow test
`02_us_covid_panel.py`	US COVID mortality (panel)	Panel SupF, common break dates, fixed effects, entity/time arrays
`03_gdp_great_moderation.py`	US GDP Great Moderation	Two-break detection, multi-regime comparison, volatility shifts

Run any example:

cd xtbreak-py
python3 examples/01_nile_river.py

Advanced Patterns

The following patterns go beyond the examples/ scripts.

Partial Structural Change (Common X)

The time-series xtbreak.select(..., common_x=..., num_fixed_regressors=common_x.shape[1]) surface implements partial structural change: every column of x is breaking, and every column of common_x is non-breaking. Because those non-breaking columns enter the denominator degrees of freedom and promoted estimate metadata, callers must set num_fixed_regressors=common_x.shape[1]. A zero-column common_x design is not accepted.

When some regressors break and others don’t:

import numpy as np
from xtbreak import select

time = np.arange(1, 13, dtype=float)
centered = (time - time.mean()) / 10.0
x = np.column_stack([np.ones(12), centered])        # Breaking regressors
common_x = np.column_stack([centered, centered**2]) # Non-breaking regressors
slope = np.where(time <= 6, 1.0, 3.0)
noise = np.array([0, 0.03, -0.02, 0.01, -0.01, 0.02,
                  -0.03, 0.01, 0.02, -0.02, 0.01, 0])
y = 2.0 + 0.7 * common_x[:, 0] - 1.2 * common_x[:, 1] + slope * x[:, 1] + noise
time_values = np.arange(201, 213, dtype=int)

selection = select(
    y, x,
    common_x=common_x,
    trimming=0.25,
    max_breaks=1,
    critical_values={1: 1.0},
    level_critical_values={1: 10.0},
    skip_h2=False,
    wdmax=True,
    time_values=time_values,
    num_fixed_regressors=common_x.shape[1],
)

print(f"Selected breaks: {selection.selected_breaks}")
est = selection.estimate
if est is not None:
    print(f"Break time: {est.break_table()[0]['time_value']}")
    print(f"Breaking regressors: {est.metadata['num_breaking_regressors']}")
    print(f"Fixed regressors: {est.metadata['num_fixed_regressors']}")

Stata Oracle Hydration

For parity verification against Stata xtbreak estimate output:

import json
from pathlib import Path
import xtbreak

payload = json.loads(
    Path("artifacts/verification/stata/phase1_us_ts_estimate_breaks2.json")
    .read_text(encoding="utf-8")
)
results = payload["results"]

result = xtbreak.EstimateResult.from_break_matrices(
    command=payload["command"],
    break_matrix=results["breaks"],
    ci_matrix=results["ci"],
    ssr=results["ssr"],
    metadata={"scenario_id": payload["scenario_id"]},
)

for row in result.break_table():
    print(f"Break {row['break_number']}: index={row['index']}, "
          f"time={row['time_value']}")
print(f"SSR: {result.ssr:.10f}")

SSR Path Visualization Data

For single-break estimates with SSR path data:

import json
from pathlib import Path
import xtbreak

payload = json.loads(
    Path("artifacts/verification/stata/"
         "phase1_us_ts_estimate_breaks1_ssr_path.json")
    .read_text(encoding="utf-8")
)
results = payload["results"]

result = xtbreak.EstimateResult.from_break_matrices(
    command=payload["command"],
    break_matrix=results["breaks"],
    ci_matrix=results["ci"],
    ssr=results["ssr"],
    ssr_path=results["ssr_path"],
    sample_time_values=results["sample_time_values"],
    metadata={"scenario_id": payload["scenario_id"]},
)

plot_payload = result.ssr_path_plot_data()
selected = plot_payload["selected_break"]
print(f"Selected break: index={selected['index']}, "
      f"time={selected['time_value']}")
print(f"SSR at selected: {plot_payload['points'][selected['index'] - 1]['ssr']:.10f}")

Panel Robust Covariance

Explicit-panel estimate with WPN covariance:

import numpy as np
import xtbreak

entity = np.array(["a"] * 6 + ["b"] * 6)
time = np.array([2000, 2001, 2002, 2003, 2004, 2005] * 2)
x = np.array([[-1.0], [0.0], [1.0], [-1.0], [0.0], [1.0]] * 2)
y = np.array([9.0, 10.2, 10.9, 7.1, 9.8, 13.0,
              19.1, 19.9, 21.05, 16.95, 20.2, 22.9])

result = xtbreak.estimate(
    y, x, breaks=1, trimming=2/6,
    entity=entity, time=time, vce="wpn",
)

print(f"Break: t={result.break_table()[0]['time_value']}")
print(f"Covariance: {result.metadata['covariance_estimator']}")
print(f"Bandwidth: {result.metadata['bandwidth']}")
print(f"SE: {result.metadata['covariance']['standard_errors']}")

Explicit-Panel Robust Estimate Result

Verify panel covariance metadata after estimation:

import numpy as np
import xtbreak

entity = np.array(["a"] * 6 + ["b"] * 6)
time = np.array([2000, 2001, 2002, 2003, 2004, 2005] * 2)
x = np.array([[-1.0], [0.0], [1.0], [-1.0], [0.0], [1.0]] * 2)
y = np.array([9.0, 10.2, 10.9, 7.1, 9.8, 13.0,
              19.1, 19.9, 21.05, 16.95, 20.2, 22.9])

result = xtbreak.estimate(
    y, x, breaks=1, trimming=2/6,
    entity=entity, time=time, vce="wpn",
)

assert result.metadata["covariance_estimator"] == "wpn"
assert result.metadata["covariance_runtime"] == "post_break_unrestricted_design"
assert result.metadata["bandwidth"] == 5
assert result.metadata["covariance"]["standard_errors"][0] == 0.017636701526667267
assert result.metadata["covariance"]["standard_errors"][1] == 0.022178776975932644

Quantitative Oracle Contracts

The following quantitative anchors lock examples to parity with the Stata oracle data in artifacts/verification/stata/:

US Time Series (2-break) — see phase1_us_ts_estimate_breaks2.json:

Global SSR: 49.0738641567431
Covered by: test_phase7_us_time_series_break_search.py

US Time Series (1-break SSR path) — see phase1_us_ts_estimate_breaks1_ssr_path.json:

SSR at selected break: 60.7883146046024

Panel (6-obs balanced) — see phase6_panel_manifest.json:

SupF statistic: 394.8586118251927

Monte Carlo — see test_phase7_monte_carlo_smoke.py:

False rejection rate 0.06 (nominal 5%)
Power 1.00 for large shift

Plotting (requires matplotlib)

from xtbreak import estimate, plot_breaks, plot_ssr_path, plot_regime_fit

result = estimate(y, x, breaks=1, trimming=0.15)

# Time series with break lines and CI shading
plot_breaks(result, y)

# Per-regime fitted values
plot_regime_fit(result, y, x)

# SSR objective path (single-break only)
plot_ssr_path(result)