#!/usr/bin/env python3 """Is gold sold down into the London PM fix? The whole run. The question, the windows, the statistic, the null, the day floor and the kill rule are in `rule.py`, sealed before any return around the fix was computed. Nothing here decides anything rule.py did not fix: the constants and the three decision functions (`verdict`, `critical_k`, `certified_floor`) are imported from it rather than restated, so a reader can diff this file against the seal and see that the rule was not quietly edited to fit the answer. python analyse.py [--bars PATH] [--out DIR] [--seed N] [--draws N] Writes results.md, results.json and days.csv into --out (this directory by default). results.md has no author but this script; it is never hand-edited. Decisions taken here that rule.py left to the implementation (there is no one to ask, so they are written down): 1. The placebo pool is the 496 start minutes 08:00..16:15 London inclusive; a start is valid when a bar exists at S and at S+15 on that London day. The day floor is "both the 14:45 and the 15:00 bar exist" AND "at least 400 of the 496 starts are valid", exactly as sealed. 2. The AM-fix descriptive (10:15 -> 10:30) is run on the same kept days and the same placebo table as the PM test, with its own fresh 10,000 draws. Days that clear the PM floor but lack an AM bar are dropped from the AM figure only, and the count is printed. Re-deriving a separate day floor for a descriptive would invent a rule the seal does not contain. 3. The shift-injection control uses 2,000 draws per control run rather than 10,000: 6 deltas x 300 runs is 1,800 nulls, and the seal's own text allows the reduction as long as it is stated. It is stated in results.md and in the JSON. The main test keeps the full 10,000 draws. 4. The null's own size at the critical count k01 is computed with rule.py's own binomial tail (`rule._binom_tail`) rather than a second implementation, for the same reason as (the rest of) the imports. It is private only by naming. """ import argparse import json import os import pathlib import sys import time import numpy as np import pandas as pd sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import rule # noqa: E402 from rule import certified_floor, critical_k, verdict # noqa: E402 BARS = "~/data/bars/XAUUSD_mt5_m1.csv" LONDON = "Europe/London" # minute-of-London-day indices M_1445, M_1500, M_1510, M_1525 = 14 * 60 + 45, 15 * 60, 15 * 60 + 10, 15 * 60 + 25 M_1015, M_1030 = 10 * 60 + 15, 10 * 60 + 30 POOL_LO, POOL_HI = 8 * 60, 16 * 60 + 15 # 08:00 .. 16:15 inclusive WINDOW = 15 # minutes POOL_N = POOL_HI - POOL_LO + 1 # 496 POOL_FLOOR = 400 WEEKDAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"] # ----------------------------------------------------------------- the data def load(path): """One-minute bars on the Europe/London clock, mid of the closes. Returns the frame and a (days x 1440) matrix of mid prices, NaN where the broker quoted no bar, so that every window below is an array slice. """ df = pd.read_csv(path, usecols=["minute", "bid_close", "ask_close"], parse_dates=["minute"]) ts = df["minute"].dt.tz_convert(LONDON) df["mid"] = (df["bid_close"] + df["ask_close"]) / 2.0 df["day"] = ts.dt.date df["m"] = ts.dt.hour * 60 + ts.dt.minute days = np.array(sorted(df["day"].unique())) day_ix = {d: i for i, d in enumerate(days)} mid = np.full((len(days), 1440), np.nan) mid[df["day"].map(day_ix).to_numpy(), df["m"].to_numpy()] = df["mid"].to_numpy() return df, days, mid def bp(a, b): """ln(b/a) x 10,000, NaN-safe.""" with np.errstate(invalid="ignore", divide="ignore"): return np.log(b / a) * 1e4 # ------------------------------------------------------------ placebo table def placebo_table(mid): """Per day, the 15-minute log returns (bp) of every valid pool start. Hands back a padded (days x 496) matrix of the *valid* returns packed to the left, and the per-day count of valid starts, which is what a draw indexes into. """ starts = mid[:, POOL_LO:POOL_HI + 1] ends = mid[:, POOL_LO + WINDOW:POOL_HI + WINDOW + 1] rets = bp(starts, ends) valid = np.isfinite(rets) sizes = valid.sum(axis=1) # pack each row's valid returns to the left without a Python loop over days order = np.argsort(~valid, axis=1, kind="stable") packed = np.take_along_axis(rets, order, axis=1) cols = np.arange(rets.shape[1])[None, :] packed[cols >= sizes[:, None]] = np.nan return packed, sizes def draw(packed, sizes, draws, rng): """`draws` placebo draws: one uniformly chosen valid start per day. Returns the (draws x days) matrix of chosen returns. Vectorised: an index matrix is drawn once, scaled per day by that day's own pool size. """ u = rng.random((draws, len(sizes))) idx = (u * sizes[None, :]).astype(np.int64) return packed[np.arange(len(sizes))[None, :], idx] def null_stats(packed, sizes, draws, rng): """Fraction-negative and mean statistics over `draws` placebo draws.""" r = draw(packed, sizes, draws, rng) neg = (r < 0).sum(axis=1) return neg / r.shape[1], neg, r.mean(axis=1) # -------------------------------------------------------------- the per-day def build_days(days, mid, packed, sizes): """One row per London calendar day, with the floor decision and its reason.""" m1445, m1500 = mid[:, M_1445], mid[:, M_1500] rows = pd.DataFrame({ "date": [str(d) for d in days], "weekday": [WEEKDAYS[d.weekday()] for d in days], "mid_1445": m1445, "mid_1500": m1500, "run_in_bp": bp(m1445, m1500), "pool_size": sizes, "am_run_in_bp": bp(mid[:, M_1015], mid[:, M_1030]), "auction_bp": bp(m1500, mid[:, M_1510]), "after_bp": bp(mid[:, M_1510], mid[:, M_1525]), }) has_bars = np.isfinite(m1445) & np.isfinite(m1500) clears = sizes >= POOL_FLOOR rows["kept"] = has_bars & clears reason = np.where(rows["kept"], "kept", np.where(~has_bars & ~clears, "no 14:45/15:00 bar and pool under 400", np.where(~has_bars, "no 14:45/15:00 bar", "placebo pool under 400 of 496"))) rows["drop_reason"] = reason return rows # ------------------------------------------------------------------ the test def sign_test(obs, packed, sizes, draws, rng): """F, C, p for a per-day return vector against the paired placebo null.""" f = float((obs < 0).mean()) frac, neg, means = null_stats(packed, sizes, draws, rng) c = float(frac.mean()) p = float((frac >= f).mean()) return f, c, p, frac, neg, means def injection(mu, packed, sizes, runs, draws, deltas, seed): """Shift injection for the descriptive mean: subtract DELTA bp from every day. Each run gets a fresh null of its own (the null does not depend on the data, so a fresh null is exactly the sampling variability the control is meant to expose) and a fresh seed derived from the run's ordinal. """ rows = [] for d in deltas: shifted = mu - d caught = 0 for i in range(runs): rng = np.random.default_rng([seed, int(d), i]) _, _, means = null_stats(packed, sizes, draws, rng) if float((means <= shifted).mean()) < rule.P_SUPPORT: caught += 1 rows.append({"delta_bp": int(d), "runs": runs, "draws": draws, "shifted_mu_bp": float(shifted), "certified": caught / runs}) floor = next((r["delta_bp"] for r in rows if r["certified"] >= 0.95), None) return rows, floor # -------------------------------------------------------------------- output def main(argv=None): t0 = time.time() ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) ap.add_argument("--bars", default=BARS) ap.add_argument("--out", default=str(pathlib.Path(__file__).parent)) ap.add_argument("--seed", type=int, default=20260912) ap.add_argument("--draws", type=int, default=10_000) ap.add_argument("--control-draws", type=int, default=2_000) args = ap.parse_args(argv) out = pathlib.Path(args.out) df, days, mid = load(pathlib.Path(args.bars).expanduser()) print(f"bars: {len(df)} rows, {df['minute'].min()} .. {df['minute'].max()}") packed, sizes = placebo_table(mid) rows = build_days(days, mid, packed, sizes) rows.to_csv(out / "days.csv", index=False, columns=["date", "mid_1445", "mid_1500", "run_in_bp", "pool_size", "am_run_in_bp", "auction_bp", "after_bp", "kept"]) kept = rows[rows["kept"]].reset_index(drop=True) k_ix = np.flatnonzero(rows["kept"].to_numpy()) kpacked, ksizes = packed[k_ix], sizes[k_ix] n = len(kept) drops = (rows.loc[~rows["kept"], "drop_reason"].value_counts().to_dict()) result = { "bars": int(len(df)), "first_bar": str(df["minute"].min()), "last_bar": str(df["minute"].max()), "days_seen": int(len(rows)), "n": n, "drops": drops, "seed": args.seed, "draws": args.draws, "floor": {"pool_min": POOL_FLOOR, "pool_size": POOL_N, "min_days": rule.MIN_DAYS}, } run_in = kept["run_in_bp"].to_numpy() rng = np.random.default_rng(args.seed) f, c, p, frac, neg, means = sign_test(run_in, kpacked, ksizes, args.draws, rng) k01 = critical_k(n, c, rule.P_SUPPORT) q_star = certified_floor(n, c, k01) size_exact = float(rule._binom_tail(n, c, k01)) # the rule's own tail size_empirical = float((neg >= k01).mean()) v = verdict(f, c, p, q_star, n) result["test"] = { "F": f, "down_days": int((run_in < 0).sum()), "C": c, "p": p, "k01": k01, "q_star": q_star, "null_size_at_k01_exact": size_exact, "null_size_at_k01_empirical": size_empirical, "null_sd": float(frac.std()), } result["verdict"] = v result["published_per_day"] = n >= rule.MIN_DAYS # --- descriptives ------------------------------------------------- mu = float(run_in.mean()) p_mu = float((means <= mu).mean()) desc = {"mu_bp": mu, "sd_bp": float(run_in.std(ddof=1)), "p_mu": p_mu, "null_mean_centre_bp": float(means.mean())} am = kept["am_run_in_bp"].to_numpy() am_ok = np.isfinite(am) rng_am = np.random.default_rng(args.seed + 1) af, ac, ap_, _, _, am_means = sign_test(am[am_ok], kpacked[am_ok], ksizes[am_ok], args.draws, rng_am) desc["am_fix"] = {"days": int(am_ok.sum()), "dropped_no_bar": int((~am_ok).sum()), "F": af, "C": ac, "p": ap_, "mu_bp": float(am[am_ok].mean()), "p_mu": float((am_means <= am[am_ok].mean()).mean())} for name, col in (("auction", "auction_bp"), ("after", "after_bp")): v_ = kept[col].to_numpy() ok = np.isfinite(v_) desc[name] = {"days": int(ok.sum()), "mean_bp": float(v_[ok].mean()), "fraction_negative": float((v_[ok] < 0).mean())} pool_abs = np.abs(kpacked[np.isfinite(kpacked)]) desc["median_abs_run_in_bp"] = float(np.median(np.abs(run_in))) desc["median_abs_placebo_bp"] = float(np.median(pool_abs)) desc["placebo_windows_pooled"] = int(pool_abs.size) wk = [] for d in WEEKDAYS: sel = kept[kept["weekday"] == d] if len(sel): wk.append({"weekday": d, "days": int(len(sel)), "F": float((sel["run_in_bp"] < 0).mean()), "mean_bp": float(sel["run_in_bp"].mean())}) desc["weekday"] = wk result["descriptive"] = desc # --- the positive control ---------------------------------------- ctrl, ctrl_floor = injection(mu, kpacked, ksizes, rule.CONTROL_RUNS, args.control_draws, rule.DELTAS, args.seed) result["injection"] = {"rows": ctrl, "smallest_certified_delta_bp": ctrl_floor, "draws_per_run": args.control_draws, "note": ("2,000 draws per control run rather than 10,000: " "1,800 nulls at full size was too slow. The main " "test keeps 10,000.")} result["runtime_seconds"] = round(time.time() - t0, 1) (out / "results.json").write_text(json.dumps(result, indent=1) + "\n") write_report(out, result) print(f"verdict: {v} — F={f:.4f} C={c:.4f} p={p:.4f} n={n} " f"k01={k01} q*={q_star} size@k01={size_exact:.4f} " f"(empirical {size_empirical:.4f})") print(f"MU={mu:+.3f} bp sd={desc['sd_bp']:.3f} p_mu={p_mu:.4f}") print(f"AM fix: F={af:.4f} C={ac:.4f} p={ap_:.4f} on {int(am_ok.sum())} days") print(f"auction 15:00->15:10 mean={desc['auction']['mean_bp']:+.3f} bp " f"neg={desc['auction']['fraction_negative']:.4f}; " f"after 15:10->15:25 mean={desc['after']['mean_bp']:+.3f} bp " f"neg={desc['after']['fraction_negative']:.4f}") print(f"median |run_in| = {desc['median_abs_run_in_bp']:.3f} bp vs " f"median |placebo| = {desc['median_abs_placebo_bp']:.3f} bp") for r in ctrl: print(f" delta {r['delta_bp']:>2} bp: certified {r['certified'] * 100:.1f} %") print(f"smallest certified delta: {ctrl_floor}") for r in wk: print(f" {r['weekday']:<9} n={r['days']:>3} F={r['F']:.4f} " f"mean={r['mean_bp']:+.3f} bp") for k, vcount in drops.items(): print(f" dropped {vcount}: {k}") print(f"runtime: {result['runtime_seconds']} s") return 0 def write_report(out, r): t, d = r["test"], r["descriptive"] n = r["n"] lines = [ "# Is gold sold down into the London PM fix?", "", f"**Verdict: {r['verdict']}** — F = {t['F']:.4f} against a placebo centre " f"C = {t['C']:.4f}, one-sided p = {t['p']:.4f}, on {n} London days.", "", f"Bars: {r['bars']} one-minute rows, {r['first_bar']} to {r['last_bar']} " f"(UTC as stored; every window below is on the Europe/London clock). " f"Price is the mid of the bid and ask closes.", "", ] if not r["published_per_day"]: lines += [f"Fewer than {rule.MIN_DAYS} days survived the floor " f"({n}), so **no per-day figure is published**, as the rule " "requires.", ""] lines += [ "## The test", "", "| | |", "|---|---:|", f"| days kept, n | {n} (of {r['days_seen']} London days seen) |", f"| down days (run_in < 0) | {t['down_days']} |", f"| F, fraction negative | {t['F']:.4f} |", f"| C, placebo centre | {t['C']:.4f} |", f"| null sd of F | {t['null_sd']:.4f} |", f"| one-sided p | {t['p']:.4f} |", f"| critical count k01 (p < 0.01 at C) | {t['k01']} |", f"| certified floor q\\* | {t['q_star']} |", f"| null's own size at k01, exact binomial | {t['null_size_at_k01_exact']:.4f} |", f"| null's own size at k01, from the draws | {t['null_size_at_k01_empirical']:.4f} |", "", f"The null is the pre-registered placebo-window null: {r['draws']} draws, " f"each picking one start minute per day uniformly from the " f"{r['floor']['pool_size']} minutes 08:00-16:15 London that have a bar " f"exactly 15 minutes later. Seed {r['seed']}.", "", "## Days dropped by the floor", "", "| reason | days |", "|---|---:|", ] for k, v in r["drops"].items(): lines.append(f"| {k} | {v} |") if not r["drops"]: lines.append("| (none) | 0 |") lines += [ "", "## Descriptive only — none of this can move the verdict", "", f"**The mean.** MU = {d['mu_bp']:+.4f} bp, per-day sd " f"{d['sd_bp']:.4f} bp, placebo p = {d['p_mu']:.4f} (the fraction of draws " f"whose mean is at or below MU; the null's own mean centre is " f"{d['null_mean_centre_bp']:+.4f} bp).", "", f"**The AM fix, 10:15 to 10:30 London.** F = {d['am_fix']['F']:.4f}, " f"C = {d['am_fix']['C']:.4f}, p = {d['am_fix']['p']:.4f}, mean " f"{d['am_fix']['mu_bp']:+.4f} bp (p = {d['am_fix']['p_mu']:.4f}) on " f"{d['am_fix']['days']} days; {d['am_fix']['dropped_no_bar']} kept days had " "no 10:15 or 10:30 bar. Its own draws, not the PM test's.", "", f"**The auction itself, 15:00 to 15:10.** Mean " f"{d['auction']['mean_bp']:+.4f} bp, negative on " f"{d['auction']['fraction_negative'] * 100:.1f} % of " f"{d['auction']['days']} days.", "", f"**After it, 15:10 to 15:25.** Mean {d['after']['mean_bp']:+.4f} bp, " f"negative on {d['after']['fraction_negative'] * 100:.1f} % of " f"{d['after']['days']} days.", "", f"**Busy, not just lower.** Median |run_in| = " f"{d['median_abs_run_in_bp']:.4f} bp against a median |placebo return| of " f"{d['median_abs_placebo_bp']:.4f} bp, pooled over all " f"{d['placebo_windows_pooled']} valid placebo windows.", "", "**By weekday** (descriptive; the rule fixes no weekday claim).", "", "| weekday | days | F | mean run-in |", "|---|---:|---:|---:|", ] for w in d["weekday"]: lines.append(f"| {w['weekday']} | {w['days']} | {w['F']:.4f} | " f"{w['mean_bp']:+.3f} bp |") inj = r["injection"] lines += [ "", "## The positive control: shift injection for the mean", "", f"DELTA bp subtracted from every day's run-in, {rule.CONTROL_RUNS} runs a " f"level, each against a fresh null, counting the runs that reach " f"p < {rule.P_SUPPORT} on the mean. {inj['note']}", "", "| DELTA | shifted mean | runs at p < 0.01 |", "|---:|---:|---:|", ] for row in inj["rows"]: lines.append(f"| {row['delta_bp']} bp | {row['shifted_mu_bp']:+.3f} bp | " f"{row['certified'] * 100:.1f} % |") floor = inj["smallest_certified_delta_bp"] lines += [ "", (f"Smallest shift this sample certifies at 95 % or better: " f"**{floor} bp**." if floor is not None else "No injected shift reached 95 % certification, which is itself the " "finding: this sample could not have seen a mean move of even " f"{max(rule.DELTAS)} bp."), "", "The column reads 0 % or 100 % and little between because the injection " "is deterministic: subtracting DELTA moves this one sample's mean by " "exactly DELTA, so the only thing varying across the " f"{rule.CONTROL_RUNS} runs at a level is the null's own Monte Carlo " "noise. The table is therefore a statement about where the p < 0.01 " "threshold sits for this sample, not a power curve over resampled data. " "That is what the sealed rule asked for and it is reported as it came.", "", # runtime deliberately not printed here: results.md must be byte-identical # on a re-run from a clean clone, and a timing would break that. "Method, windows and kill rule: `rule.py`, sealed before any return " "around the fix was computed. This file is written by `analyse.py` and " "never by hand.", "", "This is a test record, not advice.", ] (out / "results.md").write_text("\n".join(lines) + "\n") if __name__ == "__main__": sys.exit(main())