#!/usr/bin/env python3 """Generate the public synthetic benchmark for the regime-detector note. The experiment is intentionally separate from MorphIQ Labs' private market regime study. It uses generated equal-variance AR(1) segments, canonical features, a held-out test set, and no market data or production parameters. """ from __future__ import annotations import csv import json import shutil from dataclasses import dataclass from pathlib import Path import matplotlib.pyplot as plt import numpy as np from numpy.lib.stride_tricks import sliding_window_view HERE = Path(__file__).resolve().parent REPO_CANDIDATE = HERE.parents[1] if (REPO_CANDIDATE / "package.json").exists(): REPO_ROOT: Path | None = REPO_CANDIDATE RESULTS_DIR = HERE / "results" FIGURE_DIR = REPO_ROOT / "public" / "images" / "posts" / "evaluating-regime-detectors" PUBLIC_DATA_DIR: Path | None = ( REPO_ROOT / "public" / "research" / "evaluating-regime-detectors" ) else: REPO_ROOT = None RESULTS_DIR = HERE / "results" FIGURE_DIR = HERE / "figures" PUBLIC_DATA_DIR = None WINDOW = 96 MODWT_LEVELS = 4 TRAIN_PATHS = 30 TEST_PATHS = 120 SEGMENTS_PER_PATH = 9 MIN_SEGMENT = 288 MAX_SEGMENT = 480 TRAIN_SEED = 20_260_806 TEST_SEED = 20_260_807 BOOTSTRAP_SEED = 20_260_808 BOOTSTRAP_DRAWS = 4_000 CONFIRM_SPAN = 10 CONFIRM_COUNT = 8 REGIME_NAMES = ("Persistent", "Uncorrelated", "Alternating") REGIME_PHIS = np.array([0.75, 0.0, -0.75], dtype=np.float64) METHODS = ("Rolling variance", "Lag-1 autocorrelation", "Fourier band shares", "Haar MODWT shares") COLORS = { "ink": "#0B0D10", "graphite": "#3B4048", "mute": "#8E949E", "signal": "#E85A2C", "paper": "#F4F2ED", "paper2": "#EAE7DF", "rule": "#C7C3BA", } METHOD_COLORS = { "Rolling variance": COLORS["mute"], "Lag-1 autocorrelation": COLORS["ink"], "Fourier band shares": COLORS["graphite"], "Haar MODWT shares": COLORS["signal"], } @dataclass(frozen=True) class SyntheticPath: values: np.ndarray labels: np.ndarray change_points: np.ndarray @dataclass(frozen=True) class FeatureSet: features: dict[str, np.ndarray] labels: np.ndarray endpoint_indices: np.ndarray stable_mask: np.ndarray change_points: np.ndarray @dataclass(frozen=True) class NearestCentroidModel: mean: np.ndarray scale: np.ndarray centroids: np.ndarray def predict(self, features: np.ndarray) -> np.ndarray: normalized = (features - self.mean) / self.scale distances = np.sum( (normalized[:, np.newaxis, :] - self.centroids[np.newaxis, :, :]) ** 2, axis=2, ) return np.argmin(distances, axis=1) def balanced_regime_order(rng: np.random.Generator) -> np.ndarray: """Return three shuffled blocks while avoiding repeated boundaries.""" order: list[int] = [] for _ in range(SEGMENTS_PER_PATH // len(REGIME_NAMES)): block = list(map(int, rng.permutation(len(REGIME_NAMES)))) if order and block[0] == order[-1]: swap_index = next(i for i, value in enumerate(block[1:], start=1) if value != order[-1]) block[0], block[swap_index] = block[swap_index], block[0] order.extend(block) return np.asarray(order, dtype=np.int64) def generate_path(rng: np.random.Generator) -> SyntheticPath: """Generate stationary-variance AR(1) segments with known transitions.""" order = balanced_regime_order(rng) durations = rng.integers(MIN_SEGMENT, MAX_SEGMENT + 1, size=len(order)) total = int(np.sum(durations)) values = np.empty(total, dtype=np.float64) labels = np.empty(total, dtype=np.int64) change_points: list[int] = [] offset = 0 previous = float(rng.normal()) for segment_index, (label, duration) in enumerate(zip(order, durations, strict=True)): if segment_index > 0: change_points.append(offset) phi = float(REGIME_PHIS[label]) innovation_scale = float(np.sqrt(1.0 - phi * phi)) innovations = rng.normal(size=int(duration)) for local_index, innovation in enumerate(innovations): previous = phi * previous + innovation_scale * float(innovation) values[offset + local_index] = previous labels[offset : offset + duration] = label offset += int(duration) return SyntheticPath( values=values, labels=labels, change_points=np.asarray(change_points, dtype=np.int64), ) def fourier_band_shares(centered_windows: np.ndarray) -> np.ndarray: """Return normalized power in four fixed one-sided frequency bands.""" taper = np.hanning(WINDOW) spectrum = np.abs(np.fft.rfft(centered_windows * taper, axis=1)) ** 2 frequencies = np.fft.rfftfreq(WINDOW) edges = (0.0, 1.0 / 16.0, 1.0 / 8.0, 1.0 / 4.0, 0.5000001) bands = [] for low, high in zip(edges[:-1], edges[1:], strict=True): mask = (frequencies > low) & (frequencies <= high) bands.append(np.sum(spectrum[:, mask], axis=1)) energy = np.column_stack(bands) return energy / np.maximum(np.sum(energy, axis=1, keepdims=True), np.finfo(float).tiny) def haar_modwt_shares(centered_windows: np.ndarray) -> np.ndarray: """Return four Haar MODWT detail-energy shares plus the scaling share. Haar's normalized MODWT filters are [1/2, -1/2] and [1/2, 1/2]. Filters are dilated by 2**(j-1) at level j. Periodic extension is confined to each already-observed causal window. """ scaling = centered_windows.copy() detail_energies = [] for level in range(1, MODWT_LEVELS + 1): lag = 2 ** (level - 1) shifted = np.roll(scaling, lag, axis=1) detail = 0.5 * (scaling - shifted) scaling = 0.5 * (scaling + shifted) detail_energies.append(np.mean(detail * detail, axis=1)) all_energy = np.column_stack( [*detail_energies, np.mean(scaling * scaling, axis=1)] ) return all_energy / np.maximum( np.sum(all_energy, axis=1, keepdims=True), np.finfo(float).tiny ) def extract_features(path: SyntheticPath) -> FeatureSet: windows = sliding_window_view(path.values, WINDOW) centered = windows - np.mean(windows, axis=1, keepdims=True) variance = np.mean(centered * centered, axis=1, keepdims=True) left = centered[:, :-1] right = centered[:, 1:] autocovariance = np.sum(left * right, axis=1) autocorrelation = autocovariance / np.maximum( np.sqrt(np.sum(left * left, axis=1) * np.sum(right * right, axis=1)), np.finfo(float).tiny, ) endpoint_indices = np.arange(WINDOW - 1, len(path.values), dtype=np.int64) endpoint_labels = path.labels[endpoint_indices] last_change = np.zeros(len(path.values), dtype=np.int64) current_change = 0 change_set = set(map(int, path.change_points)) for index in range(len(path.values)): if index in change_set: current_change = index last_change[index] = current_change regime_age = endpoint_indices - last_change[endpoint_indices] return FeatureSet( features={ "Rolling variance": np.log(np.maximum(variance, np.finfo(float).tiny)), "Lag-1 autocorrelation": autocorrelation[:, np.newaxis], "Fourier band shares": fourier_band_shares(centered), "Haar MODWT shares": haar_modwt_shares(centered), }, labels=endpoint_labels, endpoint_indices=endpoint_indices, stable_mask=regime_age >= WINDOW, change_points=path.change_points, ) def fit_model(features: np.ndarray, labels: np.ndarray) -> NearestCentroidModel: mean = np.mean(features, axis=0) scale = np.std(features, axis=0) scale = np.where(scale > 1e-12, scale, 1.0) normalized = (features - mean) / scale centroids = np.vstack( [np.mean(normalized[labels == label], axis=0) for label in range(len(REGIME_NAMES))] ) return NearestCentroidModel(mean=mean, scale=scale, centroids=centroids) def balanced_accuracy(truth: np.ndarray, predicted: np.ndarray) -> float: recalls = [np.mean(predicted[truth == label] == label) for label in range(len(REGIME_NAMES))] return float(np.mean(recalls)) def transition_delays(feature_set: FeatureSet, predicted: np.ndarray) -> tuple[int, list[int]]: endpoint_start = WINDOW - 1 delays: list[int] = [] eligible = 0 for change_point in feature_set.change_points: start = int(change_point) - endpoint_start if start < 0 or start >= len(predicted): continue eligible += 1 target = int(feature_set.labels[start]) stop = min(start + WINDOW + 1, len(predicted) - CONFIRM_SPAN + 1) for candidate in range(start, stop): confirmed = np.count_nonzero( predicted[candidate : candidate + CONFIRM_SPAN] == target ) if confirmed >= CONFIRM_COUNT: delays.append(int(feature_set.endpoint_indices[candidate] - change_point)) break return eligible, delays def path_metrics(feature_set: FeatureSet, predicted: np.ndarray) -> dict[str, object]: stable_truth = feature_set.labels[feature_set.stable_mask] stable_predicted = predicted[feature_set.stable_mask] eligible, delays = transition_delays(feature_set, predicted) switch_indices = np.flatnonzero(predicted[1:] != predicted[:-1]) + 1 stable_switches = int(np.count_nonzero(feature_set.stable_mask[switch_indices])) stable_count = int(np.count_nonzero(feature_set.stable_mask)) return { "balanced_accuracy": balanced_accuracy(stable_truth, stable_predicted), "eligible_transitions": eligible, "detected_transitions": len(delays), "delays": delays, "false_switches_per_1000": 1_000.0 * stable_switches / max(stable_count, 1), "all_switches_per_1000": 1_000.0 * len(switch_indices) / len(predicted), } def bootstrap_ci( values: np.ndarray, rng: np.random.Generator, statistic: str = "mean", ) -> tuple[float, float, float]: clean = values[np.isfinite(values)] if clean.size == 0: return float("nan"), float("nan"), float("nan") estimate = float(np.mean(clean) if statistic == "mean" else np.median(clean)) samples = rng.choice(clean, size=(BOOTSTRAP_DRAWS, clean.size), replace=True) draws = np.mean(samples, axis=1) if statistic == "mean" else np.median(samples, axis=1) low, high = np.quantile(draws, [0.025, 0.975]) return estimate, float(low), float(high) def summarize( metrics: dict[str, list[dict[str, object]]], bootstrap_rng: np.random.Generator, ) -> list[dict[str, float | str]]: rows: list[dict[str, float | str]] = [] for method in METHODS: method_metrics = metrics[method] accuracy = np.asarray([m["balanced_accuracy"] for m in method_metrics], dtype=float) hit_rates = np.asarray( [m["detected_transitions"] / m["eligible_transitions"] for m in method_metrics], dtype=float, ) false_switches = np.asarray( [m["false_switches_per_1000"] for m in method_metrics], dtype=float ) all_switches = np.asarray( [m["all_switches_per_1000"] for m in method_metrics], dtype=float ) delays = np.asarray( [delay for m in method_metrics for delay in m["delays"]], dtype=float ) acc, acc_low, acc_high = bootstrap_ci(accuracy, bootstrap_rng) hit, hit_low, hit_high = bootstrap_ci(hit_rates, bootstrap_rng) delay, delay_low, delay_high = bootstrap_ci(delays, bootstrap_rng, "median") false, false_low, false_high = bootstrap_ci(false_switches, bootstrap_rng) switches, switches_low, switches_high = bootstrap_ci(all_switches, bootstrap_rng) rows.append( { "method": method, "balanced_accuracy": acc, "balanced_accuracy_ci_low": acc_low, "balanced_accuracy_ci_high": acc_high, "transition_hit_rate": hit, "transition_hit_rate_ci_low": hit_low, "transition_hit_rate_ci_high": hit_high, "median_delay": delay, "median_delay_ci_low": delay_low, "median_delay_ci_high": delay_high, "false_switches_per_1000": false, "false_switches_per_1000_ci_low": false_low, "false_switches_per_1000_ci_high": false_high, "all_switches_per_1000": switches, "all_switches_per_1000_ci_low": switches_low, "all_switches_per_1000_ci_high": switches_high, } ) return rows def write_summary(rows: list[dict[str, float | str]]) -> None: metadata = { "experiment": "Equal-variance AR(1) spectral-regime benchmark", "window": WINDOW, "modwt_levels": MODWT_LEVELS, "train_paths": TRAIN_PATHS, "test_paths": TEST_PATHS, "segments_per_path": SEGMENTS_PER_PATH, "segment_length_range": [MIN_SEGMENT, MAX_SEGMENT], "regimes": dict(zip(REGIME_NAMES, map(float, REGIME_PHIS), strict=True)), "train_seed": TRAIN_SEED, "test_seed": TEST_SEED, "bootstrap_seed": BOOTSTRAP_SEED, "bootstrap_draws": BOOTSTRAP_DRAWS, "transition_confirmation": { "span": CONFIRM_SPAN, "required_target_labels": CONFIRM_COUNT, }, "results": rows, } destinations = [RESULTS_DIR] if PUBLIC_DATA_DIR is not None: destinations.append(PUBLIC_DATA_DIR) for destination in destinations: destination.mkdir(parents=True, exist_ok=True) with (destination / "summary.csv").open("w", newline="", encoding="utf-8") as handle: writer = csv.DictWriter( handle, fieldnames=list(rows[0].keys()), lineterminator="\n", ) writer.writeheader() writer.writerows(rows) (destination / "summary.json").write_text( json.dumps(metadata, indent=2) + "\n", encoding="utf-8" ) def configure_matplotlib() -> None: plt.rcParams.update( { "figure.facecolor": COLORS["paper"], "axes.facecolor": COLORS["paper"], "savefig.facecolor": COLORS["paper"], "text.color": COLORS["ink"], "axes.labelcolor": COLORS["graphite"], "axes.edgecolor": COLORS["rule"], "xtick.color": COLORS["graphite"], "ytick.color": COLORS["graphite"], "font.family": "sans-serif", "font.size": 10, "axes.titleweight": "bold", "axes.titlelocation": "left", } ) def shade_regimes(axis: plt.Axes, path: SyntheticPath) -> None: starts = np.concatenate(([0], path.change_points)) stops = np.concatenate((path.change_points, [len(path.values)])) shades = ("#F2C6B7", "#E0DDD5", "#BEC1C7") for start, stop in zip(starts, stops, strict=True): label = int(path.labels[int(start)]) axis.axvspan(start, stop, color=shades[label], alpha=0.32, linewidth=0) def plot_representative(path: SyntheticPath, feature_set: FeatureSet) -> None: configure_matplotlib() x = feature_set.endpoint_indices variance = np.exp(feature_set.features["Rolling variance"][:, 0]) autocorrelation = feature_set.features["Lag-1 autocorrelation"][:, 0] fourier = feature_set.features["Fourier band shares"] modwt = feature_set.features["Haar MODWT shares"] fourier_contrast = fourier[:, 0] - fourier[:, -1] modwt_contrast = (modwt[:, -1] + modwt[:, -2]) - modwt[:, 0] fig, axes = plt.subplots(4, 1, figsize=(12, 8.6), sharex=True) series = ( (variance, "Rolling variance", 1.0), (autocorrelation, "Lag-1 autocorrelation", 0.0), (fourier_contrast, "Fourier low-minus-high band share", 0.0), (modwt_contrast, "Haar MODWT coarse-minus-fine share", 0.0), ) for axis, (values, label, reference) in zip(axes, series, strict=True): shade_regimes(axis, path) axis.plot(x, values, color=COLORS["ink"], linewidth=0.9) axis.axhline(reference, color=COLORS["mute"], linewidth=0.7, linestyle="--") axis.set_title(label, fontsize=10, fontweight="bold", pad=5) axis.grid(axis="y", color=COLORS["rule"], linewidth=0.5, alpha=0.65) axis.spines[["top", "right"]].set_visible(False) fig.suptitle( "One held-out path: identical marginal variance, different dependence", x=0.055, ha="left", fontsize=15, fontweight="bold", ) axes[-1].set_xlabel("Observation") fig.text( 0.995, 0.01, "Background: persistent (warm), uncorrelated (light), alternating (dark)", ha="right", color=COLORS["graphite"], fontsize=9, ) fig.tight_layout(rect=(0, 0.025, 1, 0.965), h_pad=1.35) FIGURE_DIR.mkdir(parents=True, exist_ok=True) fig.savefig(FIGURE_DIR / "held-out-path.png", dpi=180, bbox_inches="tight") plt.close(fig) def plot_summary(rows: list[dict[str, float | str]]) -> None: configure_matplotlib() fig, axes = plt.subplots(2, 2, figsize=(12, 8.5)) specifications = ( ("balanced_accuracy", "Stable balanced accuracy", lambda value: 100.0 * value, "%"), ("transition_hit_rate", "Transitions detected within one window", lambda value: 100.0 * value, "%"), ("median_delay", "Median confirmed detection delay", float, "observations"), ("false_switches_per_1000", "Off-transition label switches", float, "per 1,000 observations"), ) for axis, (key, title, transform, unit) in zip(axes.flat, specifications, strict=True): values = np.asarray([transform(float(row[key])) for row in rows]) low = np.asarray([transform(float(row[f"{key}_ci_low"])) for row in rows]) high = np.asarray([transform(float(row[f"{key}_ci_high"])) for row in rows]) errors = np.vstack((values - low, high - values)) positions = np.arange(len(rows)) axis.barh( positions, values, color=[METHOD_COLORS[str(row["method"])] for row in rows], height=0.58, xerr=errors, error_kw={"ecolor": COLORS["ink"], "elinewidth": 0.8, "capsize": 2}, ) axis.set_yticks(positions, [str(row["method"]) for row in rows]) axis.invert_yaxis() axis.set_title(title) axis.set_xlabel(unit) axis.grid(axis="x", color=COLORS["rule"], linewidth=0.5, alpha=0.65) axis.set_axisbelow(True) axis.spines[["top", "right", "left"]].set_visible(False) if key == "balanced_accuracy": axis.axvline(100.0 / len(REGIME_NAMES), color=COLORS["mute"], linestyle="--", linewidth=0.8) fig.suptitle("Held-out benchmark results", x=0.06, ha="left", fontsize=15, fontweight="bold") fig.tight_layout(rect=(0, 0, 1, 0.96), h_pad=2.4, w_pad=3.0) FIGURE_DIR.mkdir(parents=True, exist_ok=True) fig.savefig(FIGURE_DIR / "benchmark-summary.png", dpi=180, bbox_inches="tight") plt.close(fig) def main() -> None: train_rng = np.random.default_rng(TRAIN_SEED) test_rng = np.random.default_rng(TEST_SEED) bootstrap_rng = np.random.default_rng(BOOTSTRAP_SEED) train_sets = [extract_features(generate_path(train_rng)) for _ in range(TRAIN_PATHS)] models: dict[str, NearestCentroidModel] = {} for method in METHODS: train_features = np.vstack( [feature_set.features[method][feature_set.stable_mask] for feature_set in train_sets] ) train_labels = np.concatenate( [feature_set.labels[feature_set.stable_mask] for feature_set in train_sets] ) models[method] = fit_model(train_features, train_labels) test_pairs = [ (path, extract_features(path)) for path in (generate_path(test_rng) for _ in range(TEST_PATHS)) ] metrics: dict[str, list[dict[str, object]]] = {method: [] for method in METHODS} for _, feature_set in test_pairs: for method in METHODS: predicted = models[method].predict(feature_set.features[method]) metrics[method].append(path_metrics(feature_set, predicted)) rows = summarize(metrics, bootstrap_rng) write_summary(rows) plot_representative(*test_pairs[0]) plot_summary(rows) if PUBLIC_DATA_DIR is not None: PUBLIC_DATA_DIR.mkdir(parents=True, exist_ok=True) shutil.copy2(__file__, PUBLIC_DATA_DIR / "benchmark.py") for row in rows: print( f"{row['method']}: " f"accuracy={float(row['balanced_accuracy']):.3f}, " f"hit={float(row['transition_hit_rate']):.3f}, " f"delay={float(row['median_delay']):.1f}, " f"false_switches={float(row['false_switches_per_1000']):.1f}/1000" ) if __name__ == "__main__": main()