# Published copy. The two paths below are set to environment variables so the # script runs against any clean one-minute dataset with columns t, o, h, l. # The original reads them from the lab's own store. # -*- coding: utf-8 -*- """ Study 03: what a trade suffers from an entry, and which barrier it reaches first. Study 01 measured how far a five-minute candle travels and said plainly that this is not the chance a stop gets hit. This measures the thing itself. Rewritten 20 Sep 2026 after a methodology review raised five points, all of which were right: 1. Horizons were truncated at the close. An entry at 15:50 was given ten minutes and counted as though it had thirty. Every horizon now requires the whole window to sit inside the session; entries that do not qualify are dropped from that horizon and the count is published. 2. The first-touch race ran long only. It now runs both sides. 3. Regime terciles are cut across the whole sample. The label uses only prior-day information, but the boundaries are ex post, and the output says so instead of implying the cut was knowable at the time. 4. Minute entries overlap almost entirely, so 199,290 of them are nothing like 199,290 independent trials. Uncertainty is estimated by resampling whole sessions, which is the unit closer to independent. 5. A stop nearer than the target is reached first more often in any symmetric walk. That is geometry, not evidence about the stop, so the split conditional on a decision and the geometric expectation are published beside the raw figures. Read-only. Writes research/mae_v1.json. """ import glob, os, json import numpy as np import pandas as pd STORE = os.environ.get("BARS_DIR", "./bars") # one-minute CSVs: t,o,h,l OUT = os.path.dirname(os.path.abspath(__file__)) TZ = "America/New_York" CUTOFF = pd.Timestamp("2026-09-17", tz=TZ) + pd.Timedelta(days=1) OPEN_M, CLOSE_M = 570, 960 HORIZONS = [1, 5, 10, 15, 30] PCTS = [50, 75, 90, 95] LEVELS = {"NQ": [4, 6, 8, 10, 15, 20, 25, 30], "ES": [1, 2, 3, 4, 5, 6, 8]} PAIRS_BY_SYMBOL = {"NQ": [(10, 10), (10, 20), (10, 30), (15, 30), (20, 40), (25, 25)], "ES": [(2, 2), (2, 4), (2, 6), (3, 6), (4, 8), (5, 5)]} FT_HORIZON = 30 BOOT = 400 RNG = np.random.default_rng(20260920) def bars(symbol): files = sorted(glob.glob(os.path.join(STORE, f"{symbol}_FLOW_*.csv"))) raw = pd.concat([pd.read_csv(f, usecols=["t", "o", "h", "l"]) for f in files]) df = raw.drop_duplicates("t").sort_values("t") df["ts"] = pd.to_datetime(df.t, unit="ms", utc=True).dt.tz_convert(TZ) df = df.set_index("ts") return df[df.index < CUTOFF] def sessions(df): """Complete cash sessions only: all 390 minutes present, as in Study 02.""" m = df.index.hour * 60 + df.index.minute rth = df[(m >= OPEN_M) & (m < CLOSE_M)] grid = list(range(OPEN_M, CLOSE_M)) out = [] for day, g in rth.groupby(rth.index.date): g = g.sort_index() if sorted((g.index.hour * 60 + g.index.minute).tolist()) == grid: out.append((day, g)) return out def per_session(g, pairs): """Excursions and first-touch outcomes for one session. Full horizons only.""" o = g.o.values.astype(float) hi = g.h.values.astype(float) lo = g.l.values.astype(float) n = len(o) rec = {"minute": np.arange(n)} for H in HORIZONS: mae_l = np.full(n, np.nan) mae_s = np.full(n, np.nan) mfe_l = np.full(n, np.nan) for i in range(n - H + 1): # the whole horizon must fit j = i + H worst_low = lo[i:j].min() best_high = hi[i:j].max() mae_l[i] = o[i] - worst_low mae_s[i] = best_high - o[i] mfe_l[i] = best_high - o[i] rec["long_%d" % H] = mae_l rec["short_%d" % H] = mae_s rec["longmfe_%d" % H] = mfe_l for stop, target in pairs: for side in ("long", "short"): res = np.full(n, "", dtype=object) for i in range(n - FT_HORIZON + 1): entry = o[i] outcome = "neither" for k in range(i, i + FT_HORIZON): if side == "long": s_hit = (entry - lo[k]) >= stop t_hit = (hi[k] - entry) >= target else: s_hit = (hi[k] - entry) >= stop t_hit = (entry - lo[k]) >= target if s_hit and t_hit: outcome = "ambiguous"; break if s_hit: outcome = "stop"; break if t_hit: outcome = "target"; break res[i] = outcome rec["ft_%s_%d_%d" % (side, stop, target)] = res return pd.DataFrame(rec) def regime_map(sess): """Yesterday's range decides the label. The tercile cuts are ex post.""" days = [d for d, _ in sess] rng = {d: float(g.h.max() - g.l.min()) for d, g in sess} prev = {days[i]: rng[days[i - 1]] for i in range(1, len(days))} lo_c, hi_c = np.percentile(np.array(list(prev.values())), [33.333, 66.667]) lab = {str(d): ("quiet" if v < lo_c else ("busy" if v > hi_c else "ordinary")) for d, v in prev.items()} return lab, round(float(lo_c), 2), round(float(hi_c), 2) def describe(s, levels): s = s.dropna().values if len(s) == 0: return None return dict(n=int(len(s)), pct={str(p): round(float(np.percentile(s, p)), 2) for p in PCTS}, mean=round(float(s.mean()), 2), reached={str(x): round(float((s >= x).mean()) * 100, 2) for x in levels}) def boot_ci(by_session, stat): """Resample whole sessions; minute entries inside one are not independent.""" keys = [k for k, v in by_session.items() if len(v)] if not keys: return None vals = [] for _ in range(BOOT): pick = RNG.choice(len(keys), size=len(keys), replace=True) pooled = np.concatenate([by_session[keys[i]] for i in pick]) vals.append(stat(pooled)) return [round(float(np.percentile(vals, 2.5)), 2), round(float(np.percentile(vals, 97.5)), 2)] def build(symbol): pairs = PAIRS_BY_SYMBOL[symbol] levels = LEVELS[symbol] sess = sessions(bars(symbol)) frames = {} for day, g in sess: f = per_session(g, pairs) f["day"] = str(day) frames[str(day)] = f ex = pd.concat(frames.values(), ignore_index=True) lab, lo_c, hi_c = regime_map(sess) ex["regime"] = ex.day.map(lab) ex["half"] = (ex.minute // 30) * 30 ftcol0 = "ft_long_%d_%d" % pairs[0] out = dict( sessions=len(sess), entries=int(len(ex)), first=str(sess[0][0]), last=str(sess[-1][0]), full_horizon_entries={str(H): int(ex["long_%d" % H].notna().sum()) for H in HORIZONS}, first_touch_entries=int((ex[ftcol0] != "").sum()), regime_cuts=dict(quiet_below=lo_c, busy_above=hi_c, note="terciles of the previous session's range. The label uses only " "prior-day information; the tercile boundaries are computed across " "the whole sample and are therefore ex post."), overall={}, overall_mfe={}, by_half_hour={}, by_regime={}, first_touch={}, uncertainty={}) for side in ("long", "short"): out["overall"][side] = {str(H): describe(ex["%s_%d" % (side, H)], levels) for H in HORIZONS} out["overall_mfe"]["long"] = {str(H): describe(ex["longmfe_%d" % H], levels) for H in HORIZONS} for h, g in ex.groupby("half"): key = "%02d:%02d" % ((h + OPEN_M) // 60, (h + OPEN_M) % 60) out["by_half_hour"][key] = {side: {str(H): describe(g["%s_%d" % (side, H)], levels) for H in (5, 15)} for side in ("long", "short")} for r, g in ex.groupby("regime"): out["by_regime"][r] = {side: {str(H): describe(g["%s_%d" % (side, H)], levels) for H in (5, 15)} for side in ("long", "short")} for stop, target in pairs: for side in ("long", "short"): col = "ft_%s_%d_%d" % (side, stop, target) v = ex[col][ex[col] != ""] if not len(v): continue c = v.value_counts(normalize=True) * 100 dec = v[v.isin(["stop", "target"])] dc = dec.value_counts(normalize=True) * 100 if len(dec) else None out["first_touch"]["%s_%d_%d" % (side, stop, target)] = dict( n=int(len(v)), target_first=round(float(c.get("target", 0.0)), 2), stop_first=round(float(c.get("stop", 0.0)), 2), ambiguous=round(float(c.get("ambiguous", 0.0)), 2), neither=round(float(c.get("neither", 0.0)), 2), decided_n=int(len(dec)), decided_target_first=round(float(dc.get("target", 0.0)), 2) if dc is not None else None, decided_stop_first=round(float(dc.get("stop", 0.0)), 2) if dc is not None else None, geometry_stop_first=round(target / float(stop + target) * 100, 2)) lvl = 10 if symbol == "NQ" else 2 d1 = {k: f["long_5"].dropna().values for k, f in frames.items()} out["uncertainty"]["long_5_reached_%d" % lvl] = dict( point=round(float((ex["long_5"].dropna().values >= lvl).mean() * 100), 2), ci95_session_bootstrap=boot_ci(d1, lambda a: float((a >= lvl).mean() * 100)), resamples=BOOT, unit="whole sessions") sk = pairs[1] col = "ft_long_%d_%d" % sk d2 = {k: f[col][f[col] != ""].values for k, f in frames.items()} out["uncertainty"][col + "_stop_first"] = dict( point=out["first_touch"]["long_%d_%d" % sk]["stop_first"], ci95_session_bootstrap=boot_ci(d2, lambda a: float((a == "stop").mean() * 100)), resamples=BOOT, unit="whole sessions") return out def main(): rep = dict(generated="2026-09-20", cutoff="2026-09-17", session="09:30 to 16:00 America/New_York, complete sessions only", horizon_rule="an entry is used for a horizon only if the whole horizon fits " "inside the same session, so no window is silently shortened", independence="minute entries overlap, so uncertainty is estimated by resampling " "whole sessions rather than treating each entry as independent", definition="MAE is the worst move against the position within the horizon, " "measured from the entry price, entering at the open of a minute", symbols={}) for sym in ("NQ", "ES"): rep["symbols"][sym] = build(sym) d = rep["symbols"][sym] print("\n=== %s %d sessions, %s minute rows" % (sym, d["sessions"], format(d["entries"], ","))) print(" full-horizon entries:", d["full_horizon_entries"]) print(" first-touch entries :", format(d["first_touch_entries"], ",")) for side in ("long", "short"): for H in (5, 15): s = d["overall"][side][str(H)] print(" %-5s %2dm n=%-7d median %6.2f p90 %6.2f" % (side, H, s["n"], s["pct"]["50"], s["pct"]["90"])) print(" first touch: target / stop / ambiguous / neither | decided t vs s | geometry s") for k, v in d["first_touch"].items(): print(" %-18s %5.2f %5.2f %5.2f %5.2f | %5.2f %5.2f | %5.2f" % (k, v["target_first"], v["stop_first"], v["ambiguous"], v["neither"], v["decided_target_first"], v["decided_stop_first"], v["geometry_stop_first"])) for k, v in d["uncertainty"].items(): print(" %-26s %.2f%% 95%% %s" % (k, v["point"], v["ci95_session_bootstrap"])) path = os.path.join(OUT, "mae_v1.json") with open(path, "w", encoding="utf-8") as f: json.dump(rep, f, indent=2) print("\nwritten", path) if __name__ == "__main__": main()