""" markov_regime_backtest_v1.py ============================== Tests whether Markov and Hidden Markov regime classification improves the validated EURUSD model by: 1. SIMPLE MARKOV: Classify each trade's regime based on z-score persistence in the N bars before entry. Filter or reduce size on counter-trend trades in high-persistence regimes. 2. HIDDEN MARKOV MODEL (HMM): Use hmmlearn to fit a 2-state HMM on the z-score series (State 0 = low volatility/mean-reverting, State 1 = trending/persistent dislocation). Only trade when HMM is in the trending state. Key hypothesis from live trading observation: Trade 1: LONG, z=-3.2, TP hit cleanly Trade 2: SHORT next day, z=-3.0 (COUNTER to prior regime) → SL hit Regime was still LONG-biased — counter-trend trade was low quality Run from project root: python src/research/markov_regime_backtest_v1.py Requires: pip install hmmlearn """ 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)) # ── Model constants ─────────────────────────────────────────────────────────── NOTIONAL = 300_000 ZSCORE_BANDS = [(2.75, 3.50, 1.0), (3.50, 4.50, 1.5), (4.50, 99.0, 2.0)] AVG_MULT = 0.71*1.0 + 0.20*1.5 + 0.09*2.0 def get_mult(z_abs): for lo, hi, m in ZSCORE_BANDS: if lo <= z_abs < hi: return m return 2.0 # ── Regime classification ───────────────────────────────────────────────────── def classify_regime_markov(z_series, entry_time, lookback_bars=4): """ Simple Markov regime classifier. Looks at the last N z-score readings before entry and classifies: HIGH_PERSISTENCE_LONG — 3+ of last 4 bars strongly above 0 HIGH_PERSISTENCE_SHORT — 3+ of last 4 bars strongly below 0 CHOPPY — alternating, no clear direction NEUTRAL — z-scores near zero Returns: regime, confidence (0-1), direction_count """ # Get z-scores in the lookback window before entry mask = (z_series.index < entry_time) prior = z_series[mask].tail(lookback_bars) if len(prior) < 2: return "NEUTRAL", 0.5, 0 vals = prior.values # Count directional readings long_bars = (vals > 1.0).sum() short_bars = (vals < -1.0).sum() strong_long = (vals > 2.0).sum() strong_short = (vals < -2.0).sum() n = len(vals) confidence = max(long_bars, short_bars) / n if strong_long >= 3 or (long_bars >= 3 and n >= 4): return "HIGH_PERSIST_LONG", confidence, long_bars elif strong_short >= 3 or (short_bars >= 3 and n >= 4): return "HIGH_PERSIST_SHORT", confidence, short_bars elif long_bars + short_bars < n * 0.5: return "NEUTRAL", confidence, 0 else: # Check for alternating pattern (choppy) sign_changes = sum(1 for i in range(1, len(vals)) if np.sign(vals[i]) != np.sign(vals[i-1])) if sign_changes >= 2: return "CHOPPY", 1 - confidence, sign_changes return "NEUTRAL", confidence, 0 def is_counter_trend(regime, trade_direction): """ Returns True if trade direction conflicts with the current regime. trade_direction: 1=LONG, -1=SHORT """ if regime == "HIGH_PERSIST_LONG" and trade_direction == -1: return True if regime == "HIGH_PERSIST_SHORT" and trade_direction == 1: return True return False class SimpleGaussianHMM: """ 2-state Gaussian HMM — pure numpy, no hmmlearn required. Uses Baum-Welch EM for fitting and Viterbi for decoding. """ def __init__(self, n_states=2, n_iter=100, random_state=42): self.n_states = n_states self.n_iter = n_iter np.random.seed(random_state) def fit(self, X): X = X.flatten(); n = len(X); K = self.n_states idxs = np.argsort(X) self.means_ = np.array([X[idxs[:n//2]].mean(), X[idxs[n//2:]].mean()]) self.covars_ = np.array([X[idxs[:n//2]].var()+1e-6, X[idxs[n//2:]].var()+1e-6]) self.startprob_ = np.ones(K)/K self.transmat_ = np.full((K,K), 0.1/(K-1)) np.fill_diagonal(self.transmat_, 0.9) for _ in range(self.n_iter): le = self._log_emit(X) alpha, ls = self._forward(le) beta = self._backward(le, ls) gamma, xi = self._gamma_xi(alpha, beta, le) self.startprob_ = gamma[0]+1e-10; self.startprob_ /= self.startprob_.sum() self.transmat_ = xi.sum(0)+1e-10; self.transmat_ /= self.transmat_.sum(1,keepdims=True) w = gamma.sum(0)+1e-10 self.means_ = (gamma*X[:,None]).sum(0)/w self.covars_ = np.maximum((gamma*(X[:,None]-self.means_)**2).sum(0)/w, 1e-6) return self def predict(self, X): return self._viterbi(self._log_emit(X.flatten())) def _log_emit(self, X): out = np.zeros((len(X), self.n_states)) for k in range(self.n_states): out[:,k] = -0.5*np.log(2*np.pi*self.covars_[k]) - 0.5*(X-self.means_[k])**2/self.covars_[k] return out def _forward(self, le): n,K = le.shape; la = np.zeros((n,K)); ls = np.zeros(n) la[0] = np.log(self.startprob_+1e-300)+le[0] ls[0] = np.logaddexp.reduce(la[0]); la[0] -= ls[0] lt = np.log(self.transmat_+1e-300) for t in range(1,n): for k in range(K): la[t,k] = np.logaddexp.reduce(la[t-1]+lt[:,k])+le[t,k] ls[t] = np.logaddexp.reduce(la[t]); la[t] -= ls[t] return la, ls def _backward(self, le, ls): n,K = le.shape; lb = np.zeros((n,K)) lt = np.log(self.transmat_+1e-300) for t in range(n-2,-1,-1): for k in range(K): lb[t,k] = np.logaddexp.reduce(lt[k]+le[t+1]+lb[t+1]) lb[t] -= ls[t+1] return lb def _gamma_xi(self, la, lb, le): n,K = la.shape lg = la+lb; lg -= np.logaddexp.reduce(lg,1,keepdims=True); gamma = np.exp(lg) lt = np.log(self.transmat_+1e-300); xi = np.zeros((n-1,K,K)) for t in range(n-1): for i in range(K): for j in range(K): xi[t,i,j] = la[t,i]+lt[i,j]+le[t+1,j]+lb[t+1,j] xi[t] = np.exp(xi[t]-np.logaddexp.reduce(xi[t].flatten())) return gamma, xi def _viterbi(self, le): n,K = le.shape; v = np.zeros((n,K)); psi = np.zeros((n,K),dtype=int) v[0] = np.log(self.startprob_+1e-300)+le[0] lt = np.log(self.transmat_+1e-300) for t in range(1,n): for k in range(K): tp = v[t-1]+lt[:,k]; psi[t,k] = np.argmax(tp); v[t,k] = tp[psi[t,k]]+le[t,k] s = np.zeros(n,dtype=int); s[-1] = np.argmax(v[-1]) for t in range(n-2,-1,-1): s[t] = psi[t+1,s[t+1]] return s def fit_hmm(z_series, n_states=2): """Fit SimpleGaussianHMM to z-score series. No external dependencies.""" X = z_series.values model = SimpleGaussianHMM(n_states=n_states, n_iter=80, random_state=42) model.fit(X) states = model.predict(X) means = [abs(model.means_[s]) for s in range(n_states)] trending_state = int(np.argmax(means)) return model, states, trending_state # ── Main backtest ───────────────────────────────────────────────────────────── def run_markov_backtest(trades_df, z_series): """ Apply Markov regime filter to each trade. Tests multiple approaches: 1. Baseline (no filter) 2. Skip counter-trend trades in high-persistence regime 3. Reduce size (0.5x) on counter-trend trades 4. Skip choppy regime trades """ results = {} trade_details = [] for _, trade in trades_df.iterrows(): entry_time = pd.to_datetime(trade["entry_time"]) direction = int(trade["signal"]) z_abs = abs(float(trade.get("zscore_abs", 2.75))) pnl = float(trade["dollar_pnl_real"]) ret = float(trade["return_real"]) # Classify regime regime, conf, dir_count = classify_regime_markov( z_series, entry_time, lookback_bars=4) counter = is_counter_trend(regime, direction) mult = get_mult(z_abs) trade_details.append({ "entry_time" : entry_time, "direction" : direction, "regime" : regime, "confidence" : round(conf, 3), "counter_trend": counter, "pnl_base" : pnl, "ret" : ret, "mult" : mult, "win" : int(pnl > 0), }) df = pd.DataFrame(trade_details) # ── Approach 1: Baseline ────────────────────────────────────────────────── results["1_baseline"] = score(df, df["pnl_base"], "Baseline (no filter)") # ── Approach 2: Skip counter-trend in high-persistence ──────────────────── mask2 = ~((df["counter_trend"]) & (df["regime"].isin(["HIGH_PERSIST_LONG","HIGH_PERSIST_SHORT"]))) results["2_skip_counter"] = score( df[mask2], df[mask2]["pnl_base"], "Skip counter-trend in high-persistence regime") # ── Approach 3: Half size on counter-trend ──────────────────────────────── pnl3 = df["pnl_base"].copy() counter_mask = (df["counter_trend"]) & \ (df["regime"].isin(["HIGH_PERSIST_LONG","HIGH_PERSIST_SHORT"])) pnl3[counter_mask] = pnl3[counter_mask] * 0.5 results["3_half_size_counter"] = score(df, pnl3, "Half size on counter-trend trades") # ── Approach 4: Skip choppy regime ─────────────────────────────────────── mask4 = df["regime"] != "CHOPPY" results["4_skip_choppy"] = score( df[mask4], df[mask4]["pnl_base"], "Skip all choppy regime trades") # ── Approach 5: Only high persistence trades ────────────────────────────── mask5 = df["regime"].isin(["HIGH_PERSIST_LONG","HIGH_PERSIST_SHORT"]) results["5_persist_only"] = score( df[mask5], df[mask5]["pnl_base"], "Only high-persistence regime trades") return results, df def run_hmm_backtest(trades_df, z_series): """Apply Hidden Markov Model regime filter.""" print("\nFitting HMM to z-score series (pure numpy, no external deps)...") model, states, trending_state = fit_hmm(z_series) if model is None: return None, None # Create state lookup state_lookup = pd.Series(states, index=z_series.index) # Label states means = model.means_.flatten() state_labels = {} for s in range(2): abs_mean = abs(means[s]) if abs_mean > 1.0: state_labels[s] = "TRENDING" else: state_labels[s] = "MEAN_REV" print(f" State 0: mean={means[0]:.3f} → {state_labels[0]}") print(f" State 1: mean={means[1]:.3f} → {state_labels[1]}") print(f" Transition matrix:\n{model.transmat_.round(3)}") # Assign HMM state to each trade trade_states = [] for _, trade in trades_df.iterrows(): entry_time = pd.to_datetime(trade["entry_time"]) # Find closest z-score reading idx = state_lookup.index.searchsorted(entry_time) if idx > 0: state = int(state_lookup.iloc[idx-1]) else: state = 0 trade_states.append({ "hmm_state" : state, "hmm_label" : state_labels[state], "pnl" : float(trade["dollar_pnl_real"]), "win" : int(float(trade["dollar_pnl_real"]) > 0), "entry_time" : entry_time, }) df_hmm = pd.DataFrame(trade_states) hmm_results = {} # Baseline hmm_results["hmm_1_baseline"] = score( df_hmm, df_hmm["pnl"], "HMM Baseline") # Only trade in TRENDING state mask_t = df_hmm["hmm_label"] == "TRENDING" hmm_results["hmm_2_trending_only"] = score( df_hmm[mask_t], df_hmm[mask_t]["pnl"], "HMM: Only trade in TRENDING state") # Only trade in MEAN_REV state mask_m = df_hmm["hmm_label"] == "MEAN_REV" hmm_results["hmm_3_meanrev_only"] = score( df_hmm[mask_m], df_hmm[mask_m]["pnl"], "HMM: Only trade in MEAN_REVERTING state") return hmm_results, df_hmm def score(df, pnl_series, label): """Calculate key performance metrics for a trade set.""" if len(df) == 0: return {"label": label, "n": 0} pnl = pnl_series.values n = len(pnl) wins = (pnl > 0).sum() wr = wins / n * 100 total = pnl.sum() # Equity curve equity = 100_000 + np.cumsum(pnl) peak = np.maximum.accumulate(equity) dd = ((equity - peak) / peak * 100) max_dd = dd.min() # Sharpe years = n / 182 # ~182 trades per year rf = 0.04 per_yr = n / 22 # per year excess = pnl / 100_000 - rf / per_yr sharpe = (excess.mean() / excess.std() * np.sqrt(per_yr) if excess.std() > 0 else 0) # Profit factor gp = pnl[pnl > 0].sum() gl = abs(pnl[pnl < 0].sum()) pf = gp / gl if gl > 0 else 999 return { "label" : label, "n" : n, "wr_pct" : round(wr, 1), "total" : round(total), "sharpe" : round(sharpe, 2), "pf" : round(pf, 2), "max_dd" : round(max_dd, 2), "per_yr" : round(n / 22, 1), } def print_results(results, title="RESULTS"): print(f"\n{'='*70}") print(f"{title}") print(f"{'='*70}") print(f"{'Approach':<40}{'N':>6}{'WR%':>7}{'Sharpe':>8}" f"{'PF':>7}{'MaxDD%':>8}{'Total$':>10}{'Trades/yr':>10}") print(f"{'─'*70}") for k, r in sorted(results.items()): if r.get("n", 0) == 0: continue print(f" {r['label']:<38}{r['n']:>6}{r['wr_pct']:>7.1f}" f"{r['sharpe']:>8.2f}{r['pf']:>7.2f}" f"{r['max_dd']:>8.2f}{r['total']:>10,.0f}{r['per_yr']:>10.1f}") # ── Brownian Motion / Hurst Exponent Analysis ───────────────────────────────── def hurst_dfa(series, min_lag=10, max_lag=None): """ Hurst exponent via Detrended Fluctuation Analysis (DFA). H=0.5 Brownian Motion (random walk), H>0.5 persistent, H<0.5 anti-persistent """ ts = np.array(series, dtype=float) n = len(ts) if max_lag is None: max_lag = n // 4 Y = np.cumsum(ts - np.mean(ts)) lags = np.unique(np.logspace(np.log10(min_lag), np.log10(max_lag), 20).astype(int)) lags = lags[lags >= min_lag] flucts, valid_lags = [], [] for lag in lags: segs = len(Y) // lag if segs < 2: continue rms_list = [] for i in range(segs): seg = Y[i*lag:(i+1)*lag] x = np.arange(len(seg)) trend = np.polyval(np.polyfit(x, seg, 1), x) rms_list.append(np.sqrt(np.mean((seg - trend)**2))) if rms_list: flucts.append(np.log(np.mean(rms_list))) valid_lags.append(np.log(lag)) if len(valid_lags) < 2: return 0.5 return round(float(np.clip(np.polyfit(valid_lags, flucts, 1)[0], 0.0, 1.5)), 4) def hurst_regime_scale(H): if H >= 0.65: return 1.25, "HIGH_PERSISTENCE" if H >= 0.55: return 1.00, "MILD_PERSISTENCE" if H >= 0.45: return 1.00, "RANDOM_WALK" return 0.75, "ANTI_PERSISTENT" def run_brownian_backtest(trades_df, z_series): """Rolling Hurst exponent regime backtest (120-bar lookback).""" print("\nCalculating rolling Hurst exponents (120-bar window)...") LOOKBACK = 120 rows = [] for i, (_, trade) in enumerate(trades_df.iterrows()): if i % 500 == 0: print(f" {i/len(trades_df)*100:.0f}% — {i} trades processed") entry_time = pd.to_datetime(trade["entry_time"]) pnl = float(trade["dollar_pnl_real"]) window = z_series[z_series.index < entry_time].tail(LOOKBACK) H = hurst_dfa(window.values, min_lag=5, max_lag=max(10, len(window)//3)) \ if len(window) >= 30 else 0.5 scale, label = hurst_regime_scale(H) rows.append({"entry_time": entry_time, "H": H, "regime": label, "risk_scale": scale, "pnl_base": pnl, "pnl_scaled": pnl * scale, "win": int(pnl > 0)}) df = pd.DataFrame(rows) bm_results = {} bm_results["bm_1_baseline"] = score(df, df["pnl_base"], "BM Baseline (no scaling)") bm_results["bm_2_hurst_scaled"] = score(df, df["pnl_scaled"], "BM Hurst-scaled risk") mask_p = df["H"] >= 0.55 bm_results["bm_3_persist_only"] = score(df[mask_p], df[mask_p]["pnl_base"], "BM Only H>=0.55 (persistent)") mask_np = df["H"] >= 0.45 bm_results["bm_4_skip_antipersist"] = score(df[mask_np], df[mask_np]["pnl_base"], "BM Skip anti-persistent H<0.45") print("\n Regime distribution:") for r in ["HIGH_PERSISTENCE","MILD_PERSISTENCE","RANDOM_WALK","ANTI_PERSISTENT"]: sub = df[df["regime"]==r] if len(sub): print(f" {r:<20}: {len(sub):>4} trades " f"WR={sub['win'].mean()*100:.1f}% H_avg={sub['H'].mean():.3f}") print(f"\n Overall H: mean={df['H'].mean():.3f} " f"min={df['H'].min():.3f} max={df['H'].max():.3f}") return bm_results, df def main(): print("=" * 70) print("MARKOV REGIME BACKTEST v1") print("Testing regime-based filtering on validated 4,019 EURUSD trades") print("=" * 70) # ── Load validated trades ───────────────────────────────────────────────── trades_path = PROC_PATH / "trades_real_costs.csv" if not trades_path.exists(): trades_path = PROC_PATH / "trades" / "trades_real_costs.csv" if not trades_path.exists(): 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") print(f"Columns: {list(trades.columns)}") # Find z-score column z_col = next((c for c in trades.columns if "zscore" in c.lower()), None) if z_col: print(f"Z-score column: {z_col}") else: print("WARNING: No zscore column — using 2.75 for all") trades["zscore_abs"] = 2.75 # ── Load z-score series for regime classification ───────────────────────── print("\nLoading z-score series for regime classification...") try: from features.spot_lag_v3 import get_model_ready_spot_lag_v3 sig_df = get_model_ready_spot_lag_v3() sig_df["datetime"] = pd.to_datetime(sig_df["datetime"]) z_series = sig_df.set_index("datetime")["lag_zscore_24h_v3"] print(f" Loaded {len(z_series):,} hourly z-score readings") except Exception as e: print(f" ERROR: {e}") sys.exit(1) # ── Run Markov backtest ─────────────────────────────────────────────────── print("\nRunning Simple Markov regime classification...") markov_results, trade_df = run_markov_backtest(trades, z_series) print_results(markov_results, "SIMPLE MARKOV REGIME FILTER RESULTS") # ── Regime breakdown ────────────────────────────────────────────────────── print(f"\n{'─'*70}") print("REGIME BREAKDOWN") print(f"{'─'*70}") regime_stats = trade_df.groupby("regime").agg( trades = ("pnl_base", "count"), win_rate = ("win", lambda x: f"{x.mean()*100:.1f}%"), total_pnl = ("pnl_base", lambda x: f"${x.sum():,.0f}"), counter = ("counter_trend", "sum"), ).sort_values("trades", ascending=False) print(regime_stats.to_string()) print(f"\n{'─'*70}") print("COUNTER-TREND TRADES IN HIGH-PERSISTENCE REGIME") counter_df = trade_df[ (trade_df["counter_trend"]) & (trade_df["regime"].isin(["HIGH_PERSIST_LONG","HIGH_PERSIST_SHORT"])) ] if len(counter_df) > 0: wr = counter_df["win"].mean() * 100 tot = counter_df["pnl_base"].sum() print(f" Count : {len(counter_df):,}") print(f" Win rate : {wr:.1f}% (vs {trade_df['win'].mean()*100:.1f}% baseline)") print(f" Total PnL: ${tot:,.0f}") print(f" Per year : {len(counter_df)/22:.1f}") if wr < trade_df["win"].mean() * 100 - 5: print(f"\n ✓ Counter-trend trades significantly WORSE than baseline") print(f" → Evidence supports Markov regime filter") else: print(f"\n ✗ Counter-trend trades not significantly worse than baseline") print(f" → Markov regime filter unlikely to add value") else: print(" No counter-trend trades found in dataset") # ── Run HMM backtest ────────────────────────────────────────────────────── print(f"\n{'='*70}") print("HIDDEN MARKOV MODEL (HMM) REGIME FILTER") print(f"{'='*70}") hmm_results, hmm_df = run_hmm_backtest(trades, z_series) if hmm_results: print_results(hmm_results, "HMM REGIME FILTER RESULTS") # ── Conclusion ──────────────────────────────────────────────────────────── print(f"\n{'='*70}") print("CONCLUSION") print(f"{'='*70}") base_wr = markov_results["1_baseline"]["wr_pct"] base_sh = markov_results["1_baseline"]["sharpe"] best_key = max(markov_results, key=lambda k: markov_results[k].get("sharpe", 0)) best = markov_results[best_key] if best["sharpe"] > base_sh + 0.2 and best["wr_pct"] > base_wr: rec = f"IMPLEMENT: {best['label']}" detail = (f"Sharpe improves {base_sh:.2f}→{best['sharpe']:.2f}, " f"WR {base_wr:.1f}%→{best['wr_pct']:.1f}%, " f"loses {markov_results['1_baseline']['n']-best['n']} trades/yr") elif best["sharpe"] > base_sh: rec = f"MARGINAL: {best['label']} — small improvement, not worth complexity" else: rec = "NO IMPROVEMENT: Keep model as-is, Markov adds no edge" print(f"\n Baseline : {base_wr:.1f}% WR, Sharpe {base_sh:.2f}") print(f" Best filter: {rec}") print() # ── Brownian Motion / Hurst backtest ───────────────────────────────────── print("\n" + "="*70) print("BROWNIAN MOTION / HURST EXPONENT REGIME FILTER") print("="*70) bm_results, bm_df = run_brownian_backtest(trades, z_series) print_results(bm_results, "HURST EXPONENT RESULTS") # ── All approaches combined ─────────────────────────────────────────────── all_r = {**markov_results} if hmm_results: all_r.update(hmm_results) all_r.update(bm_results) print_results(all_r, "ALL APPROACHES — FINAL COMPARISON") print("\nINTERPRETATION: Goal = Sharpe +0.3 above baseline with >80% of trades") print("If no approach meets this threshold, keep model exactly as validated.\n") # Save out = PROC_PATH / "markov_regime_results.csv" trade_df.to_csv(out, index=False) print(f"Trade-level regime classifications saved: {out.name}") if __name__ == "__main__": main()