""" ou_process_backtest_v1.py =========================== Tests whether Ornstein-Uhlenbeck process parameter estimation can improve the validated EURUSD model. Unlike Markov/Hurst/Kalman tests, OU is theoretically justified here — the yield spread residual IS an OU process by construction. This tests whether we can exploit that structure more precisely than the fixed validated parameters. Three approaches tested against the validated 4,019 trade baseline: APPROACH A — Dynamic Threshold Estimate θ (reversion speed) rolling. Derive the theoretically optimal z-score entry level. Fast reversion (low θ) → lower threshold. Slow reversion (high θ) → higher threshold. APPROACH B — Dynamic Hold Time Set max hold = 2× estimated OU half-life (bounded 8-80h). If spread reverts fast, exit sooner. If slow, hold longer. APPROACH C — OU Fit Quality Filter Only trade when the spread series is well-described by OU (mean-reversion coefficient clearly > 0). Reject signals when the spread has broken OU behaviour (trending/non-stationary). APPROACH D — Combined All three together. Key question: Does OU-adaptive behaviour improve Sharpe by ≥ 0.3 while keeping ≥ 80% of trades? Uses trades_real_costs.csv (all 4,019 validated trades) — no re-simulation. OU parameters are applied as filters/adjustments ON TOP of existing trades. Run from project root: python src/research/ou_process_backtest_v1.py """ import pandas as pd import numpy as np from pathlib import Path import sys import warnings warnings.filterwarnings('ignore') BASE_PATH = Path(__file__).resolve().parents[2] SRC_PATH = BASE_PATH / "src" PROC_PATH = BASE_PATH / "data" / "processed" if str(SRC_PATH) not in sys.path: sys.path.insert(0, str(SRC_PATH)) # ── Validated model constants ───────────────────────────────────────────────── THRESHOLD_BASE = 2.75 HOLD_BASE = 52 # hours NOTIONAL = 300_000 DT = 1/24 # 1 hour in days OU_WINDOW = 120 # bars of z-score history for OU estimation ZSCORE_BANDS = [(2.75, 3.50, 1.0), (3.50, 4.50, 1.5), (4.50, 99.0, 2.0)] def get_mult(z_abs): for lo, hi, m in ZSCORE_BANDS: if lo <= z_abs < hi: return m return 2.0 # ── OU Parameter Estimation ─────────────────────────────────────────────────── def estimate_ou_params(series, dt=DT): """ Estimate OU parameters via OLS on discrete AR(1) representation. The OU process: dX = θ(μ - X)dt + σdW Discretised: X_{t+1} = a + b*X_t + ε Where: b = exp(-θ*dt), a = μ*(1-b) Returns dict with: theta : speed of mean reversion (per day) mu : long-run mean sigma : diffusion coefficient half_life : ln(2)/θ in days half_life_h: half-life in hours quality : mean-reversion strength (1-b), higher = stronger OU """ x = np.array(series, dtype=float) if len(x) < 20: return None x_lag = x[:-1] x_cur = x[1:] # OLS regression b, a = np.polyfit(x_lag, x_cur, 1) # Validate OU conditions: 0 < b < 1 (mean-reverting, not explosive) if b <= 0 or b >= 1.0: return None theta = -np.log(b) / dt mu = a / (1 - b) residuals = x_cur - (a + b * x_lag) sigma_sq = np.var(residuals) * 2 * theta / (1 - b**2) sigma = np.sqrt(max(sigma_sq, 1e-10)) half_life = np.log(2) / theta # days return { 'theta' : float(theta), 'mu' : float(mu), 'sigma' : float(sigma), 'half_life' : float(half_life), 'half_life_h': float(half_life * 24), 'b' : float(b), 'quality' : float(np.clip(1 - b, 0, 1)), } def expected_reversion_pct(theta, hold_hours): """Fraction of dislocation expected to resolve in hold_hours.""" return float(1 - np.exp(-theta * hold_hours / 24)) def dynamic_threshold(theta, base=THRESHOLD_BASE): """ Derive optimal entry threshold from OU speed of reversion. Fast reversion (θ > 2.0/day, HL < 8h): lower threshold viable Medium reversion (θ 0.5-2.0/day): near baseline Slow reversion (θ < 0.5/day, HL > 33h): higher threshold needed Uses continuous scaling: thresh = base * (1 + 0.5 * (1 - rev_pct)) where rev_pct = expected fraction resolved in 52h window """ rev_pct = expected_reversion_pct(theta, HOLD_BASE) # Scale: 100% reversion expected → threshold = base * 1.0 # 50% reversion expected → threshold = base * 1.25 # 25% reversion expected → threshold = base * 1.375 thresh = base * (1 + 0.5 * (1 - rev_pct)) return float(np.clip(thresh, 2.0, 4.5)) def dynamic_hold(theta, min_h=8, max_h=80): """Set hold time as 2× half-life, bounded by min/max hours.""" half_life_h = np.log(2) / theta * 24 return int(np.clip(2 * half_life_h, min_h, max_h)) # ── Build OU parameter history ──────────────────────────────────────────────── def build_ou_history(z_series, window=OU_WINDOW): """ Compute rolling OU parameters for every timestamp in z_series. Returns DataFrame with ou_theta, ou_mu, ou_sigma, ou_hl_h, ou_quality. """ print(f" Computing rolling OU parameters (window={window} bars)...") records = [] z_vals = z_series.values z_idx = z_series.index n = len(z_vals) for i in range(window, n): window_data = z_vals[i-window:i] params = estimate_ou_params(window_data) if params: records.append({ 'datetime' : z_idx[i], 'ou_theta' : params['theta'], 'ou_mu' : params['mu'], 'ou_sigma' : params['sigma'], 'ou_hl_h' : params['half_life_h'], 'ou_quality' : params['quality'], 'ou_b' : params['b'], }) else: records.append({ 'datetime' : z_idx[i], 'ou_theta' : np.nan, 'ou_mu' : np.nan, 'ou_sigma' : np.nan, 'ou_hl_h' : np.nan, 'ou_quality' : np.nan, 'ou_b' : np.nan, }) df = pd.DataFrame(records).set_index('datetime') valid = df['ou_theta'].notna().sum() print(f" Valid OU estimates: {valid:,} / {len(df):,} " f"({valid/len(df)*100:.1f}%)") print(f" θ median={df['ou_theta'].median():.3f} " f"HL median={df['ou_hl_h'].median():.1f}h") return df # ── Apply OU parameters to trades ───────────────────────────────────────────── def get_ou_at_entry(ou_history, entry_time): """Get OU parameters closest to trade entry time.""" if ou_history is None or ou_history.empty: return None idx = ou_history.index.searchsorted(entry_time) if idx == 0: return None row = ou_history.iloc[idx - 1] if pd.isna(row['ou_theta']): return None return row def apply_ou_approaches(trades_df, ou_history): """ Apply all four OU approaches to the validated trade set. Returns dict of results for each approach. """ records = [] for _, trade in trades_df.iterrows(): entry_time = pd.to_datetime(trade['entry_time']) z_abs = float(trade.get('zscore_abs', THRESHOLD_BASE)) pnl = float(trade['dollar_pnl_real']) win = int(pnl > 0) exit_r = trade.get('exit_reason', 'unknown') ou = get_ou_at_entry(ou_history, entry_time) # Default OU params if not available if ou is None: theta = 1.0 hl_h = 16.7 quality = 0.05 else: theta = float(ou['ou_theta']) hl_h = float(ou['ou_hl_h']) quality = float(ou['ou_quality']) # Derived OU values dyn_thresh = dynamic_threshold(theta) dyn_hold = dynamic_hold(theta) rev_pct = expected_reversion_pct(theta, HOLD_BASE) * 100 # Approach A: Dynamic threshold filter # Trade only if |z| >= dyn_thresh (OU-derived) take_a = z_abs >= dyn_thresh # Approach B: Dynamic hold time # We cannot re-simulate here, but we can approximate: # If dyn_hold < actual hold and trade was a hold_expiry exit, # it might have been cut short (could be better or worse). # For now: keep all trades but flag those where hold mismatch is large hold_ratio = dyn_hold / HOLD_BASE # Scale P&L for very short suggested holds (proxy for early exit) pnl_b = pnl * min(hold_ratio, 1.0) if exit_r == 'hold_expiry' else pnl take_b = True # keep all trades, adjust hold_expiry P&L # Approach C: OU quality filter # Only trade when mean-reversion is clearly present take_c = quality >= 0.04 # threshold: b <= 0.96 # Approach D: Combined (A + C) take_d = take_a and take_c records.append({ 'entry_time' : entry_time, 'z_abs' : z_abs, 'pnl' : pnl, 'win' : win, 'exit_reason': exit_r, 'ou_theta' : theta, 'ou_hl_h' : hl_h, 'ou_quality' : quality, 'dyn_thresh' : dyn_thresh, 'dyn_hold' : dyn_hold, 'rev_pct_52h': rev_pct, 'pnl_b' : pnl_b, 'take_a' : take_a, 'take_b' : take_b, 'take_c' : take_c, 'take_d' : take_d, }) return pd.DataFrame(records) # ── Scoring ─────────────────────────────────────────────────────────────────── def score(df, pnl_col, label): if len(df) == 0: return {"label": label, "n": 0, "wr_pct": 0, "sharpe": 0, "pf": 0, "max_dd": 0, "total": 0, "per_yr": 0} pnl = df[pnl_col].values n = len(pnl) wr = (pnl > 0).mean() * 100 total = pnl.sum() eq = 100_000 + np.cumsum(pnl) peak = np.maximum.accumulate(eq) max_dd = ((eq - peak) / peak * 100).min() per_yr = n / 22 excess = pnl / 100_000 - 0.04 / per_yr sharpe = (excess.mean() / excess.std() * np.sqrt(per_yr) if excess.std() > 0 else 0) gp = pnl[pnl > 0].sum() gl = abs(pnl[pnl < 0].sum()) pf = round(gp / gl, 2) if gl > 0 else 999 return {"label": label, "n": n, "wr_pct": round(wr, 1), "sharpe": round(sharpe, 2), "pf": pf, "max_dd": round(max_dd, 2), "total": round(total), "per_yr": round(per_yr, 1)} def print_results(results, title): print(f"\n{'='*72}") print(title) print(f"{'='*72}") print(f" {'Approach':<42}{'N':>5}{'WR%':>6}{'Sharpe':>8}" f"{'PF':>6}{'MaxDD%':>8}{'Total$':>10}{'Tr/yr':>7}") print(f" {'─'*72}") base_sharpe = max((r.get('sharpe', 0) for r in results if 'baseline' in r.get('label','').lower()), default=0) for r in results: if r.get('n', 0) == 0: continue gain = r['sharpe'] - base_sharpe marker = f" ◄ +{gain:.2f}" if gain > 0.05 and 'baseline' not in r['label'].lower() else "" print(f" {r['label']:<42}{r['n']:>5}{r['wr_pct']:>6.1f}" f"{r['sharpe']:>8.2f}{r['pf']:>6.2f}{r['max_dd']:>8.2f}" f"{r['total']:>10,.0f}{r['per_yr']:>7.1f}{marker}") # ── Main ────────────────────────────────────────────────────────────────────── def main(): print("=" * 72) print("ORNSTEIN-UHLENBECK PROCESS BACKTEST v1") print("Testing OU-adaptive parameters vs validated EURUSD model") print("=" * 72) # ── Load validated trades ───────────────────────────────────────────────── trades_path = next((p for p in [ PROC_PATH / "trades_real_costs.csv", PROC_PATH / "trades" / "trades_real_costs.csv", ] if p.exists()), None) if not trades_path: print("ERROR: trades_real_costs.csv not found"); sys.exit(1) trades = pd.read_csv(trades_path) trades["entry_time"] = pd.to_datetime(trades["entry_time"]) print(f"\nLoaded {len(trades):,} validated trades") # ── Load z-score series ─────────────────────────────────────────────────── print("\nLoading z-score series...") try: from features.spot_lag_v3 import get_model_ready_spot_lag_v3 sig_df = get_model_ready_spot_lag_v3().copy() sig_df["datetime"] = pd.to_datetime(sig_df["datetime"]) sig_df = sig_df.sort_values("datetime").reset_index(drop=True) z_series = sig_df.set_index("datetime")["lag_zscore_24h_v3"] print(f" Loaded {len(z_series):,} hourly z-score readings") print(f" Date range: {z_series.index[0].date()} → " f"{z_series.index[-1].date()}") except Exception as e: print(f" ERROR: {e}"); sys.exit(1) # ── Compute rolling OU parameters ───────────────────────────────────────── print("\nStep 1: Estimating rolling OU parameters...") ou_history = build_ou_history(z_series, window=OU_WINDOW) # Summary statistics valid = ou_history.dropna(subset=['ou_theta']) print(f"\n OU Parameter Distribution:") for col, label in [('ou_theta', 'θ (rev speed/day)'), ('ou_hl_h', 'Half-life (hours)'), ('ou_quality', 'Quality (1-b)')]: v = valid[col] print(f" {label:<22}: " f"p10={v.quantile(0.1):.3f} " f"median={v.median():.3f} " f"p90={v.quantile(0.9):.3f}") # ── Apply OU to trades ──────────────────────────────────────────────────── print("\nStep 2: Applying OU parameters to validated trades...") df = apply_ou_approaches(trades, ou_history) # Add z-score band multiplier for dollar scaling df['mult'] = df['z_abs'].apply(get_mult) # ── Results ─────────────────────────────────────────────────────────────── baseline = score(df, 'pnl', "Baseline (validated)") # Approach A: Dynamic threshold df_a = df[df['take_a']] res_a = score(df_a, 'pnl', f"A: Dynamic OU threshold (avg {df['dyn_thresh'].mean():.2f})") # Approach A variants: strict and loose df_a2 = df[df['z_abs'] >= df['dyn_thresh'] * 0.9] # 10% looser res_a2 = score(df_a2, 'pnl', "A2: Dynamic threshold -10% (looser)") # Approach B: Dynamic hold time (approximated via P&L adjustment) res_b = score(df, 'pnl_b', f"B: Dynamic hold time (avg {df['dyn_hold'].mean():.0f}h)") # Approach C: OU quality filter df_c = df[df['take_c']] q_thresh = df['ou_quality'].quantile(0.25) df_c2 = df[df['ou_quality'] >= q_thresh] res_c = score(df_c, 'pnl', f"C: OU quality filter (b<0.96)") res_c2 = score(df_c2, 'pnl', f"C2: OU quality top 75% (quality>{q_thresh:.3f})") # Approach D: Combined A + C df_d = df[df['take_d']] res_d = score(df_d, 'pnl', "D: Combined (dyn thresh + quality)") # Approach E: Only trade when half-life < hold time # (Only take trades where OU says reversion should complete in time) df_e = df[df['ou_hl_h'] < HOLD_BASE] res_e = score(df_e, 'pnl', f"E: Only trade when HL < {HOLD_BASE}h (expected to resolve)") # Approach F: Regime-based — high θ (fast reversion) only theta_median = df['ou_theta'].median() df_f = df[df['ou_theta'] >= theta_median] res_f = score(df_f, 'pnl', f"F: Fast reversion regime (θ≥{theta_median:.2f})") all_results = [baseline, res_a, res_a2, res_b, res_c, res_c2, res_d, res_e, res_f] print_results(all_results, "OU PROCESS BACKTEST — ALL APPROACHES") # ── Detailed OU diagnostics ─────────────────────────────────────────────── print(f"\n{'─'*72}") print("OU PARAMETER DIAGNOSTICS AT TRADE ENTRY") print(f"{'─'*72}") print(f"\n {'Metric':<30}{'Win trades':>12}{'Loss trades':>13}") print(f" {'─'*55}") wins = df[df['win'] == 1] losses = df[df['win'] == 0] for col, label in [('ou_theta', 'θ (reversion speed)'), ('ou_hl_h', 'Half-life (hours)'), ('ou_quality', 'OU quality (1-b)'), ('dyn_thresh', 'Derived threshold'), ('rev_pct_52h','Expected reversion%')]: w_med = wins[col].median() l_med = losses[col].median() diff = w_med - l_med flag = " ◄" if abs(diff) > 0.05 * abs(w_med + l_med + 0.01) else "" print(f" {label:<30}{w_med:>12.3f}{l_med:>13.3f}{flag}") print(f"\n Trades where OU threshold > 2.75 (model says wait):") stricter = (df['dyn_thresh'] > THRESHOLD_BASE + 0.1).sum() print(f" {stricter:,} trades ({stricter/len(df)*100:.1f}%)") print(f"\n Trades where OU threshold < 2.75 (model says enter earlier):") looser = (df['dyn_thresh'] < THRESHOLD_BASE - 0.1).sum() print(f" {looser:,} trades ({looser/len(df)*100:.1f}%)") print(f"\n OU-suggested hold time distribution:") h = df['dyn_hold'] for lo, hi, label in [(8, 24, "Very fast <24h"), (24, 52, "Fast 24-52h"), (52, 80, "Slow >52h")]: cnt = ((h >= lo) & (h < hi)).sum() print(f" {label:<20}: {cnt:>5,} trades ({cnt/len(df)*100:.1f}%)") # ── Half-life vs win rate ───────────────────────────────────────────────── print(f"\n Win rate by half-life bucket:") df['hl_bucket'] = pd.cut(df['ou_hl_h'], bins=[0, 8, 16, 24, 40, 60, 200], labels=['0-8h','8-16h','16-24h', '24-40h','40-60h','60h+']) hl_stats = df.groupby('hl_bucket').agg( n=('win','count'), wr=('win',lambda x: f"{x.mean()*100:.1f}%") ) print(f" {'HL bucket':<12}{'N':>6}{'WR':>8}") for bucket, row in hl_stats.iterrows(): print(f" {str(bucket):<12}{row['n']:>6,}{row['wr']:>8}") # ── Conclusion ──────────────────────────────────────────────────────────── best = max(all_results, key=lambda r: r.get('sharpe', 0)) base_sh = baseline['sharpe'] sharpe_gain = best['sharpe'] - base_sh trade_pct = best['n'] / baseline['n'] * 100 print(f"\n{'='*72}") print("CONCLUSION") print(f"{'='*72}") print(f"\n Baseline (validated) : {baseline['n']:,} trades " f"WR {baseline['wr_pct']:.1f}% Sharpe {base_sh:.2f}") print(f" Best OU approach : {best['n']:,} trades " f"WR {best['wr_pct']:.1f}% Sharpe {best['sharpe']:.2f}") print(f" Best approach label : {best['label']}") print(f" Sharpe gain : {sharpe_gain:+.2f}") print(f" Trades retained : {trade_pct:.1f}%") if sharpe_gain >= 0.3 and trade_pct >= 80: verdict = "IMPLEMENT — OU adaptation genuinely improves the model" detail = (f"Sharpe +{sharpe_gain:.2f}, " f"keeps {trade_pct:.0f}% of trades") elif sharpe_gain >= 0.1: verdict = "MARGINAL — improvement exists but below implementation threshold" detail = "Sharpe gain < 0.3, not worth added complexity" else: verdict = "NO IMPROVEMENT — fixed parameters already optimal for this signal" detail = ("The z-score of the yield spread is well-described by OU,\n" " but the fixed ±2.75/52h parameters already capture the\n" " optimal trade-off empirically.") print(f"\n VERDICT: {verdict}") print(f" DETAIL : {detail}") print(f""" THEORETICAL CONTEXT ──────────────────────────────────────────────────────────────────── The yield spread residual IS an OU process by construction. The z-score normalises it to zero mean — θ describes how quickly extreme dislocations pull back to equilibrium. If OU parameters add no value it means either: (a) The fixed 2.75/52h already aligns well with the OU dynamics (b) θ varies too much across 22 years to be estimated reliably (c) The Fibonacci entry timing already selects for fast-reverting signals regardless of θ Check the half-life table above — if wins and losses have similar θ values, the OU dynamics are consistent across all trades and the fixed parameters are genuinely optimal. If wins have significantly higher θ, there may be a regime-based threshold worth pursuing. ──────────────────────────────────────────────────────────────────── """) out = PROC_PATH / "ou_process_results.csv" df.to_csv(out, index=False) print(f"Trade-level OU classifications saved: {out.name}") if __name__ == "__main__": main()