r""" historical_breach_validation.py ================================ Counts ACTUAL threshold breaches across the full 23-year backtest, using the exact z-score computation the live monitor uses. Produces a definitive month-by-month table answering: "How often did the model fire historically?" Methodology: EURUSD: rolling-beta lag-gap z-score on canonical CSV combination (matches live monitor's OLS path exactly, BETA_WINDOW=120, SMOOTH_SPAN=20, EU_Z_WINDOW=60, threshold ±2.75, session 07-17 EET) USDJPY: spread-level z-score on rate CSVs (matches live monitor: UJ_Z_WINDOW=30 days, threshold ±2.0, session 07-17 EET) A "breach" is defined as a calendar day where at least one in-session bar satisfied |z| >= threshold. This matches how the model fires — at most one entry per day per pair regardless of how many hours z stayed above threshold. Run on home machine (takes ~2 minutes — runs OLS on 5700+ daily windows): cd C:\Users\paul_\OneDrive\fx_macro_intraday python src\research\historical_breach_validation.py Output: 1. Monthly breach count table for both pairs (printed) 2. Annual summary statistics 3. CSV file at data/research/historical_breaches.csv for further analysis Author: Bennett's FX Quant project. Added 2026-05-12. """ from __future__ import annotations import sys import warnings from pathlib import Path import numpy as np import pandas as pd warnings.filterwarnings("ignore") # ────────────────────────────────────────────────────────────────────── # Paths # ────────────────────────────────────────────────────────────────────── BASE = Path(__file__).resolve().parents[2] RATES_DIR = BASE / "data" / "raw" / "rates" PRICE_DIR = BASE / "data" / "raw" / "prices" / "EURUSD" OUT_DIR = BASE / "data" / "research" OUT_DIR.mkdir(parents=True, exist_ok=True) OUT_CSV = OUT_DIR / "historical_breaches.csv" # ────────────────────────────────────────────────────────────────────── # Model parameters (must match live monitor exactly) # ────────────────────────────────────────────────────────────────────── EU_THRESHOLD = 2.75 UJ_THRESHOLD = 2.00 BETA_WINDOW = 120 # daily bars for OLS SMOOTH_SPAN = 20 # EWM span for beta smoothing EU_Z_WINDOW = 60 # hourly bars for EU z-score UJ_Z_WINDOW = 30 # daily bars for UJ z-score # Session is 07:00-17:00 EET. Live monitor uses (UTC + 2) % 24 for EET. # So in-session UTC hours are 5-15 (inclusive). SESSION_UTC_START_HOUR = 5 SESSION_UTC_END_HOUR = 15 # ────────────────────────────────────────────────────────────────────── # Data loading (mirrors live monitor and zscore_recent_check.py) # ────────────────────────────────────────────────────────────────────── def load_rate(filename: str, col: str) -> pd.DataFrame: path = RATES_DIR / filename df = pd.read_csv(path) df.columns = [c.lower() for c in df.columns] df["date"] = pd.to_datetime(df.iloc[:, 0], errors="coerce") df[col] = pd.to_numeric(df.iloc[:, 1], errors="coerce") return df[["date", col]].dropna().sort_values("date").reset_index(drop=True) def load_eurusd_canonical() -> pd.DataFrame: """Load combined hourly EURUSD price data (matches live monitor).""" archive = PRICE_DIR / "EURUSD_1H_2003_2026.csv" recent = PRICE_DIR / "EURUSD_1H_recent.csv" TIMESTAMP_CANDIDATES = [ "datetime", "date", "time", "timestamp", "time (eet)", "time(eet)", "time_eet", "time eet", ] frames = [] for path in [archive, recent]: if not path.exists(): continue df = pd.read_csv(path) df.columns = [c.strip().lower() for c in df.columns] dt_col = next((c for c in TIMESTAMP_CANDIDATES if c in df.columns), None) if dt_col is None: print(f" WARNING: no datetime column in {path.name}, columns: {list(df.columns)}") continue df["datetime"] = pd.to_datetime(df[dt_col], format="mixed", errors="coerce") df["close"] = pd.to_numeric(df["close"], errors="coerce") df = df.dropna(subset=["datetime", "close"]) frames.append(df[["datetime", "close"]]) if not frames: raise RuntimeError("No price data could be loaded") combined = (pd.concat(frames, ignore_index=True) .drop_duplicates(subset=["datetime"], keep="last") .sort_values("datetime") .reset_index(drop=True)) return combined # ────────────────────────────────────────────────────────────────────── # EURUSD z-score computation (matches _compute_eurusd_zscore in live monitor) # ────────────────────────────────────────────────────────────────────── def compute_eurusd_zscore_series(prices_1h: pd.DataFrame) -> pd.DataFrame: """Compute EURUSD z-score for every hourly bar across history. Returns DataFrame with columns: datetime, close, zscore, in_session """ import statsmodels.api as sm prices_1h = prices_1h.copy() prices_1h["date"] = prices_1h["datetime"].dt.normalize() # Daily close daily_px = (prices_1h.groupby("date", as_index=False) .agg(close=("close", "last")) .sort_values("date").reset_index(drop=True)) daily_px["eurusd_return_1d"] = daily_px["close"].pct_change() daily_px["date"] = pd.to_datetime(daily_px["date"]) # Load rates us2y = load_rate("us2y.csv", "us2y") us10y = load_rate("us10y.csv", "us10y") de2y = load_rate("de2y.csv", "de2y") de10y = load_rate("de10y.csv", "de10y") # Build daily spread series with forward-fill date_range = pd.DataFrame({"date": pd.date_range( daily_px["date"].min(), daily_px["date"].max(), freq="D")}) def asof(rate_df, col): left = date_range.copy().astype({"date": "datetime64[us]"}) right = rate_df.copy().astype({"date": "datetime64[us]"}) return pd.merge_asof(left, right, on="date", direction="backward") spreads = (asof(us2y, "us2y") .merge(asof(us10y, "us10y"), on="date") .merge(asof(de2y, "de2y"), on="date") .merge(asof(de10y, "de10y"), on="date")) spreads["spread_2y"] = spreads["us2y"] - spreads["de2y"] spreads["spread_10y"] = spreads["us10y"] - spreads["de10y"] spreads["spread_2y_change_1d"] = spreads["spread_2y"].diff(1) spreads["spread_10y_change_1d"] = spreads["spread_10y"].diff(1) spreads = spreads[["date", "spread_2y_change_1d", "spread_10y_change_1d"]].dropna() model_df = (daily_px.merge(spreads, on="date", how="left") .dropna(subset=["eurusd_return_1d", "spread_2y_change_1d"]) .reset_index(drop=True)) print(f" Running OLS across {len(model_df)} daily bars (this takes ~60s)...") n = len(model_df) b2y_raw = np.full(n, np.nan) b10y_raw = np.full(n, np.nan) for i in range(BETA_WINDOW, n): sample = model_df.iloc[i - BETA_WINDOW : i].copy() sample = sample[(sample["spread_2y_change_1d"].abs() > 0) | (sample["spread_10y_change_1d"].abs() > 0)] if len(sample) < 30: continue X = sm.add_constant(sample[["spread_2y_change_1d", "spread_10y_change_1d"]]) try: res = sm.OLS(sample["eurusd_return_1d"], X).fit() b2y_raw[i] = res.params.get("spread_2y_change_1d", np.nan) b10y_raw[i] = res.params.get("spread_10y_change_1d", np.nan) except Exception: pass model_df["beta_2y"] = pd.Series(b2y_raw ).ewm(span=SMOOTH_SPAN).mean().shift(1).clip(-0.10, 0.10).values model_df["beta_10y"] = pd.Series(b10y_raw).ewm(span=SMOOTH_SPAN).mean().shift(1).clip(-0.08, 0.08).values model_df["predicted_return_1d"] = ( model_df["beta_2y"] * model_df["spread_2y_change_1d"] + model_df["beta_10y"] * model_df["spread_10y_change_1d"] ) # Merge back to hourly + compute z-score daily_pred = model_df[["date", "predicted_return_1d"]] df = prices_1h.merge(daily_pred, on="date", how="left") df["predicted_return_1d"] = df["predicted_return_1d"].ffill() df["eurusd_return_24h"] = df["close"].pct_change(24) df["lag_gap_24h"] = df["predicted_return_1d"] - df["eurusd_return_24h"] mean = df["lag_gap_24h"].rolling(EU_Z_WINDOW, min_periods=EU_Z_WINDOW).mean() std = df["lag_gap_24h"].rolling(EU_Z_WINDOW, min_periods=EU_Z_WINDOW).std() df["zscore"] = (df["lag_gap_24h"] - mean) / std df["utc_hour"] = df["datetime"].dt.hour df["in_session"] = ((df["utc_hour"] >= SESSION_UTC_START_HOUR) & (df["utc_hour"] <= SESSION_UTC_END_HOUR)) return df[["datetime", "close", "zscore", "in_session"]].dropna(subset=["zscore"]) # ────────────────────────────────────────────────────────────────────── # USDJPY z-score computation (matches _compute_spread_zscore) # ────────────────────────────────────────────────────────────────────── def compute_usdjpy_zscore_series() -> pd.DataFrame: """Compute UJ z-score for every day. Returns: date, spread, zscore.""" us2y = load_rate("us2y.csv", "us2y") jp2y = load_rate("jp2y.csv", "jp2y") rates = pd.merge(us2y, jp2y, on="date", how="outer").sort_values("date") rates = rates.set_index("date") rates = rates.reindex( pd.date_range(rates.index.min(), rates.index.max(), freq="D") ).ffill().dropna() rates["spread"] = rates["us2y"] - rates["jp2y"] rates["z_mean"] = rates["spread"].rolling(UJ_Z_WINDOW, min_periods=UJ_Z_WINDOW).mean() rates["z_std"] = rates["spread"].rolling(UJ_Z_WINDOW, min_periods=UJ_Z_WINDOW).std() rates["zscore"] = (rates["spread"] - rates["z_mean"]) / rates["z_std"] out = rates[["spread", "zscore"]].dropna(subset=["zscore"]).reset_index() out = out.rename(columns={"index": "date"}) return out # ────────────────────────────────────────────────────────────────────── # Breach counting — at most 1 per day per pair # ────────────────────────────────────────────────────────────────────── def count_eurusd_breaches_per_day(eu_z: pd.DataFrame) -> pd.DataFrame: """Returns DataFrame: date, eu_max_abs_z, eu_breach (0/1).""" df = eu_z.copy() df = df[df["in_session"]] df["date"] = df["datetime"].dt.normalize() daily = df.groupby("date").agg(eu_max_abs_z=("zscore", lambda s: s.abs().max())).reset_index() daily["eu_breach"] = (daily["eu_max_abs_z"] >= EU_THRESHOLD).astype(int) return daily def count_usdjpy_breaches_per_day(uj_z: pd.DataFrame) -> pd.DataFrame: """Returns DataFrame: date, uj_abs_z, uj_breach (0/1). UJ uses daily spread so 1 row per day already. """ df = uj_z.copy() df["uj_abs_z"] = df["zscore"].abs() df["uj_breach"] = (df["uj_abs_z"] >= UJ_THRESHOLD).astype(int) # UJ session filter: only count if it's a weekday (model only trades Mon-Fri) df["weekday"] = df["date"].dt.dayofweek df.loc[df["weekday"] >= 5, "uj_breach"] = 0 return df[["date", "uj_abs_z", "uj_breach"]] # ────────────────────────────────────────────────────────────────────── # Main # ────────────────────────────────────────────────────────────────────── def main(): print("=" * 72) print(" HISTORICAL BREACH VALIDATION — 23-year backtest") print("=" * 72) print() print("Step 1/4: Loading EURUSD canonical price data...") prices = load_eurusd_canonical() print(f" Loaded {len(prices):,} hourly bars | {prices['datetime'].min()} → {prices['datetime'].max()}") print() print("Step 2/4: Computing EURUSD z-score series across history...") eu_z = compute_eurusd_zscore_series(prices) print(f" Computed z-score for {len(eu_z):,} hourly bars") print() print("Step 3/4: Computing USDJPY z-score series across history...") uj_z = compute_usdjpy_zscore_series() print(f" Computed z-score for {len(uj_z):,} daily bars") print() print("Step 4/4: Counting per-day breaches and aggregating by month...") eu_daily = count_eurusd_breaches_per_day(eu_z) uj_daily = count_usdjpy_breaches_per_day(uj_z) daily = eu_daily.merge(uj_daily, on="date", how="outer").fillna(0) daily["eu_breach"] = daily["eu_breach"].astype(int) daily["uj_breach"] = daily["uj_breach"].astype(int) daily["total_breach"] = daily["eu_breach"] + daily["uj_breach"] daily = daily.sort_values("date").reset_index(drop=True) daily.to_csv(OUT_DIR / "historical_breaches_daily.csv", index=False) # Monthly aggregation daily["year_month"] = daily["date"].dt.to_period("M") monthly = (daily.groupby("year_month") .agg(eu_breaches=("eu_breach", "sum"), uj_breaches=("uj_breach", "sum"), total_breaches=("total_breach", "sum"), trading_days=("date", "count")) .reset_index()) monthly.to_csv(OUT_CSV, index=False) print(f" Saved {OUT_CSV}") print() # Display results print("=" * 72) print(" MONTHLY BREACH COUNTS") print("=" * 72) print() print(f" {'Month':<10} {'EU':>4} {'UJ':>4} {'Tot':>5} {'Days':>5}") print(" " + "-" * 32) for _, row in monthly.iterrows(): print(f" {str(row['year_month']):<10} {row['eu_breaches']:>4} {row['uj_breaches']:>4} " f"{row['total_breaches']:>5} {row['trading_days']:>5}") print() print("=" * 72) print(" ANNUAL SUMMARY") print("=" * 72) print() monthly["year"] = monthly["year_month"].dt.year annual = (monthly.groupby("year") .agg(eu=("eu_breaches", "sum"), uj=("uj_breaches", "sum"), tot=("total_breaches", "sum"), months=("year_month", "count")) .reset_index()) annual["eu_per_month"] = (annual["eu"] / annual["months"]).round(2) annual["uj_per_month"] = (annual["uj"] / annual["months"]).round(2) annual["tot_per_month"] = (annual["tot"] / annual["months"]).round(2) annual["tot_per_week"] = (annual["tot"] / annual["months"] / 4.33).round(2) print(f" {'Year':<6} {'EU':>4} {'UJ':>4} {'Total':>6} {'EU/mo':>7} {'UJ/mo':>7} {'Tot/mo':>8} {'Tot/wk':>8}") print(" " + "-" * 56) for _, row in annual.iterrows(): print(f" {row['year']:<6} {row['eu']:>4} {row['uj']:>4} {row['tot']:>6} " f"{row['eu_per_month']:>7.2f} {row['uj_per_month']:>7.2f} " f"{row['tot_per_month']:>8.2f} {row['tot_per_week']:>8.2f}") print() print("=" * 72) print(" OVERALL STATISTICS") print("=" * 72) print() total_months = len(monthly) total_eu = monthly["eu_breaches"].sum() total_uj = monthly["uj_breaches"].sum() total_all = monthly["total_breaches"].sum() print(f" Total months in study: {total_months}") print(f" Total EU breaches: {total_eu} ({total_eu/total_months:.2f}/month, {total_eu/total_months/4.33:.2f}/week)") print(f" Total UJ breaches: {total_uj} ({total_uj/total_months:.2f}/month, {total_uj/total_months/4.33:.2f}/week)") print(f" Total combined breaches: {total_all} ({total_all/total_months:.2f}/month, {total_all/total_months/4.33:.2f}/week)") print() print(f" Percentiles of monthly combined breach count:") for p in [10, 25, 50, 75, 90, 95, 99]: v = np.percentile(monthly["total_breaches"], p) print(f" p{p:>2}: {v:.1f}") print() print(f" Most active month: {monthly.loc[monthly['total_breaches'].idxmax(), 'year_month']} " f"with {monthly['total_breaches'].max()} breaches") print() print(f" Output written to: {OUT_CSV}") if __name__ == "__main__": main()