research: 3-year stake-sizing campaign + violent-move studies (nothing changes signals)

Walk-forward "second opinion" stake multiplier (meta-labelling on 65 decision-time
features, 3yr Oanda books @ $0.65) FAILED the pre-registered rule on M5 (clean null)
and M15 (lift inside the label-shuffle null, explained by BUY-vs-SELL). Nothing that
changes numbers ships; the 12-month M5 book re-run is byte-identical to the pinned one.

Behaviour-neutral plumbing:
- position_sizing.build_features (65 features; multiplier path returns 1.0, no artifact)
- SignalReport.penalty_tags / penalty_total / momentum_fading (record-keeping)
- TradeSetup.stake_multiplier = None (reserved)
- gpu_backtest --export-features (default export stays the pinned 63 columns; extras
  forced to the tail even when the first trade lacks an indicator)

Research tools (backtest/research): book_report, burst_map, burst_precursors,
open_trade_bursts, calendar_impact, burst_indicator_map, sizing_model, watch_rule_eval.
Findings recorded in backtest/research/sizing_campaign_2026-09.md: M5 stretch measures
are positive in both halves; every open-trade protection against violent bars loses on
both TFs; high-impact USD releases cause 20-31% of violent bars but trades held through
them do better; violent bars come when price is already at an extreme, not from calm.

Module 29 tests (13). Suite: 435 passed.

Claude-Session: https://claude.ai/code/session_01RWV2yLt3YXHHGVYHSQU2EE
This commit is contained in:
Harwinder Sekhon
2026-09-05 21:19:24 -04:00
parent bed8b549ac
commit 6bf58fd9d4
16 changed files with 4076 additions and 3 deletions
+581
View File
@@ -0,0 +1,581 @@
"""Full descriptive report of a backtest trade book (research only).
Usage:
python backtest/research/book_report.py --book BOOK.csv --bar-minutes N \
[--split YYYY-MM-DD|median] [--parquet-4h PATH] [--parquet-1d PATH] \
[--out DIR]
Prints (and writes to DIR/report.txt) headline metrics for the whole book
and for each half, group tables (weekday / decision hour / session / quality
/ regime / vol_regime / direction / strategy / exit_reason / indicators_agree),
win-loss streaks, a rolling win-rate summary, and a quantile-bucket screen of
every numeric decision-time feature the book carries (`f_*` / `x_*` columns
plus score, confidence, rr_planned, indicators_agree/oppose and `ind_*_value`).
The screen is deliberately half-split: a feature only "passes" when the
top-bucket-minus-bottom-bucket mean-R gap has the SAME SIGN in both halves of
the book and both extreme buckets are populated in both halves. Bin edges are
fit on the WHOLE book (so the halves share edges and stay comparable) — that
leaks a little cross-half information into the edges themselves, which is why
the sign test, not the size of the gap, is the thing to trust.
Everything here is measurement, not simulation: it never re-runs the signal
engine, it only describes a book that gpu_backtest.py already produced.
Analysis tooling only — never imported by the live code.
"""
from __future__ import annotations
import argparse
import math
import os
from typing import Iterable
import numpy as np
import pandas as pd
SESSIONS = (
("asian", 3, 6),
("london", 7, 12),
("overlap", 13, 16),
("ny", 17, 22),
)
# Numeric non-`f_`/`x_` columns that are known at DECISION time and are
# therefore legal screen candidates. bars_held / r_multiple / mfe_r / mae_r
# are post-entry outcomes and must never appear here.
BASE_SCREEN_COLS = ("score", "confidence", "rr_planned",
"indicators_agree", "indicators_oppose")
def session_of(hour: int) -> str:
for name, lo, hi in SESSIONS:
if lo <= hour <= hi:
return name
return "off"
def feature_columns(df: pd.DataFrame) -> list[str]:
"""Decision-time feature columns carried by this book, generically."""
extra = [c for c in df.columns if c.startswith(("f_", "x_"))]
base = [c for c in BASE_SCREEN_COLS if c in df.columns]
ind = [c for c in df.columns
if c.startswith("ind_") and c.endswith("_value")]
out, seen = [], set()
for c in extra + base + ind:
if c not in seen:
seen.add(c)
out.append(c)
return out
# ───────────────────────────── loading ──────────────────────────────
def load_book(csv: str, bar_minutes: int = 5,
split: str = "median") -> pd.DataFrame:
"""Read a gpu_backtest trade book and add the decision-time columns.
`decision_time` = entry_time − one bar: entries fill at the NEXT bar's
open, so the bar the signal actually decided on closed one bar earlier.
That is the clock blocked_hours_utc and every session filter evaluate.
`half` is 1 before the split instant and 2 at/after it.
"""
df = pd.read_csv(csv)
for col in ("entry_time", "exit_time"):
if col in df.columns:
df[col] = pd.to_datetime(df[col], utc=True, format="mixed")
df["decision_time"] = df["entry_time"] - pd.Timedelta(minutes=bar_minutes)
df["decision_hour"] = df["decision_time"].dt.hour
df["weekday"] = df["decision_time"].dt.weekday
df["session"] = df["decision_hour"].map(session_of)
df = df.sort_values("entry_time", kind="mergesort").reset_index(drop=True)
df["half"] = _halves(df["entry_time"], split)
return df
def split_instant(entry_time: pd.Series, split: str) -> pd.Timestamp:
if split == "median":
return entry_time.median()
return pd.Timestamp(split, tz="UTC")
def _halves(entry_time: pd.Series, split: str) -> pd.Series:
cut = split_instant(entry_time, split)
return np.where(entry_time < cut, 1, 2)
# ───────────────────────────── metrics ──────────────────────────────
def base_metrics(r: pd.Series, df: pd.DataFrame | None = None) -> dict:
"""Headline stats for a series of R multiples.
max_dd uses the gpu_backtest.compute_metrics convention: a 0.0 is
prepended to the cumulative-R curve so an opening losing streak counts
as drawdown from flat, not from an already-negative first trade.
pf/expectancy split on r > 0 / r < 0 (exact zeros are neither).
"""
r = pd.Series(r, dtype="float64").dropna()
n = int(len(r))
if n == 0:
return {"n": 0, "total_r": 0.0, "mean_r": 0.0, "median_r": 0.0,
"wr": 0.0, "pf": 0.0, "max_dd": 0.0, "sharpe": 0.0,
"expectancy": 0.0, "expectancy_alt": 0.0, "recovery": 0.0,
"avg_bars": float("nan"), "avg_mfe": float("nan"),
"avg_mae": float("nan")}
wins, losses = r[r > 0], r[r < 0]
gross_win = float(wins.sum())
gross_loss = float(-losses.sum())
cum = np.concatenate(([0.0], np.cumsum(r.to_numpy())))
max_dd = float((np.maximum.accumulate(cum) - cum).max())
std = float(r.std(ddof=0))
mean_r = float(r.mean())
wr = len(wins) / n
lr = len(losses) / n
avg_win = float(wins.mean()) if len(wins) else 0.0
avg_loss = float(-losses.mean()) if len(losses) else 0.0
expectancy_alt = wr * avg_win - lr * avg_loss
# Identity check: the textbook expectancy form must equal the plain mean
# whenever no trade is exactly flat (and it does even then, since a 0 R
# trade contributes 0 to both sides).
assert math.isclose(expectancy_alt, mean_r, rel_tol=1e-9, abs_tol=1e-9), (
f"expectancy mismatch: {expectancy_alt} vs mean {mean_r}")
out = {
"n": n,
"total_r": float(r.sum()),
"mean_r": mean_r,
"median_r": float(r.median()),
"wr": wr * 100.0,
"pf": (gross_win / gross_loss) if gross_loss > 0 else float("inf"),
"max_dd": max_dd,
"sharpe": (mean_r / std) if std > 0 else 0.0,
"expectancy": mean_r,
"expectancy_alt": expectancy_alt,
"recovery": (float(r.sum()) / max_dd) if max_dd > 0 else float("inf"),
"avg_bars": float("nan"),
"avg_mfe": float("nan"),
"avg_mae": float("nan"),
}
if df is not None:
for key, col in (("avg_bars", "bars_held"), ("avg_mfe", "mfe_r"),
("avg_mae", "mae_r")):
if col in df.columns:
out[key] = float(pd.to_numeric(df[col],
errors="coerce").mean())
return out
def _pf(r: pd.Series) -> float:
"""Profit factor, hour_analysis-compatible (loss side is r <= 0)."""
pos = r[r > 0].sum()
neg = -r[r <= 0].sum()
return float(pos / max(neg, 1e-9))
def group_table(df: pd.DataFrame, by, min_n: int = 10) -> pd.DataFrame:
"""Per-group P&L with a half-1 / half-2 sign-stability flag."""
keys = [by] if isinstance(by, str) else list(by)
rows = []
for key, sub in df.groupby(keys, dropna=False, sort=True):
r = sub["r_multiple"]
h1, h2 = sub[sub["half"] == 1]["r_multiple"], sub[sub["half"] == 2]["r_multiple"]
m1 = float(h1.mean()) if len(h1) else float("nan")
m2 = float(h2.mean()) if len(h2) else float("nan")
stable = bool(len(h1) >= min_n and len(h2) >= min_n
and np.sign(m1) != 0 and np.sign(m1) == np.sign(m2))
rec = {}
# pandas >= 2.2 hands back a 1-tuple even for a single group key.
vals = key if isinstance(key, tuple) else (key,)
for k, v in zip(keys, vals):
rec[k] = v
rec.update({
"n": int(len(sub)),
"total_r": round(float(r.sum()), 2),
"mean_r": round(float(r.mean()), 3),
"wr": round(float((r > 0).mean() * 100), 2),
"pf": round(_pf(r), 2),
"n_h1": int(len(h1)),
"total_r_h1": round(float(h1.sum()), 2),
"mean_r_h1": round(m1, 3) if len(h1) else float("nan"),
"n_h2": int(len(h2)),
"total_r_h2": round(float(h2.sum()), 2),
"mean_r_h2": round(m2, 3) if len(h2) else float("nan"),
"sign_stable": stable,
})
rows.append(rec)
out = pd.DataFrame(rows)
if not out.empty:
out = out.set_index(keys)
return out
def streaks(df: pd.DataFrame) -> dict:
"""Longest win/loss runs and the longest stretch under water (in trades)."""
d = df.sort_values("entry_time", kind="mergesort")
r = d["r_multiple"].to_numpy(dtype="float64")
best_w = best_l = cur_w = cur_l = 0
for x in r:
if x > 0:
cur_w, cur_l = cur_w + 1, 0
elif x < 0:
cur_l, cur_w = cur_l + 1, 0
else:
cur_w = cur_l = 0
best_w, best_l = max(best_w, cur_w), max(best_l, cur_l)
cum = np.concatenate(([0.0], np.cumsum(r)))
peak = np.maximum.accumulate(cum)
flat, longest = 0, 0
for i in range(1, len(cum)):
flat = 0 if cum[i] >= peak[i] else flat + 1
longest = max(longest, flat)
return {"max_win_streak": best_w, "max_loss_streak": best_l,
"longest_flat_trades": longest}
def rolling_wr(df: pd.DataFrame, window: int = 20) -> pd.Series:
d = df.sort_values("entry_time", kind="mergesort")
s = (d["r_multiple"] > 0).astype(float).rolling(window).mean() * 100
s.index = d["entry_time"].to_numpy()
return s
def equity_curve(df: pd.DataFrame) -> pd.DataFrame:
d = df.sort_values("entry_time", kind="mergesort")
cum = d["r_multiple"].cumsum()
peak = pd.concat([pd.Series([0.0]), cum]).cummax().iloc[1:]
peak.index = cum.index
return pd.DataFrame({"entry_time": d["entry_time"].to_numpy(),
"cum_r": cum.to_numpy(),
"drawdown": (peak - cum).to_numpy()})
# ─────────────────────── quantile bucket screen ─────────────────────
def quantile_buckets(df: pd.DataFrame, col: str, q: int = 5,
min_n: int = 10) -> pd.DataFrame:
"""Equal-frequency buckets of `col` with per-half mean R.
Bin edges are fit on the WHOLE book (both halves share edges).
Result carries gap_h1 / gap_h2 / sign_stable in `.attrs` and as
repeated columns so the CSV is self-describing.
"""
v = pd.to_numeric(df.get(col), errors="coerce")
d = df.loc[v.notna()].copy()
d["_v"] = v.loc[v.notna()]
empty = pd.DataFrame(columns=["bin", "lo", "hi", "n", "mean_r", "pf", "wr",
"n_h1", "mean_r_h1", "n_h2", "mean_r_h2",
"gap_h1", "gap_h2", "sign_stable"])
empty.attrs = {"feature": col, "gap_h1": float("nan"),
"gap_h2": float("nan"), "sign_stable": False, "n": 0,
"mode": "empty"}
if d.empty:
return empty
nuniq = int(d["_v"].nunique())
if nuniq < 3:
mode = "values"
d["_bin"] = d["_v"]
order = sorted(d["_v"].unique())
edges = {b: (float(b), float(b)) for b in order}
else:
mode = "qcut"
binned = pd.qcut(d["_v"], q, duplicates="drop")
d["_bin"] = binned
order = list(binned.cat.categories)
edges = {b: (float(b.left), float(b.right)) for b in order}
rows = []
for b in order:
sub = d[d["_bin"] == b]
if sub.empty:
continue
r = sub["r_multiple"]
h1 = sub[sub["half"] == 1]["r_multiple"]
h2 = sub[sub["half"] == 2]["r_multiple"]
lo, hi = edges[b]
rows.append({
"bin": str(b), "lo": round(lo, 4), "hi": round(hi, 4),
"n": int(len(sub)), "mean_r": round(float(r.mean()), 3),
"pf": round(_pf(r), 2),
"wr": round(float((r > 0).mean() * 100), 2),
"n_h1": int(len(h1)),
"mean_r_h1": round(float(h1.mean()), 3) if len(h1) else float("nan"),
"n_h2": int(len(h2)),
"mean_r_h2": round(float(h2.mean()), 3) if len(h2) else float("nan"),
})
out = pd.DataFrame(rows)
if out.empty:
return empty
bot, top = out.iloc[0], out.iloc[-1]
gap_h1 = float(top["mean_r_h1"] - bot["mean_r_h1"])
gap_h2 = float(top["mean_r_h2"] - bot["mean_r_h2"])
populated = min(bot["n_h1"], bot["n_h2"], top["n_h1"], top["n_h2"]) >= min_n
signs_ok = (not math.isnan(gap_h1) and not math.isnan(gap_h2)
and np.sign(gap_h1) != 0
and np.sign(gap_h1) == np.sign(gap_h2))
sign_stable = bool(populated and signs_ok and len(out) >= 2)
out["gap_h1"] = round(gap_h1, 3) if not math.isnan(gap_h1) else float("nan")
out["gap_h2"] = round(gap_h2, 3) if not math.isnan(gap_h2) else float("nan")
out["sign_stable"] = sign_stable
out.attrs = {"feature": col, "gap_h1": gap_h1, "gap_h2": gap_h2,
"sign_stable": sign_stable, "n": int(len(d)), "mode": mode}
return out
def feature_screen(df: pd.DataFrame, q: int = 5, min_n: int = 10,
columns: Iterable[str] | None = None) -> pd.DataFrame:
"""Bucket screen over every decision-time numeric feature in the book."""
cols = list(columns) if columns is not None else feature_columns(df)
rows = []
for c in cols:
bk = quantile_buckets(df, c, q=q, min_n=min_n)
a = bk.attrs
stable = bool(a.get("sign_stable"))
g1, g2 = a.get("gap_h1", float("nan")), a.get("gap_h2", float("nan"))
min_abs = (min(abs(g1), abs(g2)) if stable else 0.0)
rows.append({"feature": c,
"gap_h1": round(g1, 3) if g1 == g1 else float("nan"),
"gap_h2": round(g2, 3) if g2 == g2 else float("nan"),
"min_abs_gap": round(float(min_abs), 3),
"sign_stable": stable,
"n": int(a.get("n", 0)),
"bins": int(len(bk))})
out = pd.DataFrame(rows)
if out.empty:
return out
return out.sort_values(["min_abs_gap", "n"],
ascending=[False, False]).reset_index(drop=True)
# ────────────────────── higher-timeframe extension ──────────────────
def wilder_atr(bars: pd.DataFrame, period: int = 14) -> pd.Series:
h, l, c = bars["High"], bars["Low"], bars["Close"]
tr = pd.concat([h - l, (h - c.shift()).abs(), (l - c.shift()).abs()],
axis=1).max(axis=1)
return tr.ewm(alpha=1 / period, min_periods=period, adjust=False).mean()
def _ext_series(parquet: str, minutes: int) -> pd.DataFrame:
bars = pd.read_parquet(parquet).sort_index()
bars = bars[~bars.index.duplicated(keep="last")]
ext = ((bars["Close"] - bars["Close"].rolling(20).mean())
/ wilder_atr(bars, 14))
out = pd.DataFrame({
"close_time": bars.index + pd.Timedelta(minutes=minutes),
"ext": ext.to_numpy(),
}).dropna().reset_index(drop=True)
# Books parse to datetime64[us, UTC] (format="mixed"), parquet indexes are
# [ns, UTC]; merge_asof refuses mismatched resolutions.
out["close_time"] = out["close_time"].astype("datetime64[ns, UTC]")
return out
def add_htf_extension(df: pd.DataFrame, parquet_4h: str | None,
parquet_1d: str | None) -> pd.DataFrame:
"""Add direction-signed H4 / D1 stretch-from-mean in ATR units.
ext = (close − 20-bar SMA of close) / Wilder ATR(14) on that timeframe,
taken from the LAST higher-TF bar whose CLOSE time (bar label + one bar)
is <= decision_time — i.e. only information the live scan could have had.
Signed +1 for BUY, −1 for SELL, so positive always means "extended in the
direction of the trade".
"""
out = df.copy()
sign = np.where(out["direction"].astype(str).str.upper().str.startswith("B"),
1.0, -1.0)
for col, parquet, minutes in (("x_h4_ext_atr", parquet_4h, 240),
("x_d1_ext_atr", parquet_1d, 1440)):
if not parquet:
continue
right = _ext_series(parquet, minutes)
left = out[["decision_time"]].reset_index()
left["decision_time"] = left["decision_time"].astype(
"datetime64[ns, UTC]")
left = left.sort_values("decision_time", kind="mergesort")
merged = pd.merge_asof(left, right, left_on="decision_time",
right_on="close_time", direction="backward")
merged["_stale_h"] = ((merged["decision_time"] - merged["close_time"])
.dt.total_seconds() / 3600.0)
merged = merged.set_index("index").reindex(out.index)
out[col] = merged["ext"].to_numpy() * sign
out.attrs[f"{col}_max_stale_h"] = float(
merged["_stale_h"].max()) if len(merged) else float("nan")
return out
# ──────────────────────────── reporting ─────────────────────────────
def _fmt(v, nd=3):
if isinstance(v, float):
if math.isinf(v):
return "inf"
if math.isnan(v):
return "n/a"
return f"{v:.{nd}f}"
return str(v)
def _metrics_block(name: str, m: dict) -> list[str]:
return [
f" {name:<10} n={m['n']:<5} totalR={_fmt(m['total_r'], 2):>8} "
f"meanR={_fmt(m['mean_r'])} medR={_fmt(m['median_r'])} "
f"WR={_fmt(m['wr'], 2)}% PF={_fmt(m['pf'], 2)}",
f" {'':<10} maxDD={_fmt(m['max_dd'], 2)} recovery={_fmt(m['recovery'], 2)}"
f" sharpe={_fmt(m['sharpe'])} expectancy={_fmt(m['expectancy'])}"
f" (wr·avgW − lr·avgL = {_fmt(m['expectancy_alt'])})",
f" {'':<10} avg bars={_fmt(m['avg_bars'], 1)} "
f"avg MFE={_fmt(m['avg_mfe'])}R avg MAE={_fmt(m['avg_mae'])}R",
]
def build_report(df: pd.DataFrame, book: str, bar_minutes: int, split: str,
out_dir: str, q: int = 5, min_n: int = 10) -> str:
lines: list[str] = []
add = lines.append
cut = split_instant(df["entry_time"], split)
add("=" * 78)
add(f"BOOK REPORT {os.path.basename(book)}")
add("=" * 78)
add(f"trades : {len(df)}")
add(f"entry window : {df['entry_time'].min()} -> {df['entry_time'].max()}")
add(f"bar minutes : {bar_minutes} (decision_time = entry_time − 1 bar)")
add(f"split rule : {split} -> cut at {cut} "
f"(half 1 = before, half 2 = at/after)")
add(f"half sizes : h1={int((df['half'] == 1).sum())} "
f"h2={int((df['half'] == 2).sum())}")
add("session map : asian 3-6, london 7-12, overlap 13-16, ny 17-22, "
"else off (on DECISION hour, UTC)")
add("bucket edges : quantile bins are fit on the WHOLE book; the two "
"halves share the same edges.")
feats = feature_columns(df)
add(f"screen features : {len(feats)} "
f"({sum(1 for c in feats if c.startswith(('f_', 'x_')))} f_/x_ columns)")
for k, v in df.attrs.items():
add(f"htf staleness : {k} = {_fmt(v, 1)} h (max age of the matched bar)")
add("")
add("── HEADLINE METRICS ──")
lines.extend(_metrics_block("full", base_metrics(df["r_multiple"], df)))
for h in (1, 2):
sub = df[df["half"] == h]
lines.extend(_metrics_block(f"half {h}",
base_metrics(sub["r_multiple"], sub)))
group_specs = [
("weekday", "weekday"), ("decision_hour", "decision_hour"),
("session", "session"), ("quality", "quality"), ("regime", "regime"),
("vol_regime", "vol_regime"), ("direction", "direction"),
("strategy", "strategy"), ("exit_reason", "exit_reason"),
("indicators_agree", "indicators_agree"),
]
tables = {}
for name, col in group_specs:
if col not in df.columns:
continue
t = group_table(df, col, min_n=min_n)
tables[name] = t
add("")
add(f"── BY {name.upper()} ──")
add(t.to_string() if not t.empty else " (empty)")
# (weekday, decision_hour) pivots
add("")
add("── weekday x decision_hour : TOTAL R ──")
wh = group_table(df, ["weekday", "decision_hour"], min_n=min_n)
if not wh.empty:
piv_r = wh["total_r"].unstack("decision_hour")
piv_n = wh["n"].unstack("decision_hour")
add(piv_r.fillna(0).round(2).to_string())
add("")
add("── weekday x decision_hour : N ──")
add(piv_n.fillna(0).astype(int).to_string())
tables["weekday_hour"] = wh
piv_r.to_csv(os.path.join(out_dir, "groups_weekday_hour_total_r.csv"))
piv_n.to_csv(os.path.join(out_dir, "groups_weekday_hour_n.csv"))
st = streaks(df)
add("")
add("── STREAKS ──")
add(f" max consecutive wins : {st['max_win_streak']}")
add(f" max consecutive losses : {st['max_loss_streak']}")
add(f" longest flat stretch : {st['longest_flat_trades']} trades "
f"between equity highs")
rw = rolling_wr(df, 20).dropna()
add("")
add("── ROLLING WIN RATE (20 trades) ──")
if len(rw):
add(f" min={_fmt(rw.min(), 2)}% median={_fmt(rw.median(), 2)}% "
f"max={_fmt(rw.max(), 2)}% n={len(rw)}")
add(" last 5: " + ", ".join(
f"{ts:%Y-%m-%d %H:%M}={v:.2f}%" for ts, v in rw.tail(5).items()))
else:
add(" (fewer than 20 trades)")
scr = feature_screen(df, q=q, min_n=min_n)
add("")
add("── FEATURE SCREEN (top 25 by min |gap| across halves) ──")
add(" gap = mean R(top bucket) − mean R(bottom bucket), per half; "
"min_abs_gap is 0 unless the sign is stable.")
add(scr.head(25).to_string(index=False) if not scr.empty else " (none)")
stable = scr[scr["sign_stable"]].head(10) if not scr.empty else scr
add("")
add(f"── BUCKET TABLES FOR THE TOP {len(stable)} SIGN-STABLE FEATURES ──")
if stable.empty:
add(" (no feature has a sign-stable top-vs-bottom gap)")
for feat in stable["feature"]:
bk = quantile_buckets(df, feat, q=q, min_n=min_n)
add("")
add(f" · {feat} (mode={bk.attrs['mode']}, n={bk.attrs['n']}, "
f"gap_h1={_fmt(bk.attrs['gap_h1'])}, "
f"gap_h2={_fmt(bk.attrs['gap_h2'])})")
add(bk.to_string(index=False))
bk.to_csv(os.path.join(out_dir, f"buckets_{feat}.csv"), index=False)
# CSV artefacts
for name, t in tables.items():
t.to_csv(os.path.join(out_dir, f"groups_{name}.csv"))
scr.to_csv(os.path.join(out_dir, "feature_screen.csv"), index=False)
equity_curve(df).to_csv(os.path.join(out_dir, "equity.csv"), index=False)
all_buckets = []
for c in feats:
bk = quantile_buckets(df, c, q=q, min_n=min_n)
if bk.empty:
continue
bk = bk.copy()
bk.insert(0, "feature", c)
all_buckets.append(bk)
if all_buckets:
pd.concat(all_buckets, ignore_index=True).to_csv(
os.path.join(out_dir, "buckets.csv"), index=False)
return "\n".join(lines)
def main() -> None:
ap = argparse.ArgumentParser(description=__doc__.split("\n")[0])
ap.add_argument("--book", required=True)
ap.add_argument("--bar-minutes", type=int, required=True)
ap.add_argument("--split", default="median",
help="'median' (default) or YYYY-MM-DD")
ap.add_argument("--parquet-4h", default=None)
ap.add_argument("--parquet-1d", default=None)
ap.add_argument("--q", type=int, default=5)
ap.add_argument("--min-n", type=int, default=10)
ap.add_argument("--out", default="backtest/research/book_report_out")
args = ap.parse_args()
os.makedirs(args.out, exist_ok=True)
df = load_book(args.book, args.bar_minutes, args.split)
if args.parquet_4h or args.parquet_1d:
df = add_htf_extension(df, args.parquet_4h, args.parquet_1d)
text = build_report(df, args.book, args.bar_minutes, args.split,
args.out, q=args.q, min_n=args.min_n)
print(text)
with open(os.path.join(args.out, "report.txt"), "w") as fh:
fh.write(text + "\n")
print(f"\n[written] {os.path.join(args.out, 'report.txt')}")
if __name__ == "__main__":
main()
+121
View File
@@ -0,0 +1,121 @@
"""What do the system's OWN indicators read right before a violent bar? Analysis tooling only.
Uses backtest.gpu_backtest.precompute_indicators (the same series the backtester feeds the
engine) on 3 years of bars; every state is taken at t-1 (known before the bar). Reports burst
rate (range >= 3x prior ATR14) for the NEXT bar and within the NEXT HOUR by indicator state,
by regime, by volatility regime, by weekday x hour, day of month, calendar proximity — each in
both halves — and whether the state predicts the burst's DIRECTION.
python backtest/research/burst_indicator_map.py PARQUET BAR_MINUTES [--start --end]
"""
import argparse, sys
from pathlib import Path
import numpy as np, pandas as pd
ROOT = Path(__file__).resolve().parents[2]; sys.path.insert(0, str(ROOT)); sys.path.insert(0, str(ROOT / "backtest"))
from gpu_backtest import precompute_indicators
from config import TechnicalConfig
CAL = ROOT / "backtest" / "cache" / "calendar.parquet"
def wilder_atr(df, period=14):
pc = df["Close"].shift(1)
tr = pd.concat([df["High"] - df["Low"], (df["High"] - pc).abs(), (df["Low"] - pc).abs()], axis=1).max(axis=1)
return tr.ewm(alpha=1 / period, min_periods=period, adjust=False).mean()
def rate_table(d, col, target, half, order=None):
base1 = d.loc[half == 1, target].mean(); base2 = d.loc[half == 2, target].mean()
g = d.groupby(col, observed=True)
out = pd.DataFrame({"n": g.size(),
"h1_lift": d[half == 1].groupby(col, observed=True)[target].mean() / base1,
"h2_lift": d[half == 2].groupby(col, observed=True)[target].mean() / base2,
"rate%": 100 * g[target].mean()})
if order is not None: out = out.reindex([o for o in order if o in out.index])
out["stable"] = np.where((out.h1_lift > 1.15) & (out.h2_lift > 1.15), "HIGH", np.where((out.h1_lift < 0.85) & (out.h2_lift < 0.85), "low", ""))
return out.round(2)
def main():
ap = argparse.ArgumentParser(); ap.add_argument("parquet"); ap.add_argument("bar_minutes", type=int)
ap.add_argument("--start", default="2023-09-05"); ap.add_argument("--end", default="2026-09-05")
a = ap.parse_args(); bm = a.bar_minutes
px = pd.read_parquet(a.parquet)
px = px[(px.index >= pd.Timestamp(a.start, tz="UTC") - pd.Timedelta(days=30)) & (px.index < pd.Timestamp(a.end, tz="UTC"))]
ic = precompute_indicators(px, TechnicalConfig())
d = pd.DataFrame(index=px.index)
atr_prev = wilder_atr(px).shift(1)
rng = px["High"] - px["Low"]
d["burst"] = (rng >= 3.0 * atr_prev).astype(int)
H = max(1, 60 // bm)
d["next_hour"] = d["burst"][::-1].rolling(H, min_periods=1).max()[::-1].astype(int)
d["next_bar"] = d["burst"]
d["bdir"] = np.sign(px["Close"] - px["Open"])
# states at t-1
sh = lambda arr: pd.Series(np.asarray(arr, dtype=float), index=px.index).shift(1)
rsi = sh(ic.rsi); adx = sh(ic.adx); adxs = sh(ic.adx_slope); hist = sh(ic.macd_hist); bbp = sh(ic.bb_position)
bbw = sh(ic.bb_band_width); bbaw = sh(ic.bb_avg_width); kcsq = sh(ic.kc_squeeze); roc = sh(ic.roc)
ef, em, es = sh(ic.ema_fast), sh(ic.ema_medium), sh(ic.ema_slow); et = sh(ic.ema_trend) if ic.ema_trend is not None else None
close1 = px["Close"].shift(1); atr1 = sh(ic.atr)
kcu, kcl = sh(ic.kc_upper), sh(ic.kc_lower)
d["rsi"] = pd.cut(rsi, [0, 30, 40, 50, 60, 70, 100], labels=["<30", "30-40", "40-50", "50-60", "60-70", ">=70"])
d["adx"] = pd.cut(adx, [0, 20, 25, 35, 50, 200], labels=["<20", "20-25", "25-35", "35-50", ">=50"])
d["adx_slope"] = pd.cut(adxs, [-1e9, -3, 0, 3, 1e9], labels=["falling<-3", "-3..0", "0..3", "rising>3"])
hist_prev = hist.shift(1)
d["macd"] = np.select([(hist > 0) & (hist > hist_prev), (hist > 0), (hist < 0) & (hist < hist_prev), (hist < 0)],
["pos&rising", "pos&fading", "neg&falling", "neg&fading"], default="flat")
d["bb_pos"] = pd.cut(bbp, [-1e9, 0, 0.2, 0.8, 1.0, 1e9], labels=["<0 (below)", "0-0.2", "0.2-0.8", "0.8-1", ">1 (above)"])
d["bb_squeeze"] = np.where(bbw < 0.75 * bbaw, "squeeze", "normal")
d["kc_squeeze"] = np.where(kcsq > 0, "TTM squeeze", "no")
d["kc_pos"] = np.select([close1 > kcu, close1 < kcl], ["above KC", "below KC"], default="inside")
align = np.sign(close1 - ef) + np.sign(ef - em) + np.sign(em - es) + (np.sign(es - et) if et is not None else 0)
d["ema_align"] = pd.cut(align, [-5, -3.5, -0.5, 0.5, 3.5, 5], labels=["fully bear(-4)", "mostly bear", "mixed", "mostly bull", "fully bull(+4)"])
ribbon = (pd.concat([ef, em, es], axis=1).max(axis=1) - pd.concat([ef, em, es], axis=1).min(axis=1)) / atr1
d["ribbon"] = np.where(ribbon < 1.0, "compressed(<1 ATR)", "expanded")
ema_dir = np.where(align > 0, 1, np.where(align < 0, -1, 0))
d["regime"] = np.select([(adx >= 25) & (ema_dir > 0), (adx >= 25) & (ema_dir < 0), adx >= 25, adxs > 3.0],
["trending_up", "trending_down", "trending(mixed)", "transitioning"], default="ranging")
atr_pct = atr1.rolling(50).apply(lambda w: (w[-1] > w).mean() * 100, raw=True)
d["vol_regime"] = pd.cut(atr_pct, [-1, 30, 60, 80, 101], labels=["low", "normal", "high", "extreme"])
d["roc"] = pd.cut(roc, [-1e9, -1, -0.3, 0.3, 1, 1e9], labels=["<-1%", "-1..-0.3", "-0.3..0.3", "0.3..1", ">1%"])
d["ema200_dist"] = pd.cut((close1 - et) / et * 100 if et is not None else pd.Series(np.nan, index=px.index),
[-1e9, -2.5, -0.5, 0.5, 2.5, 1e9], labels=["<-2.5%", "-2.5..-0.5", "near(±0.5%)", "0.5..2.5", ">2.5%"])
# time
ct = px.index + pd.Timedelta(minutes=bm)
d["hour"] = ct.hour; d["weekday"] = ct.weekday; d["dom"] = pd.cut(ct.day, [0, 5, 10, 15, 20, 25, 31], labels=["1-5", "6-10", "11-15", "16-20", "21-25", "26-31"])
d["first_friday"] = np.where((ct.weekday == 4) & (ct.day <= 7), "first Friday", "other")
# calendar proximity (high-impact USD)
cal = pd.read_parquet(CAL); cal["ts"] = pd.to_datetime(cal.dateline, unit="s", utc=True)
hi = np.sort(cal[(cal.currency == "USD") & (cal.impactName.str.lower() == "high")].ts.values.astype("datetime64[ns]"))
t = ct.values.astype("datetime64[ns]")
nxt = np.searchsorted(hi, t) # next event at/after bar close
mins_to = np.where(nxt < len(hi), (hi[np.minimum(nxt, len(hi) - 1)] - t) / np.timedelta64(1, "m"), 1e9)
prv = nxt - 1
mins_since = np.where(prv >= 0, (t - hi[np.maximum(prv, 0)]) / np.timedelta64(1, "m"), 1e9)
d["calendar"] = np.select([mins_since <= 30, mins_to <= 15, mins_to <= 60, mins_to <= 180], ["0-30m AFTER hi-impact", "next 15m", "15-60m ahead", "1-3h ahead"], default="none within 3h")
d = d[d.index >= pd.Timestamp(a.start, tz="UTC")].dropna(subset=["rsi", "adx", "bb_pos", "vol_regime"])
half = pd.Series(np.where(np.arange(len(d)) < len(d) // 2, 1, 2), index=d.index)
print(f"bars {len(d)} base next_bar {100*d.next_bar.mean():.3f}% next_hour {100*d.next_hour.mean():.2f}% (lift = rate / base; 'HIGH' = >1.15x in BOTH halves, 'low' = <0.85x in both)")
groups = ["calendar", "regime", "vol_regime", "rsi", "adx", "adx_slope", "macd", "bb_pos", "bb_squeeze", "kc_squeeze", "kc_pos", "ema_align", "ribbon", "roc", "ema200_dist", "dom", "first_friday"]
for target in ("next_bar", "next_hour"):
print(f"\n################ TARGET: {target} ################")
for g in groups:
print(f"-- {g}"); print(rate_table(d, g, target, half).to_string())
# weekday x hour top cells (next_bar)
piv = d.pivot_table(index="weekday", columns="hour", values="next_bar", aggfunc="mean") * 100
base = 100 * d.next_bar.mean()
cells = piv.stack().sort_values(ascending=False)
print(f"\n################ TOP 15 weekday x close-hour cells (next-bar burst pct, base {base:.2f}) ################")
names = {0: "Mon", 1: "Tue", 2: "Wed", 3: "Thu", 4: "Fri", 6: "Sun"}
for (wd, h), v in cells.head(15).items():
n = int(((d.weekday == wd) & (d.hour == h)).sum())
print(f" {names.get(wd, wd)} {h:02d}:xx {v:.2f}% ({v/base:.1f}x, n={n})")
# direction: does the state predict the burst direction?
b = d[d.burst == 1]
print("\n################ BURST DIRECTION vs prior state (share of bursts going UP) ################")
for g, up_states in (("ema_align", None), ("rsi", None), ("macd", None), ("regime", None)):
print(f"-- {g}: " + ", ".join(f"{k}: {100*(v>0).mean():.0f}% up (n={len(v)})" for k, v in b.groupby(g, observed=True).bdir if len(v) >= 20))
if __name__ == "__main__":
main()
+445
View File
@@ -0,0 +1,445 @@
"""Where and when gold moves violently — burst map of the raw bars.
Usage:
python backtest/research/burst_map.py --parquet PATH --bar-minutes N \
[--k 3,4,5] [--book BOOK.csv] [--within-bars 12] \
[--start YYYY-MM-DD --end YYYY-MM-DD] [--out DIR]
A "burst" is one bar whose range is >= K x the PRIOR bar's Wilder ATR(14) —
the same definition burst_profile.py uses (that module computes it inline
inside burst_table() and exposes no reusable helper, so the formula is
replicated here bar-for-bar: TR = max(H-L, |H-Cprev|, |L-Cprev|), ATR =
TR.ewm(alpha=1/14, min_periods=14, adjust=False).mean(), ratio =
(H-L) / ATR.shift(1)).
Reports:
* burst frequency by (weekday, UTC bar-CLOSE hour), plus mean |1 h| and
|4 h| absolute moves in $ and in ATR units on the same grid;
* post-burst drift in the burst's own direction over 6/12/24/48 bars,
overall and by session bucket — does a violent bar continue or snap back;
* the largest trailing-60-minute close-to-close moves, de-overlapped;
* if a trade book is given: how trades that were decided just after a burst,
or that were OPEN when one printed, performed against the rest.
Every ATR-unit number here is normalised by the PRIOR bar's ATR — the same
denominator as the burst test, so a "4x burst that ran another 4 ATR" reads
directly. All hours are UTC and are the bar's CLOSE hour (label + one bar),
matching the clock the live gates evaluate.
Analysis tooling only — never imported by the live code.
"""
from __future__ import annotations
import argparse
import math
import os
import numpy as np
import pandas as pd
try: # package import (backtest.research.burst_map)
from backtest.research.book_report import base_metrics, load_book, session_of
except ImportError: # run directly as a script from anywhere
from book_report import base_metrics, load_book, session_of
SLOT_UNKNOWN = ""
# ─────────────────────────── bar primitives ─────────────────────────
def load_bars(parquet: str, start: str | None = None,
end: str | None = None) -> pd.DataFrame:
"""Read a cached OHLCV parquet, de-duplicated and oldest-first."""
df = pd.read_parquet(parquet).sort_index()
df = df[~df.index.duplicated(keep="last")]
if start:
df = df[df.index >= pd.Timestamp(start, tz="UTC")]
if end:
df = df[df.index <= pd.Timestamp(end, tz="UTC")]
return df
def wilder_atr(df: pd.DataFrame, period: int = 14) -> pd.Series:
"""Wilder ATR — identical to the inline formula in burst_profile.py."""
h, l, c = df["High"], df["Low"], df["Close"]
tr = pd.concat([h - l, (h - c.shift()).abs(), (l - c.shift()).abs()],
axis=1).max(axis=1)
return tr.ewm(alpha=1 / period, min_periods=period, adjust=False).mean()
def burst_ratio(df: pd.DataFrame, period: int = 14) -> pd.Series:
"""Bar range expressed in PRIOR-bar ATRs."""
return (df["High"] - df["Low"]) / wilder_atr(df, period).shift(1)
def burst_flags(df: pd.DataFrame, k: float, period: int = 14) -> pd.Series:
"""True where the bar's range is >= k x the prior bar's Wilder ATR(14)."""
return (burst_ratio(df, period) >= k).fillna(False)
def burst_dir(df: pd.DataFrame) -> pd.Series:
"""+1 up bar, −1 down bar, 0 doji (sign of Close − Open)."""
return np.sign(df["Close"] - df["Open"]).astype(float)
def close_times(df: pd.DataFrame, bar_minutes: int) -> pd.DatetimeIndex:
"""Bar CLOSE timestamps (the parquet index labels bar OPEN time)."""
return df.index + pd.Timedelta(minutes=bar_minutes)
# ───────────────────── when does gold move? tables ──────────────────
def hour_weekday_table(df: pd.DataFrame, k: float, bar_minutes: int = 5,
ks: tuple[float, ...] | None = None) -> dict:
"""Burst frequency and absolute-move size on a (weekday, close-hour) grid.
Returns a dict of DataFrames (several tables are wanted, one frame each):
burst_pct : weekday rows x close-hour cols, % of bars that are >= k x
bursts (one frame; per-k frames live in `burst_pct_by_k`)
absmove : (weekday, close-hour) rows; bars, mean |1 h move| and
|4 h move| in $ and in prior-bar ATR units
by_hour : the same stats collapsed to close-hour
by_weekday : the same stats collapsed to weekday
"""
ks = tuple(ks) if ks else (k,)
ct = close_times(df, bar_minutes)
hour = pd.Series(ct.hour, index=df.index)
wd = pd.Series(ct.weekday, index=df.index)
atr_prev = wilder_atr(df).shift(1)
per_hour = max(1, int(round(60 / bar_minutes)))
c = df["Close"]
mv1 = (c - c.shift(per_hour)).abs()
mv4 = (c - c.shift(4 * per_hour)).abs()
stats = pd.DataFrame({
"weekday": wd, "hour": hour,
"abs1h_usd": mv1, "abs4h_usd": mv4,
"abs1h_atr": mv1 / atr_prev, "abs4h_atr": mv4 / atr_prev,
})
for kk in ks:
stats[f"burst_{kk:g}x"] = burst_flags(df, kk).astype(float) * 100.0
agg = {"abs1h_usd": "mean", "abs4h_usd": "mean",
"abs1h_atr": "mean", "abs4h_atr": "mean"}
agg.update({f"burst_{kk:g}x": "mean" for kk in ks})
grid = stats.groupby(["weekday", "hour"]).agg(agg)
grid.insert(0, "bars", stats.groupby(["weekday", "hour"]).size())
by_hour = stats.groupby("hour").agg(agg)
by_hour.insert(0, "bars", stats.groupby("hour").size())
by_wd = stats.groupby("weekday").agg(agg)
by_wd.insert(0, "bars", stats.groupby("weekday").size())
pct_by_k = {kk: grid[f"burst_{kk:g}x"].unstack("hour") for kk in ks}
return {
"burst_pct": pct_by_k[k if k in pct_by_k else ks[0]].round(3),
"burst_pct_by_k": {kk: v.round(3) for kk, v in pct_by_k.items()},
"absmove": grid.round(3),
"by_hour": by_hour.round(3),
"by_weekday": by_wd.round(3),
}
def post_burst_drift(df: pd.DataFrame, k: float,
horizons: tuple[int, ...] = (6, 12, 24, 48),
bar_minutes: int = 5) -> pd.DataFrame:
"""Signed drift IN the burst's direction over the next h bars, in ATRs.
A positive mean means violent bars continue; negative means they snap
back. Split by the session bucket of the burst bar's close hour.
"""
flags = burst_flags(df, k)
d = burst_dir(df)
atr_prev = wilder_atr(df).shift(1)
ct = close_times(df, bar_minutes)
bucket = pd.Series([session_of(h) for h in ct.hour], index=df.index)
c = df["Close"]
rows = []
for h in horizons:
fwd = (c.shift(-h) - c) * d / atr_prev
sub = fwd[flags].dropna()
buckets = bucket[flags].reindex(sub.index)
rows.append({"horizon_bars": h, "bucket": "ALL", "n": int(len(sub)),
"mean_atr": round(float(sub.mean()), 4) if len(sub) else float("nan"),
"median_atr": round(float(sub.median()), 4) if len(sub) else float("nan"),
"pct_positive": round(float((sub > 0).mean() * 100), 2) if len(sub) else float("nan")})
for b, s in sub.groupby(buckets):
rows.append({"horizon_bars": h, "bucket": b, "n": int(len(s)),
"mean_atr": round(float(s.mean()), 4),
"median_atr": round(float(s.median()), 4),
"pct_positive": round(float((s > 0).mean() * 100), 2)})
return pd.DataFrame(rows)
def slot_label(weekday: int, hour: int, day: int) -> str:
"""Time-slot label for a big move.
NOTE: there is no economic-calendar file in this repo, so these are
TIME SLOTS, not confirmed events — "payroll window" means "the clock
slot US payrolls normally print in", nothing more.
"""
if weekday == 4 and day <= 7 and hour in (12, 13):
return "first-Friday 12:30 slot (payroll-report window)"
if weekday == 2 and hour in (18, 19):
return "Wednesday 18:00 slot (central-bank window)"
if hour in (12, 13):
return "12:30 slot (US data window)"
if hour == 14:
return "14:00 slot"
if hour in (21, 22, 23):
return "settlement/after-hours"
if hour in (0, 1, 2):
return "Asia open"
return SLOT_UNKNOWN
def largest_moves(df: pd.DataFrame, n: int = 40, window_minutes: int = 60,
bar_minutes: int = 5) -> pd.DataFrame:
"""The n biggest trailing-`window_minutes` close-to-close moves.
Overlapping windows are de-duplicated greedily: the largest move wins and
anything whose close time is within +/- `window_minutes` of an already
kept move is dropped, so one violent hour contributes one row.
"""
bars = max(1, int(round(window_minutes / bar_minutes)))
c = df["Close"]
move = c - c.shift(bars)
atr_prev = wilder_atr(df).shift(1)
ct = close_times(df, bar_minutes)
tbl = pd.DataFrame({
"close_time": ct, "move_usd": move.to_numpy(),
"move_atr": (move / atr_prev).to_numpy(),
}).dropna()
tbl["abs_usd"] = tbl["move_usd"].abs()
tbl = tbl.sort_values("abs_usd", ascending=False)
keep, taken = [], []
gap = pd.Timedelta(minutes=window_minutes)
for _, row in tbl.iterrows():
t = row["close_time"]
if any(abs(t - u) <= gap for u in taken):
continue
taken.append(t)
keep.append(row)
if len(keep) >= n:
break
out = pd.DataFrame(keep).reset_index(drop=True)
if out.empty:
return out
out["weekday"] = out["close_time"].dt.weekday
out["hour"] = out["close_time"].dt.hour
out["direction"] = np.where(out["move_usd"] > 0, "UP", "DOWN")
out["slot"] = [slot_label(w, h, t.day) for w, h, t
in zip(out["weekday"], out["hour"], out["close_time"])]
out["move_usd"] = out["move_usd"].round(2)
out["move_atr"] = out["move_atr"].round(2)
return out[["close_time", "weekday", "hour", "direction", "move_usd",
"move_atr", "slot"]]
# ───────────────────────── book <-> burst join ──────────────────────
def _stats(r: pd.Series) -> dict:
m = base_metrics(r)
return {"n": m["n"], "total_r": round(m["total_r"], 2),
"mean_r": round(m["mean_r"], 3), "wr": round(m["wr"], 2),
"pf": round(m["pf"], 2) if math.isfinite(m["pf"]) else float("inf")}
def join_book_to_bursts(book_df: pd.DataFrame, bars: pd.DataFrame, k: float,
within_bars: int = 12,
bar_minutes: int = 5) -> dict:
"""Split a trade book by its relationship to burst bars.
(a) `after_burst` : decision_time falls within `within_bars` bars after
the most recent burst's close (bar distance is counted POSITIONALLY on
the bar index, so weekend gaps cannot inflate it).
(b) `burst_during` : at least one burst closed inside [entry, exit).
(c) `burst_agrees` / `burst_opposes`: of (b), whether the FIRST burst that
printed during the trade pushed with or against the trade direction.
"""
ct = close_times(bars, bar_minutes)
ct_arr = ct.to_numpy()
flags = burst_flags(bars, k).to_numpy()
bdir = burst_dir(bars).to_numpy()
burst_pos = np.flatnonzero(flags)
# Last burst position at or before each bar (−1 where none yet).
last_burst = np.where(flags, np.arange(len(flags)), -1)
last_burst = np.maximum.accumulate(last_burst)
dec = book_df["decision_time"].astype("datetime64[ns, UTC]").to_numpy()
ent = book_df["entry_time"].astype("datetime64[ns, UTC]").to_numpy()
exi = book_df["exit_time"].astype("datetime64[ns, UTC]").to_numpy()
# Position of the last bar that had CLOSED by decision_time.
idx = np.searchsorted(ct_arr, dec, side="right") - 1
valid = idx >= 0
lb = np.where(valid, last_burst[np.clip(idx, 0, len(flags) - 1)], -1)
bars_since = np.where((lb >= 0) & valid, idx - lb, np.iinfo(np.int32).max)
after = pd.Series((bars_since >= 0) & (bars_since <= within_bars),
index=book_df.index)
# Bursts printing inside the open trade.
bct = ct_arr[burst_pos] if len(burst_pos) else np.array([], dtype=ct_arr.dtype)
lo = np.searchsorted(bct, ent, side="left")
hi = np.searchsorted(bct, exi, side="left")
during = pd.Series(hi > lo, index=book_df.index)
trade_sign = np.where(
book_df["direction"].astype(str).str.upper().str.startswith("B"),
1.0, -1.0)
has = hi > lo
first_dir = np.zeros(len(book_df), dtype=float)
if len(burst_pos) and has.any():
# hi > lo guarantees lo[has] indexes a real burst.
first_dir[has] = bdir[burst_pos[lo[has]]]
agrees = pd.Series(has & (first_dir * trade_sign > 0),
index=book_df.index)
opposes = pd.Series(has & (first_dir * trade_sign < 0),
index=book_df.index)
r = book_df["r_multiple"]
out = {
"k": k, "within_bars": within_bars,
"bursts_in_window": int(len(burst_pos)),
"after_burst": _stats(r[after]),
"not_after_burst": _stats(r[~after]),
"burst_during": _stats(r[during]),
"no_burst_during": _stats(r[~during]),
"burst_agrees": _stats(r[agrees]),
"burst_opposes": _stats(r[opposes]),
"median_bars_since_burst": float(np.median(
bars_since[bars_since < np.iinfo(np.int32).max])) if (
bars_since < np.iinfo(np.int32).max).any() else float("nan"),
}
return out
# ──────────────────────────── reporting ─────────────────────────────
def _fmt_stats(name: str, s: dict) -> str:
pf = "inf" if math.isinf(s["pf"]) else f"{s['pf']:.2f}"
return (f" {name:<18} n={s['n']:<5} totalR={s['total_r']:>8.2f} "
f"meanR={s['mean_r']:>7.3f} WR={s['wr']:>6.2f}% PF={pf}")
def main() -> None:
ap = argparse.ArgumentParser(description=__doc__.split("\n")[0])
ap.add_argument("--parquet", required=True)
ap.add_argument("--bar-minutes", type=int, required=True)
ap.add_argument("--k", default="3,4,5",
help="burst thresholds in prior-bar ATRs (default 3,4,5)")
ap.add_argument("--book", default=None)
ap.add_argument("--book-bar-minutes", type=int, default=None,
help="bar minutes of the BOOK (default: --bar-minutes)")
ap.add_argument("--within-bars", type=int, default=12)
ap.add_argument("--start", default=None)
ap.add_argument("--end", default=None)
ap.add_argument("--n-moves", type=int, default=40)
ap.add_argument("--out", default="backtest/research/burst_map_out")
args = ap.parse_args()
ks = [float(x) for x in args.k.split(",") if x.strip()]
os.makedirs(args.out, exist_ok=True)
bars = load_bars(args.parquet, args.start, args.end)
bm = args.bar_minutes
primary = ks[0]
lines: list[str] = []
add = lines.append
add("=" * 78)
add(f"BURST MAP {os.path.basename(args.parquet)}")
add("=" * 78)
add(f"bars : {len(bars)} ({bars.index[0]} -> {bars.index[-1]}"
f", labels = bar OPEN time)")
add(f"bar minutes : {bm} hours below are bar CLOSE hours (UTC)")
add(f"burst definition : range >= K x PRIOR bar Wilder ATR(14); "
f"K = {', '.join(f'{k:g}' for k in ks)}")
add("ATR units : always normalised by the PRIOR bar's ATR")
for k in ks:
f = burst_flags(bars, k)
add(f" K={k:g}: {int(f.sum())} burst bars "
f"({float(f.mean()) * 100:.3f}% of bars)")
tabs = hour_weekday_table(bars, primary, bm, ks=tuple(ks))
add("")
add(f"── BURST % BY (weekday, close-hour), K={primary:g} ──")
add(tabs["burst_pct"].fillna(0).round(2).to_string())
add("")
add("── MEAN ABSOLUTE MOVE BY CLOSE-HOUR (all weekdays) ──")
add(tabs["by_hour"].to_string())
add("")
add("── MEAN ABSOLUTE MOVE BY WEEKDAY ──")
add(tabs["by_weekday"].to_string())
add("")
add(f"── POST-BURST DRIFT (in the burst's own direction), K={primary:g} ──")
drift = post_burst_drift(bars, primary, bar_minutes=bm)
add(drift.to_string(index=False))
add("")
add(f"── LARGEST {args.n_moves} TRAILING-60-MIN MOVES (de-overlapped) ──")
add(" slot labels are CLOCK SLOTS — no economic calendar exists in "
"this repo, so no event is confirmed.")
moves = largest_moves(bars, n=args.n_moves, bar_minutes=bm)
add(moves.to_string(index=False))
join_rows = []
if args.book:
bbm = args.book_bar_minutes or bm
book = load_book(args.book, bbm)
add("")
add(f"── BOOK JOIN {os.path.basename(args.book)} "
f"({len(book)} trades, book bar = {bbm} min) ──")
first, last = close_times(bars, bm)[0], close_times(bars, bm)[-1]
outside = int(((book["decision_time"] < first)
| (book["exit_time"] > last)).sum())
if outside:
add(f" !! {outside} trade(s) start before / end after the bar "
f"window ({first} .. {last}) — their burst context is "
f"partial. Widen --start/--end or use a longer parquet.")
add(" CAUTION: group (b) is survivorship-shaped — a trade only sees "
"a burst if it survived long enough to still be open, so its high "
"PF is mostly 'winners last longer', not a burst edge.")
for k in ks:
j = join_book_to_bursts(book, bars, k, args.within_bars, bm)
add(f" K={k:g} ({j['bursts_in_window']} bursts in window; "
f"median bars since last burst = "
f"{j['median_bars_since_burst']:.0f})")
add(f" (a) decided within {args.within_bars} bars of a burst")
add(_fmt_stats("after_burst", j["after_burst"]))
add(_fmt_stats("not_after_burst", j["not_after_burst"]))
add(" (b) a burst printed while the trade was open")
add(_fmt_stats("burst_during", j["burst_during"]))
add(_fmt_stats("no_burst_during", j["no_burst_during"]))
add(" (c) of (b), first burst vs the trade's direction")
add(_fmt_stats("burst_agrees", j["burst_agrees"]))
add(_fmt_stats("burst_opposes", j["burst_opposes"]))
for grp in ("after_burst", "not_after_burst", "burst_during",
"no_burst_during", "burst_agrees", "burst_opposes"):
join_rows.append({"k": k, "group": grp,
"within_bars": args.within_bars,
**j[grp]})
text = "\n".join(lines)
print(text)
with open(os.path.join(args.out, "burst_map.txt"), "w") as fh:
fh.write(text + "\n")
pd.concat([v.assign(k=kk) for kk, v in tabs["burst_pct_by_k"].items()]
).to_csv(os.path.join(args.out, "hour_weekday_burst_pct.csv"))
tabs["absmove"].to_csv(os.path.join(args.out, "hour_weekday_absmove.csv"))
tabs["by_hour"].to_csv(os.path.join(args.out, "by_hour.csv"))
tabs["by_weekday"].to_csv(os.path.join(args.out, "by_weekday.csv"))
pd.concat([post_burst_drift(bars, k, bar_minutes=bm).assign(k=k)
for k in ks]).to_csv(
os.path.join(args.out, "post_burst_drift.csv"), index=False)
moves.to_csv(os.path.join(args.out, "largest_moves.csv"), index=False)
if join_rows:
pd.DataFrame(join_rows).to_csv(
os.path.join(args.out, "book_join.csv"), index=False)
print(f"\n[written] {os.path.join(args.out, 'burst_map.txt')}")
if __name__ == "__main__":
main()
+127
View File
@@ -0,0 +1,127 @@
"""Are there signs BEFORE a violent bar? Analysis tooling only — never imported by the live code.
For every bar t, features are computed from bars <= t-1 only. Targets:
next_bar : bar t itself is a burst (range >= K x prior Wilder ATR14)
next_hour : any burst in bars t .. t+H-1
Reports per-quintile burst rates (lift vs base rate) in both halves, the same
within/outside the scheduled-data hours, and an out-of-sample logistic model
(fit on half 1, scored on half 2): AUC + lift in the top 5% of scores.
python backtest/research/burst_precursors.py PARQUET BAR_MINUTES [--k 3] [--start ... --end ...]
"""
import argparse, sys
import numpy as np, pandas as pd
def wilder_atr(df, period=14):
pc = df["Close"].shift(1)
tr = pd.concat([df["High"] - df["Low"], (df["High"] - pc).abs(), (df["Low"] - pc).abs()], axis=1).max(axis=1)
return tr.ewm(alpha=1 / period, min_periods=period, adjust=False).mean()
def build(df, bar_minutes, k):
d = df.copy()
rng = d["High"] - d["Low"]
atr_prev = wilder_atr(d).shift(1) # ATR known at t-1
d["burst"] = (rng >= k * atr_prev).astype(int)
H = max(1, 60 // bar_minutes)
# burst within next H bars (t..t+H-1)
d["next_hour"] = d["burst"][::-1].rolling(H, min_periods=1).max()[::-1].astype(int)
d["next_bar"] = d["burst"]
close_t = d.index + pd.Timedelta(minutes=bar_minutes)
d["hour"] = close_t.hour
r1 = rng.shift(1)
d["squeeze6"] = r1.rolling(6).mean() / atr_prev
d["squeeze12"] = r1.rolling(12).mean() / atr_prev
d["last_bar_range"] = r1 / atr_prev
sma = d["Close"].shift(1).rolling(20).mean(); sd = d["Close"].shift(1).rolling(20).std(ddof=0)
bbw = (4 * sd / sma)
d["bbw_pct500"] = bbw.rolling(500).rank(pct=True)
v1 = d["Volume"].shift(1)
d["vol_ratio"] = v1.rolling(3).mean() / v1.rolling(50).mean()
d["drift12"] = (d["Close"].shift(1) - d["Close"].shift(13)).abs() / atr_prev
d["atr_pct500"] = atr_prev.rolling(500).rank(pct=True)
b = d["burst"].shift(1).fillna(0).astype(int)
since = np.zeros(len(d)); cnt = 200
bv = b.values
for i in range(len(d)):
cnt = 0 if bv[i] == 1 else min(200, cnt + 1)
since[i] = cnt
d["bars_since_burst"] = since
return d.dropna(subset=["squeeze6", "squeeze12", "bbw_pct500", "vol_ratio", "drift12", "atr_pct500"])
FEATS = ["squeeze6", "squeeze12", "last_bar_range", "bbw_pct500", "vol_ratio", "drift12", "atr_pct500", "bars_since_burst"]
def quintile_table(d, feat, target, half):
q = pd.qcut(d[feat], 5, labels=False, duplicates="drop")
base1 = d.loc[half == 1, target].mean(); base2 = d.loc[half == 2, target].mean()
rows = []
for b in sorted(q.dropna().unique()):
m = q == b
lo, hi = d.loc[m, feat].min(), d.loc[m, feat].max()
r1 = d.loc[m & (half == 1), target].mean(); r2 = d.loc[m & (half == 2), target].mean()
rows.append((int(b), lo, hi, int(m.sum()), 100 * r1, r1 / base1, 100 * r2, r2 / base2))
t = pd.DataFrame(rows, columns=["bin", "lo", "hi", "n", "h1_rate%", "h1_lift", "h2_rate%", "h2_lift"])
return t.round(3)
def main():
ap = argparse.ArgumentParser(); ap.add_argument("parquet"); ap.add_argument("bar_minutes", type=int)
ap.add_argument("--k", type=float, default=3.0); ap.add_argument("--start"); ap.add_argument("--end")
a = ap.parse_args()
df = pd.read_parquet(a.parquet)
if a.start: df = df[df.index >= pd.Timestamp(a.start, tz="UTC")]
if a.end: df = df[df.index < pd.Timestamp(a.end, tz="UTC")]
d = build(df, a.bar_minutes, a.k)
half = pd.Series(np.where(np.arange(len(d)) < len(d) // 2, 1, 2), index=d.index)
print(f"bars {len(d)} K={a.k} base rate next_bar {100*d.next_bar.mean():.3f}% next_hour {100*d.next_hour.mean():.2f}%")
print(f"half 1 ends {d.index[len(d)//2 - 1]}")
# 1. the clock
hr = d.groupby("hour")[["next_bar", "next_hour"]].mean() * 100
hr["lift_next_bar"] = hr["next_bar"] / (100 * d.next_bar.mean())
print("\n== THE CLOCK (bar-close hour UTC): burst rate % and lift ==")
print(hr.round(2).T.to_string())
sched = d["hour"].isin([1, 12, 13])
print(f"\nscheduled hours (01,12,13) hold {100*d.loc[sched,'next_bar'].sum()/d.next_bar.sum():.1f}% of bursts in {100*sched.mean():.1f}% of bars")
# 2. per-feature quintiles, both halves
for target in ("next_bar", "next_hour"):
print(f"\n== TARGET {target}: quintile burst rate (lift vs base) — ALL hours ==")
for f in FEATS:
t = quintile_table(d, f, target, half)
print(f"-- {f}: lifts h1 {list(t.h1_lift)} h2 {list(t.h2_lift)} (bins lo..hi {t.lo.iloc[0]:.2f}..{t.hi.iloc[-1]:.2f})")
# 3. same, OUTSIDE the scheduled hours (does anything predict the unscheduled ones?)
u = d[~sched]; hu = half[~sched]
print(f"\n== next_bar lifts OUTSIDE scheduled hours (base {100*u.next_bar.mean():.3f}%) ==")
for f in FEATS:
t = quintile_table(u, f, "next_bar", hu)
print(f"-- {f}: h1 {list(t.h1_lift)} h2 {list(t.h2_lift)}")
# 4. OOS logistic model: fit half 1, score half 2
try:
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
X = d[FEATS].copy()
for h in range(24): X[f"h{h}"] = (d["hour"] == h).astype(float)
Xc = d[FEATS].copy() # clock-free variant
for name, XX in (("features+clock", X), ("features only", Xc), ("clock only", X[[f"h{h}" for h in range(24)]])):
mu, sd = XX[half == 1].mean(), XX[half == 1].std().replace(0, 1)
Z = (XX - mu) / sd
for target in ("next_bar", "next_hour"):
m = LogisticRegression(C=0.1, max_iter=2000, class_weight="balanced").fit(Z[half == 1], d.loc[half == 1, target])
p = m.predict_proba(Z[half == 2])[:, 1]; y = d.loc[half == 2, target]
auc = roc_auc_score(y, p)
top = p >= np.quantile(p, 0.95)
print(f"OOS {name:15s} {target:9s}: AUC {auc:.3f} | top-5% of scores: burst rate {100*y[top].mean():.2f}% vs base {100*y.mean():.2f}% (lift {y[top].mean()/y.mean():.2f}x), captures {100*y[top].sum()/y.sum():.1f}% of bursts")
except ImportError:
print("sklearn missing — OOS model skipped")
# 5. direction: does the burst continue the prior 12-bar drift?
bt = d[d.burst == 1]
prior = np.sign(bt["Close"].shift(1) - bt["Close"].shift(13)) if False else np.sign(d["Close"].shift(1) - d["Close"].shift(13))[d.burst == 1]
bdir = np.sign(bt["Close"] - bt["Open"])
agree = (prior == bdir)[(prior != 0) & (bdir != 0)]
print(f"\nburst direction == prior-hour drift direction: {100*agree.mean():.1f}% of {len(agree)} bursts (50% = no information)")
if __name__ == "__main__":
main()
+111
View File
@@ -0,0 +1,111 @@
"""Do scheduled economic events move gold violently, and do they hurt OPEN trades? Analysis tooling only.
Needs backtest/cache/calendar.parquet (ff_calendar_scraper.py). Uses USD events by impact.
python backtest/research/calendar_impact.py PARQUET_BARS BAR_MINUTES BOOK.csv ARM DIST [--start --end]
"""
import argparse, sys
from pathlib import Path
import numpy as np, pandas as pd
sys.path.insert(0, str(Path(__file__).resolve().parent))
from open_trade_bursts import wilder_atr, simulate, r_of, stats, K, COST
CAL = Path(__file__).resolve().parent.parent / "cache" / "calendar.parquet"
def main():
ap = argparse.ArgumentParser(); ap.add_argument("parquet"); ap.add_argument("bar_minutes", type=int)
ap.add_argument("book"); ap.add_argument("arm", type=float); ap.add_argument("dist", type=float)
ap.add_argument("--start", default="2023-09-05"); ap.add_argument("--end", default="2026-09-05")
a = ap.parse_args(); bm = a.bar_minutes
cal = pd.read_parquet(CAL)
cal["ts"] = pd.to_datetime(cal["dateline"], unit="s", utc=True)
cal = cal[(cal.ts >= pd.Timestamp(a.start, tz="UTC")) & (cal.ts < pd.Timestamp(a.end, tz="UTC"))]
usd = cal[cal.currency == "USD"].copy()
usd["impact"] = usd["impactName"].str.lower()
print(f"calendar: {len(cal)} events in window, USD {len(usd)}: " + ", ".join(f"{k}={v}" for k, v in usd.impact.value_counts().items()))
px = pd.read_parquet(a.parquet); px = px[(px.index >= pd.Timestamp(a.start, tz="UTC")) & (px.index < pd.Timestamp(a.end, tz="UTC"))]
atr_prev = wilder_atr(px).shift(1)
rng = px["High"] - px["Low"]; burst = (rng >= K * atr_prev).fillna(False)
close = px["Close"]; idx = px.index
H1 = 60 // bm
move60 = (close.shift(-H1) - close).abs() / atr_prev # |move| over the next hour from bar OPEN-ish (label)
move15 = (close.shift(-max(1, 15 // bm)) - close).abs() / atr_prev
burst60 = burst[::-1].rolling(H1, min_periods=1).max()[::-1].astype(bool)
# event -> bar containing it (label <= ts < label + bm)
def bar_pos(ts):
p = idx.searchsorted(ts, side="right") - 1
return p if 0 <= p < len(idx) and (ts - idx[p]) < pd.Timedelta(minutes=bm) else -1
usd["pos"] = [bar_pos(t) for t in usd.ts]
usd = usd[usd.pos >= 0]
usd["bar_range_atr"] = rng.values[usd.pos] / atr_prev.values[usd.pos]
usd["burst_in_60m"] = burst60.values[usd.pos]
usd["move60_atr"] = move60.values[usd.pos]; usd["move60_usd"] = (close.shift(-H1) - close).abs().values[usd.pos]
usd["move15_atr"] = move15.values[usd.pos]
# baseline: all bars, and bars in the same close-hours without any USD high/medium event within the hour
chour = (idx + pd.Timedelta(minutes=bm)).hour
ev_bars = set()
for p in usd[usd.impact.isin(["high", "medium"])].pos: ev_bars.update(range(max(0, p - H1), min(len(idx), p + H1)))
quiet = ~np.isin(np.arange(len(idx)), list(ev_bars))
def base_row(mask, label):
return dict(group=label, n=int(mask.sum()), bar_range_atr=float((rng / atr_prev)[mask].mean()), burst_in_60m_pct=100 * float(burst60[mask].mean()),
move15_atr=float(move15[mask].mean()), move60_atr=float(move60[mask].mean()), move60_usd=float((close.shift(-H1) - close).abs()[mask].mean()))
rows = [base_row(np.ones(len(idx), bool), "ALL bars"), base_row(quiet & np.isin(chour, [12, 13, 14]), "12-14 UTC bars with NO USD event")]
for imp in ("high", "medium", "low"):
g = usd[usd.impact == imp]
rows.append(dict(group=f"USD {imp}-impact events", n=len(g), bar_range_atr=g.bar_range_atr.mean(), burst_in_60m_pct=100 * g.burst_in_60m.mean(),
move15_atr=g.move15_atr.mean(), move60_atr=g.move60_atr.mean(), move60_usd=g.move60_usd.mean()))
print("\n== HOW MUCH GOLD MOVES AROUND USD EVENTS (event bar = the bar containing the release) ==")
print(pd.DataFrame(rows).round(2).to_string(index=False))
hi = usd[usd.impact == "high"]
top = hi.groupby("name").agg(n=("ts", "size"), bar_range_atr=("bar_range_atr", "mean"), burst60_pct=("burst_in_60m", lambda s: 100 * s.mean()),
move60_usd=("move60_usd", "mean"), move60_atr=("move60_atr", "mean")).sort_values("n", ascending=False)
print("\n== HIGH-IMPACT USD EVENTS BY NAME (n >= 6) ==")
print(top[top.n >= 6].round(2).to_string())
# what share of ALL bursts sit within ±30 min of a high/medium USD event?
near = np.zeros(len(idx), bool)
for p in usd[usd.impact.isin(["high", "medium"])].pos: near[max(0, p - 30 // bm):min(len(idx), p + 30 // bm + 1)] = True
nearh = np.zeros(len(idx), bool)
for p in hi.pos: nearh[max(0, p - 30 // bm):min(len(idx), p + 30 // bm + 1)] = True
print(f"\nshare of ALL violent bars within ±30 min of a USD high-impact event: {100*burst[nearh].sum()/burst.sum():.1f}% (those bars are {100*nearh.mean():.1f}% of time)")
print(f"share within ±30 min of a USD high- or medium-impact event: {100*burst[near].sum()/burst.sum():.1f}% (bars {100*near.mean():.1f}% of time)")
# open trades through high-impact events
book = pd.read_csv(a.book); book["entry_time"] = pd.to_datetime(book.entry_time, utc=True, format="mixed")
bars = dict(O=px.Open.values, H=px.High.values, L=px.Low.values, C=px.Close.values, burst=burst.values, chour=chour.values)
pos = idx.get_indexer(book.entry_time); ok = pos >= 0; book = book[ok].reset_index(drop=True); pos = pos[ok]
ev_pos = np.sort(hi.pos.values)
# flat mask: exit at the open of the bar that CONTAINS a high-impact event... use the bar before (exit ~5-15 min before)
flat_mask = np.zeros(len(idx), bool); tight_mask = np.zeros(len(idx), bool)
for p in ev_pos:
if p - 1 >= 0: flat_mask[p - 1] = True
tight_mask[max(0, p - 1):min(len(idx), p + H1)] = True
res = {}
for v, kw in (("base", {}), ("mask_flat", dict(flat_mask=flat_mask)), ("mask_tight", dict(tight_mask=tight_mask))):
rs, fl, thru = [], [], []
for i, t in book.iterrows():
xp, held, reason, flag = simulate(bars, pos[i], t.direction, t.entry, t.sl, a.arm, a.dist, v, **kw)
rs.append(r_of(t.direction, t.entry, t.sl, xp)); fl.append(flag)
if v == "base":
j0, j1 = pos[i], pos[i] + held
k = ev_pos[(ev_pos > j0) & (ev_pos <= j1)]
thru.append(k[0] if len(k) else -1)
res[v] = (np.array(rs), np.array(fl))
if v == "base": book["ev_pos"] = thru
base = res["base"][0]; half = np.where(np.arange(len(book)) < len(book) // 2, 1, 2)
hit = book.ev_pos.values >= 0
print(f"\n== OPEN TRADES THAT HELD THROUGH A HIGH-IMPACT USD EVENT ({a.book.split('/')[-1]}) ==")
print(f"trades through an event: {hit.sum()} of {len(book)} ({100*hit.mean():.1f}%) | mean R {base[hit].mean():+.3f} vs {base[~hit].mean():+.3f} for the rest | WR {100*(base[hit]>0).mean():.0f}% vs {100*(base[~hit]>0).mean():.0f}%")
if hit.sum():
d = np.where(book.direction == "BUY", 1.0, -1.0)
m15 = np.array([ (close.shift(-max(1,15//bm)).values[p] - close.values[p]) if p >= 0 else np.nan for p in book.ev_pos])
with_ = hit & (d * m15 > 0); against = hit & (d * m15 < 0)
print(f" event's first 15-min move WITH the trade: n={with_.sum()} mean R {base[with_].mean():+.3f} | AGAINST: n={against.sum()} mean R {base[against].mean():+.3f} (WR {100*(base[against]>0).mean():.0f}%)")
print("\n== PROTECTION AROUND HIGH-IMPACT USD EVENTS ==")
b = stats(base)
for v, (r, f) in res.items():
s = stats(r); aff = f == 1
d1 = r[half == 1].sum() - base[half == 1].sum(); d2 = r[half == 2].sum() - base[half == 2].sum()
print(f"{v:11s} total {s['total']:+8.2f} Δ {s['total']-b['total']:+7.2f} PF {s['pf']:.2f} DD {s['dd']:.2f} affected {aff.sum():4d} (affected: base {base[aff].sum():+.2f} -> {r[aff].sum():+.2f}) | Δ h1 {d1:+.2f} / h2 {d2:+.2f}")
if __name__ == "__main__":
main()
+150
View File
@@ -0,0 +1,150 @@
"""What do violent bars do to OPEN trades, and can we protect them? Analysis tooling only.
Replays every trade of a book on the raw bars with the backtester's exit model (SL-first,
gap-through at open, close-only trail ratchet, 300-bar timeout, $0.65 cost) and compares:
base : reproduce the book (validation: must match r_multiple for almost every trade)
flat_data : exit at the OPEN of the first bar whose close-hour enters {12,13} UTC (US data)
flat_asia : same for close-hour {0,1} UTC (Asian open)
tight_data : entering {12,13}, hold a protective stop 0.5R below the best price for the window
tight_asia : same for {0,1}
react_exit : a burst bar (range >= 3x prior ATR14) closes AGAINST the trade -> exit next open
react_tight : after an opposing burst, protective stop 0.5R below the best price for the rest
Plus a damage table: what happens to a trade when the first burst hits it (with vs against).
python backtest/research/open_trade_bursts.py BOOK.csv PARQUET BAR_MINUTES ARM DIST
"""
import sys
import numpy as np, pandas as pd
COST = 0.65; MAX_BARS = 300; K = 3.0; PROT = 0.5
def wilder_atr(df, period=14):
pc = df["Close"].shift(1)
tr = pd.concat([df["High"] - df["Low"], (df["High"] - pc).abs(), (df["Low"] - pc).abs()], axis=1).max(axis=1)
return tr.ewm(alpha=1 / period, min_periods=period, adjust=False).mean()
def simulate(bars, i0, direction, entry, sl, arm, dist, variant, hours_flat=(), hours_tight=(), flat_mask=None, tight_mask=None):
"""bars: dict of numpy arrays; i0 = index of the entry bar. Returns (exit_price, bars_held, reason, flag)."""
O, H, L, C, burst, chour = bars["O"], bars["H"], bars["L"], bars["C"], bars["burst"], bars["chour"]
d = 1.0 if direction == "BUY" else -1.0
risk = abs(entry - sl); stop = sl; peak = entry; flag = 0
n = len(O); tight = False; tight_until = -1; exit_next = False
for j in range(i0, min(n, i0 + MAX_BARS)):
if exit_next:
return O[j], j - i0, "REACT", flag
h = chour[j]
in_flat = h in hours_flat; in_tight = h in hours_tight
if variant.startswith("flat") and in_flat:
return O[j], j - i0, "FLAT", 1
if variant == "mask_flat" and flat_mask is not None and flat_mask[j]:
return O[j], j - i0, "FLAT", 1
if variant.startswith("tight") and in_tight:
tight = True; tight_until = j; flag = 1
if variant == "mask_tight" and tight_mask is not None and tight_mask[j]:
tight = True; tight_until = j; flag = 1
# stop check (SL first, gap-through at open)
if d > 0:
if O[j] <= stop: return O[j], j - i0, "SL", flag
if L[j] <= stop: return stop, j - i0, "SL", flag
else:
if O[j] >= stop: return O[j], j - i0, "SL", flag
if H[j] >= stop: return stop, j - i0, "SL", flag
# update peak / trail at bar close
peak = max(peak, H[j]) if d > 0 else min(peak, L[j])
fav = d * (peak - entry)
if fav >= arm * risk:
stop = max(stop, peak - dist * risk) if d > 0 else min(stop, peak + dist * risk)
if tight and (variant == "react_tight" or j <= tight_until or in_tight):
stop = max(stop, peak - PROT * risk) if d > 0 else min(stop, peak + PROT * risk)
if variant.startswith("react") and burst[j] and d * (C[j] - O[j]) < 0:
flag = 1
if variant == "react_exit": exit_next = True
else: tight = True
j = min(n - 1, i0 + MAX_BARS - 1)
return C[j], j - i0, "TIMEOUT", flag
def r_of(direction, entry, sl, exit_price):
d = 1.0 if direction == "BUY" else -1.0
return (d * (exit_price - entry) - COST) / abs(entry - sl)
def stats(r):
r = np.asarray(r); eq = np.concatenate([[0.0], np.cumsum(r)]); dd = (np.maximum.accumulate(eq) - eq).max()
pf = r[r > 0].sum() / max(-r[r < 0].sum(), 1e-9)
return dict(n=len(r), total=r.sum(), mean=r.mean(), wr=100 * (r > 0).mean(), pf=pf, dd=dd)
def main():
book_csv, pq, bm, arm, dist = sys.argv[1], sys.argv[2], int(sys.argv[3]), float(sys.argv[4]), float(sys.argv[5])
book = pd.read_csv(book_csv); book["entry_time"] = pd.to_datetime(book["entry_time"], utc=True, format="mixed")
px = pd.read_parquet(pq)
atr_prev = wilder_atr(px).shift(1)
bars = dict(O=px["Open"].values, H=px["High"].values, L=px["Low"].values, C=px["Close"].values,
burst=((px["High"] - px["Low"]) >= K * atr_prev).fillna(False).values,
chour=(px.index + pd.Timedelta(minutes=bm)).hour.values)
idx = pd.Index(px.index)
pos = idx.get_indexer(book["entry_time"])
ok = pos >= 0
print(f"{book_csv}: {len(book)} trades, {ok.sum()} matched to bars; arm {arm} dist {dist}")
book = book[ok].reset_index(drop=True); pos = pos[ok]
half = np.where(np.arange(len(book)) < len(book) // 2, 1, 2)
variants = {"base": {}, "flat_data": dict(hours_flat=(12, 13)), "flat_asia": dict(hours_flat=(0, 1)),
"tight_data": dict(hours_tight=(12, 13)), "tight_asia": dict(hours_tight=(0, 1)),
"react_exit": {}, "react_tight": {}}
results = {}
for v, kw in variants.items():
rs, flags = [], []
for i, t in book.iterrows():
xp, held, reason, flag = simulate(bars, pos[i], t.direction, t.entry, t.sl, arm, dist, v, **kw)
rs.append(r_of(t.direction, t.entry, t.sl, xp)); flags.append(flag)
results[v] = (np.array(rs), np.array(flags))
base = results["base"][0]
match = np.isclose(base, book["r_multiple"].values, atol=0.02).mean()
print(f"validation: replay reproduces the book's r_multiple (±0.02) on {100*match:.1f}% of trades; "
f"sum replay {base.sum():+.2f} vs book {book.r_multiple.sum():+.2f}")
print("\n== PROTECTION VARIANTS (whole book; 'affected' = trades where the rule fired) ==")
print(f"{'variant':12s} {'n':>5s} {'affected':>8s} {'total R':>9s} {'Δ vs base':>9s} {'PF':>5s} {'maxDD':>6s} {'WR%':>5s} | affected only: base R -> variant R | Δ h1 / Δ h2")
b = stats(base)
for v, (r, f) in results.items():
s = stats(r); aff = f == 1
d1 = r[half == 1].sum() - base[half == 1].sum(); d2 = r[half == 2].sum() - base[half == 2].sum()
print(f"{v:12s} {s['n']:5d} {aff.sum():8d} {s['total']:+9.2f} {s['total']-b['total']:+9.2f} {s['pf']:5.2f} {s['dd']:6.2f} {s['wr']:5.1f} | "
f"{base[aff].sum():+8.2f} -> {r[aff].sum():+8.2f} (n={aff.sum()}) | {d1:+.2f} / {d2:+.2f}")
# damage table: first burst while open (base path), with vs against
print("\n== DAMAGE: first violent bar that hits an OPEN trade (base exits) ==")
rows = []
for i, t in book.iterrows():
d = 1.0 if t.direction == "BUY" else -1.0; risk = abs(t.entry - t.sl)
xp, held, reason, _ = simulate(bars, pos[i], t.direction, t.entry, t.sl, arm, dist, "base")
j0 = pos[i]; j1 = j0 + held
for j in range(j0 + 1, j1 + 1): # bursts after the entry bar, up to and including the exit bar
if bars["burst"][j]:
against = d * (bars["C"][j] - bars["O"][j]) < 0
mtm_before = (d * (bars["C"][j - 1] - t.entry) - COST) / risk
mtm_after = (d * (bars["C"][j] - t.entry) - COST) / risk
rows.append(dict(against=against, bars_open=j - j0, mtm_before=mtm_before, mtm_after=mtm_after,
final=r_of(t.direction, t.entry, t.sl, xp), ended_here=(j == j1), hour=bars["chour"][j],
sched=bars["chour"][j] in (0, 1, 12, 13)))
break
dmg = pd.DataFrame(rows)
if len(dmg):
g = dmg.groupby("against").agg(n=("final", "size"), bars_open_med=("bars_open", "median"),
mtm_before=("mtm_before", "mean"), mtm_after=("mtm_after", "mean"),
final=("final", "mean"), ended_on_that_bar_pct=("ended_here", lambda s: 100 * s.mean()),
in_sched_hours_pct=("sched", lambda s: 100 * s.mean()))
print(g.round(3).to_string())
ag = dmg[dmg.against]
print(f"\nopposing bursts: {len(ag)} trades ({100*len(ag)/len(book):.1f}% of the book); "
f"mean R just before {ag.mtm_before.mean():+.3f} -> just after {ag.mtm_after.mean():+.3f} -> final {ag.final.mean():+.3f}; "
f"{100*(ag.final<0).mean():.0f}% end as losers; {100*ag.ended_here.mean():.0f}% are stopped ON that bar")
print(f" those with an open PROFIT before the burst (mtm>0): n={int((ag.mtm_before>0).sum())}, final mean {ag[ag.mtm_before>0].final.mean():+.3f}, "
f"{100*(ag[ag.mtm_before>0].final<0).mean():.0f}% end as losers")
print(f" by hour of the opposing burst (n, mean final R): " + ", ".join(f"{h:02d}:{n}/{m:+.2f}" for h, n, m in
ag.groupby("hour").final.agg(["size", "mean"]).reset_index().itertuples(index=False)))
if __name__ == "__main__":
main()
@@ -0,0 +1,348 @@
# Stake-sizing campaign — 2026-09-05 (3-year Oanda window, $0.65)
## Purpose and owner decisions
Owner asked (2026-09-05): predict direction and choose order size better; protect against
trend exhaustion at signal time; map hours/days where institutions move gold violently;
backtests ≤ 3 years. Decisions taken in planning:
- Order size = a per-signal RELATIVE stake ("Stake: 0.7× normal"), not account-based units.
- Exhaustion protection = RISK LESS, DON'T SKIP (every hard exhaustion block tested since
2026-07 removed net-positive trades — the winners are the stretched runs).
- Model only M5 and M15; H1/H4 measured descriptively.
- scikit-learn / lightgbm allowed in ./venv for research only; the shipped artefact is
numpy-only JSON. Neither library enters requirements.txt.
- No tuned threshold, gate or exit changes (rule 1). Sizing is a post-decision overlay computed
by `position_sizing.build_features` / `stake_multiplier` INSIDE `analyze_from_dataframes`, so
live and backtest use the same call on the same objects (rule 4).
## Harness
`backtest/gpu_backtest.py --pair XAUUSD --start 2023-09-05 --end 2026-09-05 --spread 0.65
--workers 12`, pair-factory config (no --set/--quality/--strategies/--confirm-tfs), NumPy CPU
backend (CuPy cannot locate CUDA_PATH on this box; the 2026-09-03 pinned books were also CPU).
Cache: `oanda_XAU_USD_{5min,15min,1h,4h}_20230705_20260905.parquet`, `…_1d_20221110_20260905`.
Feature books add `--export-features` (65 `f_*` columns + `stake_multiplier` + `r_sized`; the
default export stays the pinned 63-column shape). Outputs in `backtest/research/sizing_2026-09/`
(gitignored).
## Baseline (control books, unmodified backtester, 2026-09-05 23:36–23:41 UTC)
| TF | trades | total R | WR | PF | max DD |
|---|---|---|---|---|---|
| M5 | 857 | +161.32 | 53.7% | 1.53 | 9.85 |
| M15 | 1086 | +148.69 | 59.1% | 1.37 | 15.47 |
| H1 | 255 | +70.52 | 64.3% | 1.98 | 4.43 |
| H4 | 113 | +22.68 | 58.4% | 1.63 | 5.39 |
### Behaviour-neutrality checks (2026-09-05 23:57 UTC)
- 12-month M5 canonical window (2025-08-14 → 2026-08-14, $0.65) re-run on the modified tree:
`book_diff` vs `backtest/oanda_m5_12mo_065_nopb_20260903.csv` = ΔR +0.0, Δtrades +0, and the
two CSVs have the SAME MD5 (`a28fed4becd056466c26d6a947e1a0e5`). Record-keeping fields,
`SizingConfig`, and the `--export-features` hook change nothing when the flag is off.
- Feature books: the 63 shared columns are value-identical to the control books on every TF
(`assert_frame_equal`). A column-ORDER quirk was found and fixed the same day: pandas orders
export columns by first appearance, so when the first trade of a book lacks an indicator
(the 3-year M5 book's first trade had no candlestick pattern) that indicator's columns were
appended after the extras. `trades_to_dataframe` now forces the extras to the tail; test
extended with that exact case (435 tests green).
## Pre-registered acceptance rule (written 2026-09-05, BEFORE any model fit)
A per-timeframe stake-multiplier model may ship only if, evaluated on OUT-OF-SAMPLE quarters
only (all calendar quarters after a 12-month burn-in, each predicted by a model trained
exclusively on trades whose entry time is ≥ 300 primary bars before that quarter's start,
expanding window, standardisation fit on train only), ALL of the following hold:
1. Sized total R > flat total R in BOTH halves of the OOS period (split at the median OOS
entry time).
2. Sized total R > flat total R over the full OOS period.
3. Sized max drawdown (0-prepended cumulative sized-R curve, as `compute_metrics`) ≤ flat max
drawdown × 1.10.
4. Sized OOS book has PF ≥ 1.3 and positive total R in both halves (the rule-1 floor).
5. Real OOS lift > the 95th percentile of the lift distribution from 50 label-shuffle refits
(R permuted within training windows only; test labels and walk-forward geometry unchanged).
6. Mean multiplier over OOS trades within 1.00 ± 0.05; every multiplier in [0.5, 1.5].
Plus the mechanical check: the `--sizing` full-window replay reproduces the flat book's
entry/exit/`r_multiple` columns exactly (sizing changes no entries).
Model family and size map (fixed): primary = ridge regression on winsorised R (clip to
[−1.5, +4]); secondary = L2 logistic on R>0 (ships only if the primary fails and it passes);
simplicity baseline = one-feature model on `f_is_buy` (the shipped model must beat its OOS
lift, else the simpler one ships if it passes); upper bound = lightgbm on all `f_*` (never
ships; if its OOS lift is inside the shuffle null, stop); ablation = exhaustion-block-only model.
Map: `m = 0.5 + 1.0 × clip((z − q10)/(q90 − q10), 0, 1)` with q10/q90 from TRAIN predictions,
times a normaliser making the mean train multiplier 1.0, clipped to [0.5, 1.5], rounded to 0.05.
`SHIPPED_FEATURES` (≤ 12 per TF, ≥ 2 from the exhaustion block) are chosen from the Stage-B
sign-stable screen on the 3-year feature books and recorded below BEFORE the first fit; any
later change restarts the shuffle control. If any clause fails on a timeframe, that timeframe
ships nothing. No widening the clip, changing the split, dropping the purge, re-picking
features or adding quarters to rescue a result.
## Pre-registered feature lists
Derived MECHANICALLY (script in the session log, saved as `sizing_2026-09/preregistered_features.json`):
sign-stable `f_*` features from the Stage-B screen (top-minus-bottom quintile mean-R gap has the
same sign in both halves, both halves ≥ 10 trades per bucket), zero-variance columns removed,
ranked by min(|gap_h1|, |gap_h2|), top 12; the ≥ 2-from-the-exhaustion-block clause was already
satisfied on both TFs. Written 2026-09-06 00:05 UTC, before any Stage-D fit.
- **M5 (29 sign-stable candidates):** f_sl_atr, f_confidence, f_rsi, f_mtf_c1, f_swing_dist_atr,
f_ret_48_atr, f_tag_no_pullback, f_rsi_dir, f_first_friday, f_mtf_agreement,
f_exhaustion_count, f_ret_24_atr. Exhaustion-block members: f_swing_dist_atr,
f_tag_no_pullback, f_exhaustion_count (all three run POSITIVE in both halves).
- **M15 (37 sign-stable candidates):** f_rsi, f_bb_pos, f_day_of_month, f_tag_ema200_near,
f_atr_pct, f_is_buy, f_confidence, f_close_pos_in_bar, f_tag_weak_pullback, f_tag_exhaustion,
f_penalty_total, f_ret_6_atr. Exhaustion-block members: f_tag_exhaustion, f_penalty_total
(both NEGATIVE in both halves; f_tag_weak_pullback and f_tag_ema200_near also negative).
Exhaustion-only ablation list (fixed by the plan): f_swing_dist_atr, f_ema21_dist_atr,
f_ema200_dist_atr, f_exhaustion_count, f_penalty_total, f_tag_bb_overextended,
f_tag_rsi_extended, f_tag_no_pullback, f_tag_exhaustion.
## Stage B — book breakdown (3-year feature books)
Tool: `backtest/research/book_report.py` (outputs in `sizing_2026-09/report_{5min,15min,1h,4h}/`).
Halves = median TRADE split (M5 cut 2025-12-12 because trade density doubled in the last year;
M15 cut 2025-04-25; H1 2025-04-09).
**M5 (857 trades, +161.3R, PF 1.53, DD 9.85; halves +51.0R PF 1.33 / +110.3R PF 1.74; max 11
wins / 8 losses in a row; longest flat stretch 88 trades).** Mix: GENERIC 469, XAU_CP 388, BUY
765 / SELL 92 (SELLs +11.8R, positive both halves). Grade: A +0.26R/trade (227) > B +0.17 (608)
> C −0.01 (21). Session (decision hour): NY 17–22 best (+0.31R/trade, PF 1.88), overlap 13–16
weakest (+0.08, PF 1.22) — both stable. Weekday: Wed and Fri best and stable; Mon/Tue flip sign.
Feature screen: 29 sign-stable `f_*` features. **Every stretch measure is POSITIVE**: wider
stop in ATR (+0.22/+0.23 gap), higher RSI (+0.20/+0.19), farther from the last swing
(+0.18/+0.27), 48-bar run-up (+0.18/+0.20), the `no_pullback` tag (+0.17/+0.28), more exhaustion
tags (+0.12/+0.24). Only ATR-relative-to-500-bar-mean is negative (−0.11/−0.12: entries during
a volatility spike do worse). "Risk less when stretched" has NO support on M5 — a model that
honours the data will size stretched M5 trades UP, not down.
**M15 (1086 trades, +148.7R, PF 1.37, DD 15.47; halves +65.8R / +82.9R; max 12 wins / 8 losses;
flat stretch 131 trades).** All GENERIC; BUY 739 (+158.4R, PF 1.62) vs SELL 347 (−9.7R, PF 0.93;
−19.2R in half 1, +9.4R in half 2). Grade NOT monotone: C +0.45R/trade (82, PF 2.41) > A +0.14
(481) > B +0.08 (516) — stable in both halves (matches the 2026-08-04 quant review). Session:
Asian and NY ≈ +0.18–0.20R/trade, London and overlap ≈ +0.09–0.10 (all stable). Weekday:
Wednesday flat/negative in both halves (229 trades, −0.6R); Thu/Fri flip sign between halves.
Feature screen: 37 sign-stable features. Here the exhaustion block runs NEGATIVE: `exhaustion`
tag (−0.19/−0.10), `weak_pullback` tag (−0.24/−0.11), `ema200_near` tag (−0.17/−0.32), penalty
total (−0.16/−0.10), ATR percentile (−0.20/−0.15); while raw RSI (+0.31/+0.24) and Bollinger
position (+0.36/+0.22) are positive. So on M15 "stretched by the aggregator's own penalty
logic" is bad while "strong momentum reading" is good — the sizing model gets both.
**H1 (255 trades, +70.5R, PF 1.98, DD 4.43; halves +28.7R PF 1.79 / +41.9R PF 2.17)** and
**H4 (113 trades, +22.7R, PF 1.63, DD 5.39)** are descriptive only. H1 is all-BUY by
construction (shorts blocked 2026-09-03); Mon/Thu strongest, Friday the only losing weekday
(−1.4R/50, flips between halves); 03:00 UTC decision hour is the best single hour (+20.0R/40,
stable) and 04:00–05:00 negative.
**Labelling quirk (pre-existing, not fixed):** 33 M5 trades carry the backtester `strategy`
tag GENERIC while their setup strategy is XAU_ATR (`f_strat_atr = 1`); the CSV tag and the
setup's own string disagree for the ATR path. Labels only — R is unaffected.
**Stop-trigger check: PASSED** (both TFs have many sign-stable features) → proceed to Stage D.
## Stage C — institutional burst / time map
Tool: `backtest/research/burst_map.py` on the 3-year raw bars (`burst_5min/`, `burst_15min/`
under `sizing_2026-09/`). Burst = bar range ≥ K × prior Wilder ATR(14). Base rate at K=3:
0.86 % of M5 bars (1840), 1.19 % of M15 bars (846).
**Where the violent bars cluster (bar-CLOSE hour, UTC, K=3 burst % of bars):**
- 12:00–13:00 (US data window): M5 2.4 % / 2.3 %, M15 3.7 % / 5.0 % — 3–4× the base rate;
Friday 12:00–13:00 is the densest weekday cell on M5 (3.4 % / 3.1 %). Largest mean 1-hour
dollar moves are at close-hours 13–15 (M5: $7.8 / $9.4 / $8.5 vs ≈$5 elsewhere).
- 01:00 (Asia open): M5 2.9 %, M15 3.3 %. Also 05:00 (Tokyo/London hand-over) ≈1.2 % on M5.
- 22:00–23:00 (post-settlement thin book): ≈1.3–1.4 % on M5; Sunday reopen 22–23 h is the most
burst-prone slot of all (2.4–3.0 % of bars) but the calendar already blocks Sunday entries
before 18:00 NY.
- Quietest hours: 03:00–04:00 and 16:00–17:00 (0.14–0.3 %).
**Post-burst drift (in the burst's own direction, ATR units):** ≈ 0 at every horizon on both
TFs (M5 6/12/24/48-bar means −0.01/−0.01/+0.09/+0.05; M15 −0.08/+0.14/+0.18/+0.17; % positive
46–51 %). London-session bursts on M5 show a mild short-horizon reversal (42 % positive at 6
bars). No continuation or reversal edge worth a rule.
**Do bursts hurt the system's entries? No.** Trades DECIDED within 12 bars after a K=3 burst:
M5 149 trades, mean +0.199R, PF 1.62 vs +0.186R / PF 1.52 for the rest; M15 242 trades,
+0.160R / PF 1.47 vs +0.130R / PF 1.34. At K=4/5 the after-burst group is better still
(M5 +0.276R PF 1.98; M15 +0.259R PF 1.80). This re-confirms the 2026-07-30 big-candle-gate
rejection with 3 years of data: expansion bars are part of the entry signal, not a danger.
(Group (b) "burst while open" is survivorship-shaped — winners live longer — and is not evidence.)
**The 40 largest 1-hour moves** are almost all one episode (2026-01-28 → 2026-02-06, incl. a
−$394 hour at 15:45 UTC on 2026-01-29, and 2026-03-23), not recurring clock slots; only 3 of 40
fall in the 12:30 data slot and 1 in the Wednesday 18:00 central-bank slot (slot labels only —
no economic calendar file exists in the repo). Violent moves are episodic regime events; the
recurring, predictable clustering is the 12:00–15:00 UTC US-data window and the 01:00 Asia open.
**Use in this campaign:** `f_bars_since_burst` and the hour/session features are already in the
feature set, so Stage D can weigh them as a stake input; no entry gate is proposed (rule 1).
## Stage D — walk-forward results (2026-09-06 00:20 UTC)
Driver: `backtest/research/sizing_model.py` (numpy closed-form ridge + Newton logistic; lightgbm
4.7 + scikit-learn 1.9 installed in ./venv for the upper bound only; NOT in requirements.txt).
Outputs: `sizing_2026-09/sm_{m5,m15}_{ridge,logistic}/` (summary.txt, oos_trades.csv,
quarters.csv, coefs_by_quarter.csv, shuffle_lifts.csv). Validation of the driver: planted
synthetic effect recovered (lift +116.7R, shuffle rank 100 %, all clauses PASS); pure-noise
book correctly FAILS (rank 32 %); ridge matches `lstsq` to 2e-16; logistic gradient < 3e-15.
**M5 (5 OOS quarters 2025Q3→2026Q3, 619 trades — 2024Q4–2025Q2 skipped by the pre-registered
min-train rule because only 96/119/142 trades precede them; 2026Q1 alone is +92R of the +125R
flat OOS):**
| model | flat R | sized R | lift | PF | DD | mean m | shuffle p95 | rank | verdict |
|---|---|---|---|---|---|---|---|---|---|
| ridge (primary) | +125.41 | +121.98 | −3.43 | 1.63 | 10.69 | 0.891 | +13.69 | 48 % | FAIL (1,2,5,6) |
| logistic | +125.41 | +128.40 | +2.99 | 1.66 | 8.70 | 0.900 | +14.39 | 62 % | FAIL |
| isbuy baseline | +125.41 | +123.22 | −2.20 | 1.55 | 11.92 | 1.000 | +7.05 | 4 % | FAIL |
| exhaustion-only | +125.41 | +118.03 | −7.38 | 1.59 | 10.88 | 0.918 | +14.36 | 28 % | FAIL |
| lightgbm (all 65) | +125.41 | +116.05 | −9.36 | 1.51 | 11.64 | 1.007 | +7.60 | 10 % | FAIL |
Ridge weights (standardised, 5 quarters): `f_tag_no_pullback` +0.018 sign-stable 5/5 — the
"no pullback" flag INCREASES the stake; `f_swing_dist_atr` and `f_exhaustion_count` ≈ 0 and
sign-unstable. Exhaustion-only sizing LOSES 7.4R. **A clean null: nothing in the 65 features
predicts M5 trade outcomes beyond noise, linear or non-linear.**
**M15 (8 OOS quarters 2024Q4→2026Q3, 762 trades; flat +115.93R PF 1.42 DD 13.73):**
| model | sized R | lift | PF | DD | mean m | shuffle p95 | rank | clauses |
|---|---|---|---|---|---|---|---|---|
| ridge (primary) | +142.04 | +26.11 | 1.53 | 14.46 | 0.980 | +28.50 | 88 % | 1,2,3,4,6 PASS · **5 FAIL** |
| logistic | +144.78 | +28.84 | 1.53 | 13.97 | 1.002 | +31.89 | 92 % | 1,2,3,4,6 PASS · **5 FAIL** |
| isbuy baseline | +139.38 | +23.45 | 1.50 | 16.44 | 1.007 | +24.47 | 58 % | FAIL (3, 5) |
| exhaustion-only | +123.83 | +7.90 | 1.45 | 14.85 | 1.010 | +19.77 | 70 % | FAIL |
| lightgbm (all 65) | +123.87 | +7.94 | 1.44 | 14.74 | 0.996 | +17.66 | 60 % | FAIL |
Ridge weights sign-stable in 8/8 quarters: `f_tag_exhaustion` −0.022, `f_tag_ema200_near`
−0.033 (both SHRINK the stake — the "risk less when stretched" direction), `f_rsi` +0.028,
`f_bb_pos` +0.030, `f_is_buy` +0.029, `f_close_pos_in_bar` +0.022, `f_ret_6_atr` +0.022;
λ = 1000 every quarter (heavy shrinkage — the score is weak and the map is essentially a rank
transform). Sized R beat flat in 6 of 8 quarters and in both halves, BUT the label-shuffle
null has p95 = +28.5R: a random re-weighting of these 12 pre-screened features produces the
same lift 12 % of the time, and the one-feature BUY/SELL baseline already delivers +23.5R of
the +26.1R. **The model's information beyond "buys are better than sells" (≈ +2.7R over 762
trades) is indistinguishable from noise.**
## Verdict (pre-registered rule applied, no rescue attempted)
**FAIL on both timeframes → no stake multiplier ships.** No artifact, no `sizing.enabled`, no
alert line, no `--sizing` backtester flag. What remains in the tree is behaviour-neutral
research plumbing: `position_sizing.build_features` (65 decision-time features; the multiplier
path returns 1.0 and no artifact exists), `SignalReport.penalty_tags/penalty_total/
momentum_fading` (record-keeping), `TradeSetup.stake_multiplier` (None; reserved), the
`--export-features` flag, `book_report.py`, `burst_map.py`, `sizing_model.py`, and Module 29
tests. `SizingConfig` was REMOVED again before merge (rule 8: nothing reads it).
**What the owner should take from this (plain language):**
1. On the 5-minute chart there is nothing to size on. Stretched trades are the good ones; the
system's exhaustion flags carry no usable information about the outcome; a flexible model
given all 65 facts does worse than flat sizing out of sample.
2. On the 15-minute chart the only robust sizing signal is BUY vs SELL, which the gates already
act on. The exhaustion flag does point the "risk less" way in every one of 8 unseen
quarters, but the money it would make (≈ +3R per year) is inside noise at 1,086 trades.
3. "Predict direction better" with a learned second opinion is closed at this sample size:
the strongest non-linear model (all features, gradient boosting) fails on both charts.
**Open avenues (not rescues — each needs a NEW pre-registration):**
- Re-run the M15 exhaustion-tag sizing when the book has ≈ 2× the trades (≈ 2028) or on a
sealed forward sample; the sign is stable, the magnitude is not proven.
- The shuffle null is wide because the size map rank-stretches a weak score to 0.5–1.5×; a
narrower map (e.g. 0.75–1.25×) would shrink both lift and null — pre-register before testing.
- Session and weekday structure (M5: NY best / overlap weakest; M15: Wednesday flat) is stable
across halves but was rejected as a gate repeatedly; as a *size* input it was inside the
12-feature model and did not survive the null.
## Addendum 2026-09-06 — violent bars: signs before them, and protecting OPEN trades
Owner follow-up: "we might have an order running for hours before a violent move; we want
protection." Tools: `backtest/research/burst_precursors.py` (signs before a burst) and
`backtest/research/open_trade_bursts.py` (exact replay of every book trade on raw bars with the
backtester's exit model — reproduces `r_multiple` on 100 % of trades — plus protection variants).
Outputs: `sizing_2026-09/burst_precursors_{m5,m15}.txt`, `open_trade_bursts_{m5,m15}.txt`.
**Signs before a burst (both TFs, both halves stable):** the clock (close-hours 01, 12, 13 UTC
hold ~40 % of bursts in 13 % of bars); tick activity in the last 3 bars vs 50-bar mean (top
quintile ≈ 2× rate); the last bars already wider than normal (top quintile ≈ 2.2×; the
tightest quintile ≈ 0.3–0.5× — NOT calm-before-storm at bar scale); a move already under way
over the past hour (≈ 1.6×); a burst minutes ago (≈ 1.4×). Out of sample (fit half 1, score
half 2) a logistic model reaches AUC 0.75–0.77 for the next bar; the top 5 % of moments carry a
3.4–5.5 % burst chance vs 0.8–1.1 % base, and 24 % vs 9 % (M5) for "within the next hour".
Burst direction continues the prior hour's drift 61 % of the time. Signs say "ripe", not "now".
**Damage to open trades (base exits):** M5 — 164 trades (19 %) meet an opposing burst while
open: mean open R +0.42 just before → −0.06 just after → final +0.01; 43 % are stopped on that
bar; 54 % end as losers. M15 — 246 trades (23 %): +0.27 → −0.20 → final −0.12; 39 % stopped on
the bar. Trades already in profit before the burst still finish positive on average (M5 +0.53,
M15 +0.27). Bursts WITH the trade (n 176 / 246) finish at +1.11 / +0.83R — they are where the
winners come from.
**Protection variants — ALL negative on both TFs, both halves:**
| variant | M5 ΔR (h1/h2) | M15 ΔR (h1/h2) |
|---|---|---|
| flatten before US-data hours (12–13 UTC) | −31.9 (−13.6/−18.3) | −103.0 (−67.4/−35.6) |
| flatten before Asian open (00–01 UTC) | −38.7 (−22.0/−16.7) | −73.0 (−42.0/−30.9) |
| protective 0.5R stop through the data hours | −29.6 (−12.3/−17.3) | −82.3 (−55.5/−26.8) |
| protective 0.5R stop through the Asian open | −19.2 (−9.5/−9.7) | −35.1 (−17.3/−17.8) |
| exit next open after an opposing burst | −19.1 (−6.5/−12.5) | −6.9 (−11.9/+5.0) |
| 0.5R protective stop after an opposing burst | −9.0 (−3.0/−6.0) | −8.3 (−15.5/+7.3) |
The existing stop/trail already contains the opposing-burst damage (≈ 0.4–0.5R per hit); every
pre-emptive rule removes far more from the with-trade bursts than it saves. **Verdict: no
open-trade protection rule ships; do not re-propose flatten/tighten-before-window or
react-to-burst exits without a new mechanism.** Display-only "violent-move conditions" note on
alerts is the only harmless use of the precursor signs.
## Addendum 2026-09-06 — calendar events, indicator states and time slots before violent bars
Tools: `backtest/research/calendar_impact.py` (needs `backtest/cache/calendar.parquet`, scraped
2023-09-01→2026-09-10 with `ff_calendar_scraper.py`: 15,124 events, 1,055 high-impact USD in
window) and `backtest/research/burst_indicator_map.py` (states from the backtester's own
`precompute_indicators`, taken at t−1). Outputs `sizing_2026-09/calendar_impact_{m5,m15}.txt`,
`burst_indicator_map_{m5,m15}.txt`. Burst = range ≥ 3× prior ATR14; lifts quoted are stable in
BOTH halves on both TFs unless noted.
**Calendar.** The M5 bar containing a high-impact USD release averages 4.4× normal range; 54 %
of such releases produce a burst within the hour (base 9 %). By event: Non-Farm Payrolls (with
Unemployment Rate / Avg Hourly Earnings) 97 %, 9.3× bar, $11.9 mean 1-h move; CPI 91 %, 8.3×,
$9.0; FOMC decision/statement 87.5 %, 7.1×, $18.2 (press conference a further $20 hour); FOMC
projections 92 %; Retail Sales 64 %; PPI 54 %; Core PCE 50 %; Unemployment Claims 44 %; GDP
45–58 %; ISM PMIs 26 %; Powell speaks 33 %; Trump speaks 21 %; UoM sentiment ≈ 5 %. 20 % (M5)
/ 31 % (M15) of ALL bursts fall within ±30 min of a high-impact USD release (3.7–4.2 % of time);
in the 30 min AFTER one, next-bar burst odds are 8–17 % (7–16× base); 15 min before, 2–3 %.
Open trades that held through a high-impact release did BETTER (M5 +0.41R vs +0.16R; M15 +0.27
vs +0.07), even when the first 15-min move went against them (M5 +0.09, M15 +0.22). Flattening
5 min before a release: −5.4R (M5) / −59.3R (M15); 0.5R protective stop through the hour:
−2.4R / −43.1R. **Events are the biggest single cause of violent bars and still not a reason
to touch an open trade.**
**Indicator states before a burst (next-bar lift, M5 / M15).** Price OUTSIDE the Bollinger
bands: below 2.9× / 1.8–2.9×, above 2.1× / 2.2×, middle 0.2–0.8 band 0.7×. RSI < 30: 2.6–2.9× /
1.8–2.2×; RSI ≥ 70: 1.9–2.1× / 1.6–1.9×; RSI 40–60: 0.7× (calmest). Close outside the Keltner
channel: 1.4–2.2×; inside 0.76×. ADX RISING > 3: 1.4–1.5× / 1.2×; ADX FALLING < −3: 0.5–0.6×
(safest state); the ADX level itself carries nothing. MACD histogram negative and still falling
1.3× / 1.2–1.4×; histogram fading (either sign) 0.8×. Fully bearish 4-EMA stack 1.3–1.4× /
1.3×; "mostly bull" 0.8× on M5. Regime: transitioning 1.4–1.6× / 1.5–1.6×, ranging 0.88×.
Volatility regime: extreme 1.5–1.6× / 1.4–1.9×, low 0.6–0.8×. ROC beyond ±0.3 %: 1.2–1.9×
(>1 %: 1.3–10×, n small). Bollinger-width squeeze: 0.59× on M5 (LESS likely next bar — narrow
bands are not a warning); TTM squeeze (bands inside Keltner) 1.6× on M15 only. EMA200 distance,
day of month, first-Friday-as-a-day: nothing stable (first Friday matters only through its
12:30 slot).
**Time slots (M5 next-bar burst %, base 0.86 %):** Fri 01:xx 3.4 %, Fri 12:xx 3.4 %, Tue/Wed
01:xx 3.2 %, Thu 12:xx 3.1 %, Fri 13:xx 3.1 %, Sun 22:xx reopen 3.0 %, Thu 01:xx 2.7 %, Wed
12:xx 2.6 %, Sun 23:xx 2.4 %, Tue 00:xx 2.3 % (all 2.5–3.9× base). Quietest: 03–04 and 16–17.
**Direction of the burst (share going UP):** fully-bull stack 62 % (M5) / 50 % (M15); fully-bear
stack 26 % / 29 %; RSI < 30 → 22 % / 30 % up (i.e. 70–78 % DOWN); RSI ≥ 70 → 64 % / 52 % up;
MACD neg & falling → 30 % / 34 % up; trending_down regime → 26 % / 26 % up; trending_up → 60 % /
46 % up. Bearish states foretell a DOWN burst reliably on both TFs; bullish states foretell an
UP burst on M5 only.
**Reading:** violent bars arrive when price is already at an extreme with momentum accelerating
(outside the bands, RSI extreme, ADX rising, MACD still falling, transitioning regime, extreme
volatility) and/or a release is minutes away — "the storm is already brewing", not "calm before
the storm". Calm mid-band bars with fading momentum and falling ADX are the safest. None of it
changes trade handling (protection tests above); it is display material for a warning line.
## Addendum 2026-09-06 — violent-move WATCH rule (display-only feature, `violent_watch.py`)
Owner asked for a warning. Delivery: a line on signal alerts + a standalone heads-up to the owner
(once per hour per trigger set, daily cap), dark by default (`violent_watch` config block).
The first draft rule ("hot" = calendar | ≥3 indicator readings | clock & ≥2 readings) was
evaluated bar-by-bar on 3 years (`backtest/research/watch_rule_eval.py`, scratch variant
sweep): it flagged 27 % of bars hot at only ≈2× odds and would have sent ≈13 heads-ups/day on
M5 — rejected. Sweep of 11 variants (M5 / M15, both halves):
| hot definition | bars hot | next-bar lift h1/h2 | next-hour rate | hot hours/day |
|---|---|---|---|---|
| calendar only (hi-impact USD ±30 min) | 3.4 % | 6.7×/4.4× · 9.1×/6.4× | 31 % · 28 % | 2.6 |
| clock slot only | 18.8 % | 3.1×/2.0× · 3.6×/2.6× | 16 % · 13 % | 5.3 |
| ≥2 price-extreme readings | 14.6 % | 2.3×/2.2× · 1.8×/2.0× | 11 % · 6 % | 10.7 |
| **calendar ∨ (clock ∧ ≥1 price-extreme reading)** | **8.6 %** | **4.9×/3.1× · 5.9×/4.2×** | **21 % · 17 %** | **4.1 · 3.5** |
| draft (cal ∨ readings≥3 ∨ clock∧readings≥2) | 26.9 % | 2.4×/2.0× · 2.4×/2.1× | 13 % · 8 % | 13.0 |
**Shipped rule:** price-extreme readings = outside the Bollinger bands, RSI < 30 or ≥ 70, close
outside the Keltner channel (the neutral readings — ADX rising, extreme volatility, transitioning
— are dropped from the level logic; they are common and only ≈1.5×). `hot` = high-impact USD
release within 30 min or 30 min after, OR a hot clock slot with ≥ 1 price-extreme reading; `watch`
= clock slot alone or ≥ 2 price-extreme readings (≈2× odds). Direction sentence from the bear/bull
polarity of the readings plus MACD acceleration and trend regime ("DOWN" on both TFs, "UP" on M5
only). Heads-up only at `hot`. Expected ≈4 hot hours/day → ≤4 heads-ups/day after per-hour dedup.
+961
View File
@@ -0,0 +1,961 @@
"""Stage-D stake-sizing model driver (2026-09 campaign).
Walk-forward fit + evaluation of a per-signal stake multiplier on a trade
book exported by ``gpu_backtest.py --export-features`` (63 base columns +
``stake_multiplier``/``r_sized`` + 65 ``f_*`` decision-time features).
The protocol implemented here is the one pre-registered in
``backtest/research/sizing_campaign_2026-09.md`` and is NOT to be varied to
rescue a result:
* expanding walk-forward, 12-month burn-in, calendar-quarter refits, a
``purge_bars × bar_minutes`` gap between the end of train and the start of
the test quarter, standardisation fit on train only;
* models: ridge on winsorised R (primary), L2 logistic on ``R > 0``
(secondary), a one-feature ``f_is_buy`` simplicity baseline, an
exhaustion-block-only ablation, and an optional lightgbm upper bound;
* size map ``m = round_0.05(clip(normaliser × (0.5 + clip((z−lo)/(hi−lo),0,1)),
0.5, 1.5))`` with ``lo``/``hi``/``normaliser`` calibrated on TRAIN
predictions only;
* a label-shuffle null (training labels permuted, test labels untouched);
* six acceptance clauses, printed PASS/FAIL individually.
Usage:
python backtest/research/sizing_model.py --book BOOK.csv --bar-minutes 5 \
--features f_a,f_b,... [--model ridge|logistic] [--out DIR] \
[--n-shuffle 50] [--skip-lgbm] \
[--write-artifact PATH --tf 5min --confirm-tfs 15min,1h]
Analysis tooling only — never imported by the live code.
"""
from __future__ import annotations
import argparse
import hashlib
import json
import math
import subprocess
import sys
import warnings
from pathlib import Path
from typing import Dict, List, Optional, Sequence, Tuple
import numpy as np
import pandas as pd
PROJECT_ROOT = Path(__file__).resolve().parents[2]
if str(PROJECT_ROOT) not in sys.path:
sys.path.insert(0, str(PROJECT_ROOT))
from position_sizing import FEATURE_ORDER, FEATURE_VERSION # noqa: E402 (read-only import)
# ── Pre-registered constants (flag any change loudly) ──────────────────
WINSOR_LO, WINSOR_HI = -1.5, 4.0
LAMBDA_GRID: Tuple[float, ...] = (1.0, 10.0, 100.0, 1000.0)
CV_FOLDS = 5
MIN_TRAIN = 150
MAP_Q_LO, MAP_Q_HI = 0.10, 0.90
MAP_OUT_LO, MAP_OUT_HI = 0.5, 1.5
MAP_CLIP_LO, MAP_CLIP_HI = 0.5, 1.5
MAP_ROUND_TO = 0.05
MEAN_MULT_TOL = 0.05
DD_TOL = 1.10
PF_FLOOR = 1.3
NEWTON_MAX_ITER = 50
NEWTON_TOL = 1e-8
#: The exhaustion-block ablation's fixed feature list.
EXHAUSTION_FEATURES: Tuple[str, ...] = (
"f_swing_dist_atr",
"f_ema21_dist_atr",
"f_ema200_dist_atr",
"f_exhaustion_count",
"f_penalty_total",
"f_tag_bb_overextended",
"f_tag_rsi_extended",
"f_tag_no_pullback",
"f_tag_exhaustion",
)
ISBUY_FEATURES: Tuple[str, ...] = ("f_is_buy",)
#: lightgbm never ships; running the full shuffle null through 200-tree fits
#: costs minutes for a number nobody acts on, so its null is capped here.
#: DEVIATION FROM THE PROTOCOL CONSTANT — printed in the summary when used.
LGBM_MAX_SHUFFLE = 10
LGBM_PARAMS = dict(
n_estimators=200,
learning_rate=0.03,
num_leaves=7,
min_child_samples=30,
subsample=0.8,
colsample_bytree=0.8,
reg_lambda=10,
)
_LINES: List[str] = []
def say(msg: str = "") -> None:
"""Print and tee into the summary buffer."""
print(msg)
_LINES.append(msg)
# ══════════════════════════════════════════════════════════════════════
# Book loading
# ══════════════════════════════════════════════════════════════════════
def load_book(path: Path) -> pd.DataFrame:
"""Read a feature book, parse ``entry_time`` UTC, drop NaN-R rows."""
df = pd.read_csv(path)
if "entry_time" not in df.columns or "r_multiple" not in df.columns:
raise SystemExit(f"{path}: not a trade book (no entry_time/r_multiple)")
df["entry_time"] = pd.to_datetime(df["entry_time"], utc=True, format="mixed")
n0 = len(df)
bad = ~np.isfinite(pd.to_numeric(df["r_multiple"], errors="coerce").to_numpy(float))
if bad.any():
say(f"WARNING: dropping {int(bad.sum())}/{n0} rows with non-finite r_multiple")
df = df.loc[~bad].copy()
df["r_multiple"] = df["r_multiple"].astype(float)
df = df.sort_values("entry_time", kind="mergesort").reset_index(drop=True)
return df
def resolve_features(df: pd.DataFrame, names: Sequence[str], label: str) -> List[str]:
"""Validate a requested feature list against the book and FEATURE_ORDER."""
missing = [n for n in names if n not in df.columns]
if missing:
raise SystemExit(f"{label}: columns absent from book: {', '.join(missing)}")
unknown = [n for n in names if n.startswith("f_") and n not in FEATURE_ORDER]
if unknown:
say(f"NOTE: {label}: {', '.join(unknown)} not in position_sizing.FEATURE_ORDER")
return list(names)
def sha256_file(path: Path) -> str:
h = hashlib.sha256()
with open(path, "rb") as fh:
for chunk in iter(lambda: fh.read(1 << 20), b""):
h.update(chunk)
return h.hexdigest()
def git_commit() -> str:
try:
return subprocess.check_output(
["git", "rev-parse", "HEAD"], cwd=PROJECT_ROOT, text=True).strip()
except Exception: # pragma: no cover - research tool
return "unknown"
# ══════════════════════════════════════════════════════════════════════
# Walk-forward geometry
# ══════════════════════════════════════════════════════════════════════
def quarter_starts(first: pd.Timestamp, last: pd.Timestamp) -> List[pd.Timestamp]:
"""Every Jan/Apr/Jul/Oct 1 00:00 UTC in ``[first, last]`` inclusive."""
out: List[pd.Timestamp] = []
year = first.year
while True:
for month in (1, 4, 7, 10):
ts = pd.Timestamp(year=year, month=month, day=1, tz="UTC")
if ts < first:
continue
if ts > last:
return out
out.append(ts)
year += 1
def quarter_label(ts: pd.Timestamp) -> str:
return f"{ts.year}Q{(ts.month - 1) // 3 + 1}"
def build_folds(times: np.ndarray, n_folds: int) -> List[np.ndarray]:
"""Contiguous-by-time index blocks (rows are already time-sorted)."""
n = len(times)
edges = np.linspace(0, n, n_folds + 1).astype(int)
return [np.arange(edges[i], edges[i + 1]) for i in range(n_folds)
if edges[i + 1] > edges[i]]
# ══════════════════════════════════════════════════════════════════════
# Standardisation
# ══════════════════════════════════════════════════════════════════════
def fit_standardiser(X: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Return (mean, scale, keep) from a training matrix.
NaNs are ignored when computing the mean (they are later filled with it).
``scale`` is the population std with zeros replaced by 1.0; ``keep`` is
False for zero-variance-in-train columns, which are dropped from the fit
and carry a coefficient of 0.
"""
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=RuntimeWarning)
mean = np.nanmean(X, axis=0)
std = np.nanstd(X, axis=0)
mean = np.where(np.isfinite(mean), mean, 0.0)
std = np.where(np.isfinite(std), std, 0.0)
keep = std > 0
scale = np.where(keep, std, 1.0)
return mean, scale, keep
def apply_standardiser(X: np.ndarray, mean: np.ndarray, scale: np.ndarray) -> np.ndarray:
"""Fill NaN with the TRAIN mean, then standardise."""
Z = np.array(X, dtype=float, copy=True)
bad = ~np.isfinite(Z)
if bad.any():
Z[bad] = np.take(mean, np.where(bad)[1])
return (Z - mean) / scale
# ══════════════════════════════════════════════════════════════════════
# Models (numpy only)
# ══════════════════════════════════════════════════════════════════════
def _solve(A: np.ndarray, b: np.ndarray) -> np.ndarray:
try:
return np.linalg.solve(A, b)
except np.linalg.LinAlgError: # pragma: no cover - degenerate design
return np.linalg.lstsq(A, b, rcond=None)[0]
def fit_ridge(X: np.ndarray, y: np.ndarray, lam: float) -> Tuple[np.ndarray, float]:
"""Closed-form ridge with an UNPENALISED intercept.
Centring X and y is exactly equivalent to appending a 1-column carrying
zero penalty. Returns ``(coef, intercept)`` with ``pred = intercept + X @ coef``.
"""
if X.shape[1] == 0:
return np.zeros(0), float(y.mean()) if len(y) else 0.0
xm = X.mean(axis=0)
ym = float(y.mean())
Xc = X - xm
w = _solve(Xc.T @ Xc + lam * np.eye(X.shape[1]), Xc.T @ (y - ym))
return w, ym - float(xm @ w)
def fit_logistic(X: np.ndarray, y: np.ndarray, lam: float,
max_iter: int = NEWTON_MAX_ITER,
tol: float = NEWTON_TOL) -> Tuple[np.ndarray, float, float]:
"""L2 logistic via Newton-Raphson, intercept unpenalised.
Returns ``(coef, intercept, grad_norm)`` where ``grad_norm`` is the
infinity-norm of the penalised gradient at the returned point.
"""
n, p = X.shape
A = np.hstack([np.ones((n, 1)), X])
pen = np.zeros(p + 1)
pen[1:] = lam
beta = np.zeros(p + 1)
# Warm-start the intercept at the empirical log-odds.
frac = float(np.clip(y.mean(), 1e-6, 1 - 1e-6))
beta[0] = math.log(frac / (1 - frac))
grad = np.zeros(p + 1)
for _ in range(max_iter):
eta = A @ beta
prob = 1.0 / (1.0 + np.exp(-np.clip(eta, -35, 35)))
grad = A.T @ (y - prob) - pen * beta
W = np.clip(prob * (1 - prob), 1e-10, None)
H = A.T @ (A * W[:, None]) + np.diag(pen)
step = _solve(H, grad)
beta = beta + step
if np.max(np.abs(step)) < tol:
break
eta = A @ beta
prob = 1.0 / (1.0 + np.exp(-np.clip(eta, -35, 35)))
grad = A.T @ (y - prob) - pen * beta
return beta[1:], float(beta[0]), float(np.max(np.abs(grad)))
def predict_linear(X: np.ndarray, coef: np.ndarray, intercept: float) -> np.ndarray:
if X.shape[1] == 0:
return np.full(len(X), intercept, dtype=float)
return intercept + X @ coef
def sigmoid(z: np.ndarray) -> np.ndarray:
return 1.0 / (1.0 + np.exp(-np.clip(z, -35, 35)))
def log_loss(y: np.ndarray, p: np.ndarray) -> float:
p = np.clip(p, 1e-12, 1 - 1e-12)
return float(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))
def choose_lambda(X: np.ndarray, y: np.ndarray, kind: str,
grid: Sequence[float] = LAMBDA_GRID,
n_folds: int = CV_FOLDS) -> Tuple[float, Dict[float, float]]:
"""5-fold time-contiguous CV inside the training window."""
folds = build_folds(np.arange(len(y)), n_folds)
if X.shape[1] == 0 or len(folds) < 2:
return float(grid[0]), {}
scores: Dict[float, float] = {}
for lam in grid:
vals: List[float] = []
for hold in folds:
mask = np.ones(len(y), dtype=bool)
mask[hold] = False
if mask.sum() < 2:
continue
if kind == "logistic":
ytr = y[mask]
if ytr.min() == ytr.max(): # single-class fold
p = np.full(len(hold), float(np.clip(ytr.mean(), 1e-6, 1 - 1e-6)))
else:
c, b, _ = fit_logistic(X[mask], ytr, lam)
p = sigmoid(predict_linear(X[hold], c, b))
vals.append(log_loss(y[hold], p))
else:
c, b = fit_ridge(X[mask], y[mask], lam)
resid = y[hold] - predict_linear(X[hold], c, b)
vals.append(float(np.mean(resid ** 2)))
scores[float(lam)] = float(np.mean(vals)) if vals else float("inf")
best = min(scores, key=lambda k: scores[k])
return best, scores
def lgbm_available() -> bool:
try:
import lightgbm # noqa: F401
return True
except Exception:
return False
def fit_lgbm(X: np.ndarray, y: np.ndarray, seed: int):
import lightgbm as lgb
model = lgb.LGBMRegressor(random_state=seed, verbose=-1, **LGBM_PARAMS)
model.fit(X, y)
return model
# ══════════════════════════════════════════════════════════════════════
# Size map
# ══════════════════════════════════════════════════════════════════════
def round_to(x, step: float = MAP_ROUND_TO) -> np.ndarray:
"""Round to the nearest ``step`` and kill float dust (0.05 → 2 dp)."""
decimals = max(0, int(round(-math.log10(step)))) + 1
return np.round(np.round(np.asarray(x, dtype=float) / step) * step, decimals)
def calibrate_map(z_train: np.ndarray) -> Dict[str, float]:
"""lo/hi/normaliser from TRAIN predictions only."""
lo = float(np.quantile(z_train, MAP_Q_LO))
hi = float(np.quantile(z_train, MAP_Q_HI))
if not np.isfinite(lo) or not np.isfinite(hi) or hi <= lo:
return {"lo": lo, "hi": hi, "normalizer": 1.0, "degenerate": True}
raw = MAP_OUT_LO + (MAP_OUT_HI - MAP_OUT_LO) * np.clip((z_train - lo) / (hi - lo), 0, 1)
mean_raw = float(np.mean(raw))
norm = 1.0 / mean_raw if mean_raw > 0 else 1.0
return {"lo": lo, "hi": hi, "normalizer": norm, "degenerate": False}
def apply_map(z: np.ndarray, cal: Dict[str, float]) -> np.ndarray:
z = np.asarray(z, dtype=float)
if cal.get("degenerate"):
return np.ones(len(z))
lo, hi, norm = cal["lo"], cal["hi"], cal["normalizer"]
raw = MAP_OUT_LO + (MAP_OUT_HI - MAP_OUT_LO) * np.clip((z - lo) / (hi - lo), 0, 1)
return round_to(np.clip(raw * norm, MAP_CLIP_LO, MAP_CLIP_HI))
# ══════════════════════════════════════════════════════════════════════
# Metrics
# ══════════════════════════════════════════════════════════════════════
def max_drawdown(r: np.ndarray) -> float:
"""0-prepended cumulative-R drawdown (gpu_backtest.compute_metrics)."""
r = np.asarray(r, dtype=float)
if len(r) == 0:
return 0.0
cum = np.concatenate(([0.0], np.cumsum(r)))
return float((np.maximum.accumulate(cum) - cum).max())
def profit_factor(r: np.ndarray) -> float:
r = np.asarray(r, dtype=float)
gp = float(r[r > 0].sum()) if (r > 0).any() else 0.0
gl = float(-r[r < 0].sum()) if (r < 0).any() else 0.0
return gp / gl if gl > 0 else float("inf")
def multiplier_hist(m: np.ndarray) -> pd.DataFrame:
"""Counts in 0.1-wide bins from 0.5 to 1.5.
``np.histogram`` already closes the final bin on the right, so an
exactly-1.5 multiplier lands in [1.4, 1.5] once and the counts sum to
``len(m)``.
"""
edges = np.round(np.arange(MAP_OUT_LO, MAP_OUT_HI + 1e-9, 0.1), 2)
counts, _ = np.histogram(np.asarray(m, dtype=float), bins=edges)
counts = counts.astype(int)
return pd.DataFrame({
"bin_lo": edges[:-1], "bin_hi": edges[1:], "count": counts,
"pct": np.where(len(m), counts / max(len(m), 1) * 100, 0.0),
})
# ══════════════════════════════════════════════════════════════════════
# Walk-forward driver
# ══════════════════════════════════════════════════════════════════════
def walk_forward(df: pd.DataFrame, features: Sequence[str], model: str,
bar_minutes: int, burn_in_months: int, purge_bars: int,
seed: int, y_train_override: Optional[np.ndarray] = None,
verbose: bool = True) -> Dict:
"""One full expanding walk-forward pass.
``y_train_override`` (full-book length, aligned to ``df``) replaces the
TRAINING labels only — this is how the shuffle null is run. Test labels
always come from the real ``r_multiple``.
"""
times = df["entry_time"].to_numpy()
r_all = df["r_multiple"].to_numpy(float)
r_train_src = r_all if y_train_override is None else np.asarray(y_train_override, float)
X_all = df[list(features)].apply(pd.to_numeric, errors="coerce").to_numpy(float) \
if features else np.zeros((len(df), 0))
book_start = df["entry_time"].iloc[0]
book_end = df["entry_time"].iloc[-1]
burn_end = book_start + pd.DateOffset(months=burn_in_months)
purge = pd.Timedelta(minutes=purge_bars * bar_minutes)
qs = [q for q in quarter_starts(book_start, book_end) if q >= burn_end]
rows: List[Dict] = []
per_q: List[Dict] = []
coef_rows: List[Dict] = []
skipped: List[str] = []
dropped_seen: set = set()
for q_start in qs:
q_end = q_start + pd.DateOffset(months=3)
train_mask = df["entry_time"] < (q_start - purge)
test_mask = (df["entry_time"] >= q_start) & (df["entry_time"] < q_end)
n_tr, n_te = int(train_mask.sum()), int(test_mask.sum())
if n_tr < MIN_TRAIN or n_te == 0:
skipped.append(f"{quarter_label(q_start)} (train={n_tr}, test={n_te})")
continue
tr = np.flatnonzero(train_mask.to_numpy())
te = np.flatnonzero(test_mask.to_numpy())
Xtr_raw, Xte_raw = X_all[tr], X_all[te]
mean, scale, keep = fit_standardiser(Xtr_raw)
dropped = [features[i] for i in np.flatnonzero(~keep)]
if verbose and dropped:
tag = (quarter_label(q_start), tuple(dropped))
if tag not in dropped_seen:
dropped_seen.add(tag)
say(f" WARNING {quarter_label(q_start)}: zero variance in train, "
f"dropped from fit (coef 0): {', '.join(dropped)}")
Xtr = apply_standardiser(Xtr_raw, mean, scale)[:, keep]
Xte = apply_standardiser(Xte_raw, mean, scale)[:, keep]
r_tr = r_train_src[tr]
lam: Optional[float] = None
if model == "logistic":
y_tr = (r_tr > 0).astype(float)
lam, _ = choose_lambda(Xtr, y_tr, "logistic")
coef, b, _ = fit_logistic(Xtr, y_tr, lam)
z_tr = sigmoid(predict_linear(Xtr, coef, b))
z_te = sigmoid(predict_linear(Xte, coef, b))
elif model == "lgbm":
y_tr = np.clip(r_tr, WINSOR_LO, WINSOR_HI)
gbm = fit_lgbm(Xtr, y_tr, seed)
z_tr = np.asarray(gbm.predict(Xtr), dtype=float)
z_te = np.asarray(gbm.predict(Xte), dtype=float)
coef, b = np.zeros(int(keep.sum())), 0.0
else: # ridge / isbuy / exhaustion
y_tr = np.clip(r_tr, WINSOR_LO, WINSOR_HI)
lam, _ = choose_lambda(Xtr, y_tr, "ridge")
coef, b = fit_ridge(Xtr, y_tr, lam)
z_tr = predict_linear(Xtr, coef, b)
z_te = predict_linear(Xte, coef, b)
cal = calibrate_map(z_tr)
m_te = apply_map(z_te, cal)
if model != "lgbm":
full_coef = np.zeros(len(features))
full_coef[np.flatnonzero(keep)] = coef
for name, c in zip(features, full_coef):
coef_rows.append({"quarter": quarter_label(q_start), "feature": name,
"coef_std": float(c), "lambda": lam,
"dropped": name in dropped})
coef_rows.append({"quarter": quarter_label(q_start), "feature": "_intercept",
"coef_std": float(b), "lambda": lam, "dropped": False})
r_te = r_all[te]
for i, idx in enumerate(te):
rows.append({
"entry_time": df["entry_time"].iloc[idx],
"r_multiple": float(r_te[i]),
"mult": float(m_te[i]),
"r_sized": float(r_te[i] * m_te[i]),
"quarter": quarter_label(q_start),
"z": float(z_te[i]),
})
per_q.append({
"quarter": quarter_label(q_start),
"q_start": q_start, "n_train": n_tr, "n_test": n_te,
"flat_r": float(r_te.sum()), "sized_r": float((r_te * m_te).sum()),
"mean_m": float(m_te.mean()), "lambda": lam,
"map_lo": cal["lo"], "map_hi": cal["hi"], "normalizer": cal["normalizer"],
})
oos = pd.DataFrame(rows)
return {"oos": oos, "quarters": pd.DataFrame(per_q),
"coefs": pd.DataFrame(coef_rows), "skipped": skipped,
"quarter_starts": qs}
def evaluate(oos: pd.DataFrame) -> Dict:
"""Headline OOS numbers for one model's walk-forward output."""
if oos.empty:
return {"n": 0, "flat_r": 0.0, "sized_r": 0.0, "lift": 0.0}
r = oos["r_multiple"].to_numpy(float)
rs = oos["r_sized"].to_numpy(float)
m = oos["mult"].to_numpy(float)
med = oos["entry_time"].median()
first = (oos["entry_time"] < med).to_numpy()
second = ~first
return {
"n": int(len(oos)),
"flat_r": float(r.sum()), "sized_r": float(rs.sum()),
"lift": float(rs.sum() - r.sum()),
"flat_pf": profit_factor(r), "sized_pf": profit_factor(rs),
"flat_dd": max_drawdown(r), "sized_dd": max_drawdown(rs),
"mean_m": float(m.mean()), "min_m": float(m.min()), "max_m": float(m.max()),
"median_time": med,
"n_h1": int(first.sum()), "n_h2": int(second.sum()),
"flat_h1": float(r[first].sum()), "flat_h2": float(r[second].sum()),
"sized_h1": float(rs[first].sum()), "sized_h2": float(rs[second].sum()),
}
def acceptance(ev: Dict, real_lift: float, shuffle_p95: float) -> Dict:
"""The six pre-registered clauses."""
c1 = ev["sized_h1"] > ev["flat_h1"] and ev["sized_h2"] > ev["flat_h2"]
c2 = ev["sized_r"] > ev["flat_r"]
c3 = ev["sized_dd"] <= DD_TOL * ev["flat_dd"]
c4 = (ev["sized_pf"] >= PF_FLOOR) and ev["sized_h1"] > 0 and ev["sized_h2"] > 0
c5 = real_lift > shuffle_p95 if np.isfinite(shuffle_p95) else False
c6 = (abs(ev["mean_m"] - 1.0) <= MEAN_MULT_TOL
and ev["min_m"] >= MAP_CLIP_LO - 1e-9 and ev["max_m"] <= MAP_CLIP_HI + 1e-9)
clauses = {"1_both_halves": bool(c1), "2_overall": bool(c2), "3_drawdown": bool(c3),
"4_pf_and_halves": bool(c4), "5_beats_shuffle": bool(c5),
"6_multiplier_bounds": bool(c6)}
clauses["overall"] = bool(all(clauses.values()))
return clauses
def shuffle_null(df: pd.DataFrame, features: Sequence[str], model: str,
bar_minutes: int, burn_in_months: int, purge_bars: int,
n_shuffle: int, seed: int) -> np.ndarray:
"""Lift distribution under permuted TRAINING labels (test labels real)."""
rng = np.random.default_rng(seed)
r_all = df["r_multiple"].to_numpy(float)
lifts: List[float] = []
for _ in range(n_shuffle):
r_perm = r_all[rng.permutation(len(r_all))]
res = walk_forward(df, features, model, bar_minutes, burn_in_months,
purge_bars, seed, y_train_override=r_perm, verbose=False)
ev = evaluate(res["oos"])
lifts.append(ev["lift"])
return np.asarray(lifts, dtype=float)
# ══════════════════════════════════════════════════════════════════════
# Final fit + artifact
# ══════════════════════════════════════════════════════════════════════
def final_fit(df: pd.DataFrame, features: Sequence[str], model: str,
seed: int) -> Dict:
"""Fit on ALL trades; calibrate the map on in-sample predictions."""
X_raw = df[list(features)].apply(pd.to_numeric, errors="coerce").to_numpy(float)
r = df["r_multiple"].to_numpy(float)
mean, scale, keep = fit_standardiser(X_raw)
X = apply_standardiser(X_raw, mean, scale)[:, keep]
if model == "logistic":
y = (r > 0).astype(float)
lam, _ = choose_lambda(X, y, "logistic")
coef, b, gnorm = fit_logistic(X, y, lam)
z = sigmoid(predict_linear(X, coef, b))
target = "r_positive"
else:
y = np.clip(r, WINSOR_LO, WINSOR_HI)
lam, _ = choose_lambda(X, y, "ridge")
coef, b = fit_ridge(X, y, lam)
gnorm = float("nan")
z = predict_linear(X, coef, b)
target = "r_winsorised"
full_coef = np.zeros(len(features))
full_coef[np.flatnonzero(keep)] = coef
cal = calibrate_map(z)
return {"mean": mean, "scale": scale, "keep": keep, "coef": full_coef,
"intercept": float(b), "lambda": float(lam), "z": z,
"mult": apply_map(z, cal), "map": cal, "target": target,
"grad_norm": gnorm}
def predict_multiplier(feature_dict: Dict[str, float], artifact: Dict) -> float:
"""Runtime formula: standardise → linear score → size map → 0.05 grid."""
names = artifact["features"]
mean = artifact["mean"]
scale = artifact["scale"]
coef = artifact["coef"]
z = float(artifact["intercept"])
for i, name in enumerate(names):
raw = feature_dict.get(name, mean[i])
if raw is None or not np.isfinite(float(raw)):
raw = mean[i]
z += float(coef[i]) * (float(raw) - float(mean[i])) / float(scale[i])
if artifact["model"] == "logistic":
z = 1.0 / (1.0 + math.exp(-max(-35.0, min(35.0, z))))
mp = artifact["map"]
if mp["hi"] <= mp["lo"]:
return 1.0
frac = min(1.0, max(0.0, (z - mp["lo"]) / (mp["hi"] - mp["lo"])))
raw_m = mp["out_lo"] + (mp["out_hi"] - mp["out_lo"]) * frac
val = min(mp["clip_hi"], max(mp["clip_lo"], raw_m * mp["normalizer"]))
step = mp["round_to"]
return float(round(round(val / step) * step, 2))
def build_artifact(fit: Dict, features: Sequence[str], model: str, df: pd.DataFrame,
ev: Dict, clauses: Dict, real_lift: float, shuffle: np.ndarray,
n_shuffle: int, book: Path, tf: str, confirm_tfs: List[str],
burn_in_months: int, purge_bars: int, spread: float) -> Dict:
return {
"schema_version": 1,
"feature_version": FEATURE_VERSION,
"pair": "XAU/USD",
"primary_tf": tf,
"confirmation_tfs": confirm_tfs,
"model": model,
"target": fit["target"],
"features": list(features),
"mean": [float(v) for v in fit["mean"]],
"scale": [float(v) for v in fit["scale"]],
"coef": [float(v) for v in fit["coef"]],
"intercept": float(fit["intercept"]),
"map": {
"kind": "clipped_linear_on_quantiles",
"lo": float(fit["map"]["lo"]), "hi": float(fit["map"]["hi"]),
"out_lo": MAP_OUT_LO, "out_hi": MAP_OUT_HI,
"normalizer": float(fit["map"]["normalizer"]),
"clip_lo": MAP_CLIP_LO, "clip_hi": MAP_CLIP_HI,
"round_to": MAP_ROUND_TO,
},
"training": {
"window": [df["entry_time"].iloc[0].isoformat(),
df["entry_time"].iloc[-1].isoformat()],
"spread": spread,
"n_train_trades": int(len(df)),
"n_oos_trades": int(ev.get("n", 0)),
"burn_in_months": burn_in_months,
"purge_bars": purge_bars,
"refit": "calendar_quarter",
"final_fit": "full_window",
"lambda": float(fit["lambda"]),
"git_commit": git_commit(),
"book_sha256": sha256_file(book),
},
"acceptance": {
**{k: v for k, v in clauses.items()},
"flat_r": ev.get("flat_r"), "sized_r": ev.get("sized_r"),
"sized_pf": ev.get("sized_pf"), "flat_dd": ev.get("flat_dd"),
"sized_dd": ev.get("sized_dd"),
"sized_h1": ev.get("sized_h1"), "sized_h2": ev.get("sized_h2"),
"flat_h1": ev.get("flat_h1"), "flat_h2": ev.get("flat_h2"),
"real_lift": float(real_lift),
"shuffle_p95_lift": float(np.percentile(shuffle, 95)) if len(shuffle) else None,
"n_shuffle": int(n_shuffle),
"mean_multiplier": ev.get("mean_m"),
},
}
# ══════════════════════════════════════════════════════════════════════
# Reporting
# ══════════════════════════════════════════════════════════════════════
def fmt(x, nd=2, width=0) -> str:
if x is None:
return "n/a".rjust(width)
if isinstance(x, float) and not np.isfinite(x):
return ("inf" if x > 0 else "-inf").rjust(width)
s = f"{x:.{nd}f}" if isinstance(x, (int, float, np.floating)) else str(x)
return s.rjust(width) if width else s
def print_quarters(q: pd.DataFrame) -> None:
say(f"{'quarter':<9}{'n_tr':>6}{'n_te':>6}{'flat R':>10}{'sized R':>10}"
f"{'mean m':>9}{'lambda':>9}")
for _, row in q.iterrows():
say(f"{row['quarter']:<9}{int(row['n_train']):>6}{int(row['n_test']):>6}"
f"{fmt(row['flat_r'], 2, 10)}{fmt(row['sized_r'], 2, 10)}"
f"{fmt(row['mean_m'], 3, 9)}"
f"{('n/a' if row['lambda'] is None else fmt(row['lambda'], 0)):>9}")
def print_hist(m: np.ndarray) -> None:
h = multiplier_hist(m)
say(f"{'bin':<12}{'count':>8}{'pct':>8}")
for _, row in h.iterrows():
bar = "#" * int(round(row["pct"] / 2))
say(f"{row['bin_lo']:.1f}–{row['bin_hi']:.1f}{'':<4}{int(row['count']):>8}"
f"{row['pct']:>7.1f}% {bar}")
def print_coef_signs(coefs: pd.DataFrame, label: str) -> None:
"""Mean standardised coefficient + sign stability across quarters."""
if coefs.empty:
return
d = coefs[coefs["feature"] != "_intercept"]
say(f" standardised coefficients across quarters ({label}):")
say(f" {'feature':<24}{'mean':>10}{'min':>10}{'max':>10}{'+q':>5}{'-q':>5} sign-stable")
for name, grp in d.groupby("feature", sort=False):
c = grp["coef_std"].to_numpy(float)
pos, neg = int((c > 0).sum()), int((c < 0).sum())
stable = "yes" if (pos == 0 or neg == 0) and (pos + neg) > 0 else "no"
say(f" {name:<24}{c.mean():>10.4f}{c.min():>10.4f}{c.max():>10.4f}"
f"{pos:>5}{neg:>5} {stable}")
def print_comparison(table: List[Dict]) -> None:
say(f"{'model':<12}{'n':>6}{'flat R':>10}{'sized R':>10}{'lift':>9}"
f"{'PF':>7}{'DD':>8}{'mean m':>8}{'p95 null':>10}{'rank%':>7} pass")
for row in table:
say(f"{row['model']:<12}{row['n']:>6}{fmt(row['flat_r'], 2, 10)}"
f"{fmt(row['sized_r'], 2, 10)}{fmt(row['lift'], 2, 9)}"
f"{fmt(row['sized_pf'], 2, 7)}{fmt(row['sized_dd'], 2, 8)}"
f"{fmt(row['mean_m'], 3, 8)}{fmt(row['shuffle_p95'], 2, 10)}"
f"{fmt(row['rank_pct'], 0, 7)} {'PASS' if row['pass'] else 'FAIL'}")
# ══════════════════════════════════════════════════════════════════════
# main
# ══════════════════════════════════════════════════════════════════════
def main() -> None:
ap = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("--book", required=True, type=Path)
ap.add_argument("--bar-minutes", required=True, type=int)
ap.add_argument("--features", required=True,
help="comma-separated f_* column names for the main model")
ap.add_argument("--model", default="ridge", choices=("ridge", "logistic"))
ap.add_argument("--burn-in-months", type=int, default=12)
ap.add_argument("--purge-bars", type=int, default=300)
ap.add_argument("--n-shuffle", type=int, default=50)
ap.add_argument("--seed", type=int, default=20260905)
ap.add_argument("--skip-lgbm", action="store_true")
ap.add_argument("--out", type=Path, default=None)
ap.add_argument("--write-artifact", type=Path, default=None)
ap.add_argument("--tf", default=None, help="primary TF label for the artifact")
ap.add_argument("--confirm-tfs", default="15min,1h")
ap.add_argument("--spread", type=float, default=0.65,
help="cost the book was produced at (recorded in the artifact)")
args = ap.parse_args()
if args.write_artifact and not args.tf:
raise SystemExit("--write-artifact requires --tf")
say("=" * 96)
say("Stage-D stake-sizing walk-forward — RESEARCH ONLY")
say("=" * 96)
say(f"book : {args.book}")
df = load_book(args.book)
say(f"trades : {len(df)} "
f"({df['entry_time'].iloc[0]} → {df['entry_time'].iloc[-1]})")
say(f"bar minutes : {args.bar_minutes} purge = {args.purge_bars} bars = "
f"{args.purge_bars * args.bar_minutes} min "
f"({args.purge_bars * args.bar_minutes / 1440:.2f} d)")
say(f"burn-in : {args.burn_in_months} months "
f"(train-only until {df['entry_time'].iloc[0] + pd.DateOffset(months=args.burn_in_months)})")
say(f"model / seed : {args.model} / {args.seed} shuffles: {args.n_shuffle}")
say(f"winsor / lambda : [{WINSOR_LO}, {WINSOR_HI}] / {list(LAMBDA_GRID)} "
f"CV folds {CV_FOLDS} (time-contiguous)")
say(f"map : m = round_{MAP_ROUND_TO}(clip(norm × (0.5 + clip((z−lo)/(hi−lo),0,1)"
f"), {MAP_CLIP_LO}, {MAP_CLIP_HI})) lo/hi = train q{MAP_Q_LO}/q{MAP_Q_HI}")
say("")
main_feats = resolve_features(df, [f.strip() for f in args.features.split(",") if f.strip()],
"--features")
say(f"main features ({len(main_feats)}): {', '.join(main_feats)}")
specs: List[Tuple[str, List[str], str]] = [(args.model, main_feats, args.model)]
for label, feats in (("isbuy", ISBUY_FEATURES), ("exhaustion", EXHAUSTION_FEATURES)):
have = [f for f in feats if f in df.columns]
if len(have) != len(feats):
say(f"NOTE: {label} baseline missing "
f"{', '.join(f for f in feats if f not in df.columns)} — using {len(have)}")
if have:
specs.append((label, have, "ridge"))
use_lgbm = (not args.skip_lgbm) and lgbm_available()
if not args.skip_lgbm and not use_lgbm:
say("lightgbm not installed — upper bound skipped")
if use_lgbm:
all_f = [c for c in df.columns if c.startswith("f_")]
specs.append(("lgbm", all_f, "lgbm"))
say(f"lightgbm upper bound: {len(all_f)} f_* columns")
say("")
results: Dict[str, Dict] = {}
table: List[Dict] = []
for label, feats, kind in specs:
say("-" * 96)
say(f"MODEL: {label} ({kind}, {len(feats)} feature(s))")
say("-" * 96)
res = walk_forward(df, feats, kind, args.bar_minutes, args.burn_in_months,
args.purge_bars, args.seed, verbose=True)
if res["skipped"]:
say(f" skipped quarters: {'; '.join(res['skipped'])}")
if res["oos"].empty:
say(" no OOS quarters — nothing to evaluate")
continue
ev = evaluate(res["oos"])
n_sh = args.n_shuffle if kind != "lgbm" else min(args.n_shuffle, LGBM_MAX_SHUFFLE)
if kind == "lgbm" and n_sh != args.n_shuffle:
say(f" NOTE (protocol deviation, flagged): lgbm shuffle null capped at "
f"{n_sh} repeats (never ships)")
sh = shuffle_null(df, feats, kind, args.bar_minutes, args.burn_in_months,
args.purge_bars, n_sh, args.seed) if n_sh > 0 else np.array([])
p95 = float(np.percentile(sh, 95)) if len(sh) else float("nan")
rank = float((sh < ev["lift"]).mean() * 100) if len(sh) else float("nan")
cl = acceptance(ev, ev["lift"], p95)
say(f" OOS quarters : {len(res['quarters'])} "
f"({res['quarters']['quarter'].iloc[0]} → {res['quarters']['quarter'].iloc[-1]}), "
f"{ev['n']} trades")
say(f" flat R : {ev['flat_r']:+.2f} PF {fmt(ev['flat_pf'])} "
f"maxDD {ev['flat_dd']:.2f}")
say(f" sized R : {ev['sized_r']:+.2f} PF {fmt(ev['sized_pf'])} "
f"maxDD {ev['sized_dd']:.2f}")
say(f" lift : {ev['lift']:+.2f} R")
say(f" halves (split {ev['median_time']}): "
f"h1 flat {ev['flat_h1']:+.2f} / sized {ev['sized_h1']:+.2f} (n={ev['n_h1']}) | "
f"h2 flat {ev['flat_h2']:+.2f} / sized {ev['sized_h2']:+.2f} (n={ev['n_h2']})")
say(f" multiplier : mean {ev['mean_m']:.3f} min {ev['min_m']:.2f} "
f"max {ev['max_m']:.2f}")
say(f" shuffle null : n={len(sh)} mean {fmt(float(sh.mean()) if len(sh) else None)} "
f"p95 {fmt(p95)} real-lift percentile rank {fmt(rank, 0)}%")
say(" acceptance:")
for k, v in cl.items():
if k == "overall":
continue
say(f" {k:<22} {'PASS' if v else 'FAIL'}")
say(f" {'OVERALL':<22} {'PASS' if cl['overall'] else 'FAIL'}")
say(" per-quarter:")
print_quarters(res["quarters"])
print_coef_signs(res["coefs"], label)
say("")
results[label] = {"res": res, "ev": ev, "shuffle": sh, "p95": p95,
"rank": rank, "clauses": cl, "features": feats, "kind": kind}
table.append({"model": label, "n": ev["n"], "flat_r": ev["flat_r"],
"sized_r": ev["sized_r"], "lift": ev["lift"],
"sized_pf": ev["sized_pf"], "sized_dd": ev["sized_dd"],
"mean_m": ev["mean_m"], "shuffle_p95": p95, "rank_pct": rank,
"pass": cl["overall"]})
if args.model not in results:
raise SystemExit("the main model produced no OOS quarters — check burn-in / "
"book window / MIN_TRAIN")
say("=" * 96)
say("COMPARISON")
say("=" * 96)
print_comparison(table)
say("")
say(f"OOS multiplier histogram ({args.model}):")
print_hist(results[args.model]["res"]["oos"]["mult"].to_numpy())
say("")
main_label = args.model
beats = None
if main_label in results and "isbuy" in results:
beats = results[main_label]["ev"]["lift"] > results["isbuy"]["ev"]["lift"]
say(f"beats isbuy baseline: {'yes' if beats else 'no'} "
f"({results[main_label]['ev']['lift']:+.2f} R vs "
f"{results['isbuy']['ev']['lift']:+.2f} R)")
say(f"OVERALL ACCEPTANCE ({main_label}): "
f"{'PASS' if results[main_label]['clauses']['overall'] else 'FAIL'}")
say("")
# ── artifact ──────────────────────────────────────────────────────
artifact = None
if args.write_artifact:
say("-" * 96)
say("FINAL FIT (all trades) + ARTIFACT")
say("-" * 96)
fit = final_fit(df, main_feats, args.model, args.seed)
say(f" lambda (CV on full window): {fit['lambda']:g}")
if np.isfinite(fit["grad_norm"]):
say(f" logistic |grad|_inf at convergence: {fit['grad_norm']:.3e}")
say(f" map lo/hi/normaliser: {fit['map']['lo']:.5f} / {fit['map']['hi']:.5f} / "
f"{fit['map']['normalizer']:.5f}")
artifact = build_artifact(
fit, main_feats, args.model, df, results[main_label]["ev"],
results[main_label]["clauses"], results[main_label]["ev"]["lift"],
results[main_label]["shuffle"], args.n_shuffle, args.book, args.tf,
[t.strip() for t in args.confirm_tfs.split(",") if t.strip()],
args.burn_in_months, args.purge_bars, args.spread)
# Self-check: the runtime formula must reproduce the in-sample map exactly.
rows = df[main_feats].to_dict("records")
repro = np.array([predict_multiplier(r, artifact) for r in rows])
exact = bool(np.array_equal(repro, fit["mult"]))
say(f" predict_multiplier reproduces in-sample multipliers exactly: "
f"{'YES' if exact else 'NO'} (max |diff| "
f"{float(np.max(np.abs(repro - fit['mult']))):.3e})")
if not exact:
raise SystemExit("predict_multiplier does not reproduce the final fit — abort")
args.write_artifact.parent.mkdir(parents=True, exist_ok=True)
args.write_artifact.write_text(json.dumps(artifact, indent=2, default=str))
say(f" wrote {args.write_artifact}")
say("")
# ── outputs ───────────────────────────────────────────────────────
if args.out:
out = args.out
out.mkdir(parents=True, exist_ok=True)
main_res = results[main_label]
main_res["res"]["oos"].to_csv(out / "oos_trades.csv", index=False)
main_res["res"]["quarters"].to_csv(out / "quarters.csv", index=False)
coefs = pd.concat(
[r["res"]["coefs"].assign(model=lbl) for lbl, r in results.items()
if not r["res"]["coefs"].empty],
ignore_index=True) if results else pd.DataFrame()
coefs.to_csv(out / "coefs_by_quarter.csv", index=False)
pd.DataFrame({
"model": np.repeat([lbl for lbl in results], [len(r["shuffle"]) for r in results.values()]),
"lift": np.concatenate([r["shuffle"] for r in results.values()])
if results else np.array([]),
}).to_csv(out / "shuffle_lifts.csv", index=False)
multiplier_hist(main_res["res"]["oos"]["mult"].to_numpy()).to_csv(
out / "multiplier_hist.csv", index=False)
payload = {
"book": str(args.book), "book_sha256": sha256_file(args.book),
"git_commit": git_commit(), "bar_minutes": args.bar_minutes,
"model": args.model, "features": main_feats, "seed": args.seed,
"n_shuffle": args.n_shuffle, "burn_in_months": args.burn_in_months,
"purge_bars": args.purge_bars,
"beats_isbuy": beats,
"models": {
lbl: {
"kind": r["kind"], "features": r["features"],
"evaluation": {k: (None if isinstance(v, float) and not np.isfinite(v)
else (str(v) if isinstance(v, pd.Timestamp) else v))
for k, v in r["ev"].items()},
"shuffle": {"n": int(len(r["shuffle"])),
"mean": float(r["shuffle"].mean()) if len(r["shuffle"]) else None,
"p95": None if not np.isfinite(r["p95"]) else r["p95"],
"real_lift_percentile": None if not np.isfinite(r["rank"])
else r["rank"]},
"acceptance": r["clauses"],
} for lbl, r in results.items()
},
"artifact": artifact,
}
(out / "results.json").write_text(json.dumps(payload, indent=2, default=str))
(out / "summary.txt").write_text("\n".join(_LINES) + "\n")
print(f"\nwrote outputs to {out}")
if __name__ == "__main__":
main()
+78
View File
@@ -0,0 +1,78 @@
"""Evaluate the violent-move WATCH rule (as shipped in violent_watch.py) bar by bar on history. Analysis tooling only.
Level per bar from: calendar (high-impact USD within 30 min ahead / 30 min after), clock slots, indicator readings.
Reports share of bars per level, burst rate next bar / next hour per level (both halves), and hot-hours per day.
python backtest/research/watch_rule_eval.py PARQUET BAR_MINUTES
"""
import sys
from pathlib import Path
import numpy as np, pandas as pd
ROOT = Path(__file__).resolve().parents[2]; sys.path.insert(0, str(ROOT)); sys.path.insert(0, str(ROOT / "backtest"))
from gpu_backtest import precompute_indicators
from config import TechnicalConfig
def wilder_atr(df, period=14):
pc = df["Close"].shift(1)
tr = pd.concat([df["High"] - df["Low"], (df["High"] - pc).abs(), (df["Low"] - pc).abs()], axis=1).max(axis=1)
return tr.ewm(alpha=1 / period, min_periods=period, adjust=False).mean()
def main():
pq, bm = sys.argv[1], int(sys.argv[2])
px = pd.read_parquet(pq); px = px[(px.index >= pd.Timestamp("2023-08-05", tz="UTC")) & (px.index < pd.Timestamp("2026-09-05", tz="UTC"))]
ic = precompute_indicators(px, TechnicalConfig())
sh = lambda a: pd.Series(np.asarray(a, dtype=float), index=px.index).shift(1)
rsi, adxs, hist, bbp, kcu, kcl = sh(ic.rsi), sh(ic.adx_slope), sh(ic.macd_hist), sh(ic.bb_position), sh(ic.kc_upper), sh(ic.kc_lower)
ef, em, es, et, adx, atr1 = sh(ic.ema_fast), sh(ic.ema_medium), sh(ic.ema_slow), sh(ic.ema_trend), sh(ic.adx), sh(ic.atr)
c1 = px["Close"].shift(1); hp = hist.shift(1)
atr_prev = wilder_atr(px).shift(1); rng = px["High"] - px["Low"]
burst = (rng >= 3.0 * atr_prev).astype(int); H = max(1, 60 // bm)
next_hour = burst[::-1].rolling(H, min_periods=1).max()[::-1].astype(int)
align = np.sign(c1 - ef) + np.sign(ef - em) + np.sign(em - es) + np.sign(es - et)
regime = np.select([(adx >= 25) & (align > 0), (adx >= 25) & (align < 0), adx >= 25, adxs > 3.0], ["trending_up", "trending_down", "trending", "transitioning"], default="ranging")
atr_pct = atr1.rolling(50).apply(lambda w: (w[-1] > w).mean() * 100, raw=True); extreme = atr_pct > 80
# readings
bear = ((bbp < 0).astype(int) + (rsi < 30).astype(int) + (c1 < kcl).astype(int) + ((hist < 0) & (hist < hp) & (hp < 0)).astype(int))
bull = ((bbp > 1).astype(int) + (rsi >= 70).astype(int) + (c1 > kcu).astype(int) + ((hist > 0) & (hist > hp) & (hp > 0)).astype(int))
neutral = (adxs > 3).astype(int) + extreme.astype(int) + (regime == "transitioning").astype(int)
readings = bear + bull + neutral
# clock (now = bar close time)
ct = px.index + pd.Timedelta(minutes=bm); wd = ct.weekday; hr = ct.hour
clock = ((wd <= 4) & np.isin(hr, [12, 13, 14])) | (hr == 1) | ((wd == 6) & np.isin(hr, [22, 23]))
# calendar
cal = pd.read_parquet(ROOT / "backtest/cache/calendar.parquet"); cal["ts"] = pd.to_datetime(cal.dateline, unit="s", utc=True)
hi = np.sort(cal[(cal.currency == "USD") & (cal.impactName.str.lower() == "high")].ts.values.astype("datetime64[ns]"))
t = ct.values.astype("datetime64[ns]"); nx = np.searchsorted(hi, t)
to_next = np.where(nx < len(hi), (hi[np.minimum(nx, len(hi) - 1)] - t) / np.timedelta64(1, "m"), 1e9)
since = np.where(nx - 1 >= 0, (t - hi[np.maximum(nx - 1, 0)]) / np.timedelta64(1, "m"), 1e9)
calendar = ((to_next > 0) & (to_next <= 30)) | (since <= 30)
# FINAL RULE (2026-09-06, chosen from the variant sweep recorded in sizing_campaign_2026-09.md):
# extreme = price-at-extreme readings only (outside Bollinger, RSI beyond 30/70, outside Keltner)
# hot = calendar | (clock & extreme >= 1) -> ~8.6% of bars, next-bar odds 3-5x normal, ~4 hot hours/day
# watch = clock | extreme >= 2 -> ~2x normal
extreme = ((bbp < 0) | (bbp > 1)).astype(int) + ((rsi < 30) | (rsi >= 70)).astype(int) + ((c1 < kcl) | (c1 > kcu)).astype(int)
level = np.where(calendar | (clock & (extreme >= 1)), "hot", np.where(clock | (extreme >= 2), "watch", "none"))
d = pd.DataFrame({"level": level, "burst": burst.values, "next_hour": next_hour.values, "readings": readings.values, "clock": clock, "calendar": calendar}, index=px.index)
d = d[d.index >= pd.Timestamp("2023-09-05", tz="UTC")]
half = np.where(np.arange(len(d)) < len(d) // 2, 1, 2)
b1, b2 = d.burst[half == 1].mean(), d.burst[half == 2].mean(); h1, h2 = d.next_hour[half == 1].mean(), d.next_hour[half == 2].mean()
print(f"{Path(pq).name}: {len(d)} bars; base next-bar {100*d.burst.mean():.2f}%, next-hour {100*d.next_hour.mean():.1f}%")
rows = []
for lv in ("none", "watch", "hot"):
m = d.level == lv
rows.append(dict(level=lv, share_of_bars=100 * m.mean(), next_bar_pct=100 * d.burst[m].mean(), lift_h1=d.burst[m & (half == 1)].mean() / b1, lift_h2=d.burst[m & (half == 2)].mean() / b2,
next_hour_pct=100 * d.next_hour[m].mean(), hour_lift_h1=d.next_hour[m & (half == 1)].mean() / h1, hour_lift_h2=d.next_hour[m & (half == 2)].mean() / h2,
share_of_all_bursts=100 * d.burst[m].sum() / d.burst.sum()))
print(pd.DataFrame(rows).round(2).to_string(index=False))
hot = d[d.level == "hot"]
kinds = pd.Series(np.where(hot.calendar, "calendar", "clock+extreme reading"), index=hot.index)
print("hot bars by cause:", (100 * kinds.value_counts(normalize=True)).round(1).to_dict())
hot_hours = hot.index.floor("h").unique()
days = pd.Series(hot_hours.date).value_counts()
print(f"distinct HOT hours per calendar day: mean {days.mean():.1f}, median {days.median():.0f}, max {days.max()} (≈ heads-up messages/day before the daily cap)")
hb = d.burst.groupby(d.index.floor("h")).max(); hh = pd.Series(1, index=hot_hours)
print(f"of all hours with a violent bar, the hour was flagged HOT at its start or during: {100*hb[hb==1].index.isin(hot_hours).mean():.1f}%")
if __name__ == "__main__":
main()