#!/usr/bin/env python3
"""Half-hour profile of gold movement and spread, on the New York clock.

Usage:  ~/.venv/bin/python profile.py [bars.csv]
Reads ~/data/bars/XAUUSD_mt5_m1.csv by default (built by research/dukascopy/bars_mt5.py).
Writes bins.csv, weekdays.csv and results.md next to this script. Definitions: README.md.
"""
import os
import sys

import numpy as np
import pandas as pd

HERE = os.path.dirname(os.path.abspath(__file__))
BARS = sys.argv[1] if len(sys.argv) > 1 else os.path.expanduser("~/data/bars/XAUUSD_mt5_m1.csv")
TZ = "America/New_York"

df = pd.read_csv(BARS, parse_dates=["minute"])
df["minute"] = df["minute"].dt.tz_convert(TZ)
df = df.sort_values("minute").reset_index(drop=True)

# Returns in USD/oz between consecutive minute closes. A gap larger than one minute (the daily
# break, weekends, holidays, missing quotes) is not a minute's movement: mark it and keep it apart.
df["gap_min"] = df["minute"].diff().dt.total_seconds().div(60)
df["dclose"] = df["bid_close"].diff()
df["is_gap"] = df["gap_min"].ne(1.0)
df.loc[df.index[0], "is_gap"] = True

df["ny_date"] = df["minute"].dt.date
df["dow"] = df["minute"].dt.dayofweek  # 0 = Monday
df["bin"] = df["minute"].dt.hour * 2 + (df["minute"].dt.minute >= 30).astype(int)  # 0..47
df["range"] = df["bid_high"] - df["bid_low"]
df["spread"] = df["ask_close"] - df["bid_close"]

wk = df[df["dow"] < 5].copy()  # weekdays on the New York clock
n_days = wk["ny_date"].nunique()
first, last = df["minute"].iloc[0], df["minute"].iloc[-1]

# Realised variance per (day, bin) from non-gap minute returns, then averaged over days that have
# the bin at all, so a bin's typical move is not diluted by days when the market was shut.
mv = wk[~wk["is_gap"]].assign(sq=lambda d: d["dclose"] ** 2)
per_day_bin = mv.groupby(["ny_date", "bin"])["sq"].sum()
bin_var = per_day_bin.groupby("bin").mean()
bin_var_med = per_day_bin.groupby("bin").median()      # robust to one wild day
top_day_share = per_day_bin.groupby("bin").max() / per_day_bin.groupby("bin").sum()
bin_days = per_day_bin.groupby("bin").size()
day_var = per_day_bin.groupby("ny_date").sum()
share = bin_var / bin_var.sum()

g = wk.groupby("bin")
bins = pd.DataFrame({
    "bin": range(48),
}).set_index("bin")
bins["ny_time"] = [f"{b // 2:02d}:{'30' if b % 2 else '00'}" for b in bins.index]
bins["days"] = bin_days.reindex(bins.index).fillna(0).astype(int)
bins["coverage"] = (g["minute"].count() / (30 * n_days)).reindex(bins.index).fillna(0).round(3)
bins["typical_move_usd"] = np.sqrt(bin_var.reindex(bins.index)).round(2)
bins["variance_share"] = share.reindex(bins.index).round(4)
bins["median_day_move_usd"] = np.sqrt(bin_var_med.reindex(bins.index)).round(2)
bins["largest_day_share"] = top_day_share.reindex(bins.index).round(2)
bins["median_range_usd"] = g["range"].median().reindex(bins.index).round(3)
bins["p90_range_usd"] = g["range"].quantile(0.9).reindex(bins.index).round(3)
bins["median_spread_usd"] = g["spread"].median().reindex(bins.index).round(3)
bins["p90_spread_usd"] = g["spread"].quantile(0.9).reindex(bins.index).round(3)
bins["spread_over_range"] = (bins["median_spread_usd"] / bins["median_range_usd"]).round(2)
bins["median_ticks_per_min"] = g["ticks"].median().reindex(bins.index)
bins.to_csv(os.path.join(HERE, "bins.csv"))

overlap = list(range(16, 24))            # 08:00–11:59 New York
asian = list(range(38, 48)) + list(range(0, 6))  # 19:00–02:59 New York

# Per-minute profile, for the shape inside the half hours.
mm = wk[~wk["is_gap"]].assign(hm=wk["minute"].dt.strftime("%H:%M"), sq=wk["dclose"] ** 2)
minutes = mm.groupby("hm").agg(rms_move_usd=("sq", lambda x: np.sqrt(x.mean())),
                               median_abs_move_usd=("dclose", lambda x: x.abs().median()),
                               median_range_usd=("range", "median"), days=("sq", "size")).round(3)
minutes["median_spread_usd"] = wk.assign(hm=wk["minute"].dt.strftime("%H:%M")).groupby("hm")["spread"].median().round(3)
minutes.to_csv(os.path.join(HERE, "minutes.csv"))
top_minutes = minutes[minutes["days"] > n_days * 0.5].sort_values("median_range_usd", ascending=False).head(8)

# The spread around the break, minute by minute: the pre-registered rule measured the 17:00 hour,
# which on this broker is the break itself for all but three weeks. Measure the edges instead.
hm = wk["minute"].dt.strftime("%H:%M")
day_median_spread = wk["spread"].median()
before = wk[hm.between("16:45", "16:59")].groupby(hm)["spread"].median()
after = wk[hm.between("18:00", "18:14")].groupby(hm)["spread"].median()
wide_minutes = pd.concat([before, after])
wide_minutes = wide_minutes[wide_minutes >= 2 * day_median_spread]
asian_p90 = wk[wk["bin"].isin(asian)]["spread"].quantile(0.9)
overlap_p90 = wk[wk["bin"].isin(overlap)]["spread"].quantile(0.9)

# Day of week: typical whole-day move and median spread.
dd = wk[~wk["is_gap"]].assign(sq=lambda d: d["dclose"] ** 2).groupby(["dow", "ny_date"])["sq"].sum()
weekdays = pd.DataFrame({
    "typical_day_move_usd": np.sqrt(dd.groupby("dow").mean()).round(2),
    "days": dd.groupby("dow").size(),
    "median_spread_usd": wk.groupby("dow")["spread"].median().round(3),
})
weekdays.index = ["Mon", "Tue", "Wed", "Thu", "Fri"]
weekdays.to_csv(os.path.join(HERE, "weekdays.csv"))

# The daily break: where the gaps are, and how big the jump across the main one is.
gaps = df[df["is_gap"] & df["gap_min"].between(30, 180)]
gap_bins = gaps["bin"].value_counts().head(3)
break_jump = gaps["dclose"].abs()

# Verdicts as defined in README.md.
day_median_var = bin_var.median()
overlap_share = share.reindex(overlap).sum()
busiest = int(share.idxmax())
c1_overlap = "confirmed" if (overlap_share > 1 / 3 and busiest in overlap) else (
    "refuted" if (overlap_share < 0.25 or busiest not in overlap) else "partly")
asian_var = bin_var.reindex(asian).dropna().sort_values()
c1_asian = "confirmed" if (asian_var.head(4) < day_median_var).all() else "refuted"

spread_med = bins["median_spread_usd"]
roll_bins = [34, 35]                     # 17:00–17:59 New York
roll_ratio = spread_med.reindex(roll_bins).max() / spread_med.median()
asian_ratio = spread_med.reindex(asian).median() / spread_med.reindex(overlap).median()
c2_hits = int(roll_ratio >= 2) + int(asian_ratio >= 1.5)
c2 = {2: "confirmed", 1: "partly", 0: "refuted"}[c2_hits]

full = bins[bins["days"] > n_days * 0.5]
top = full.sort_values("typical_move_usd", ascending=False).head(6)
top_robust = full.sort_values("median_day_move_usd", ascending=False).head(6)
low = bins[bins["days"] > n_days * 0.5].sort_values("typical_move_usd").head(6)
best_ratio = bins[bins["days"] > n_days * 0.5].sort_values("spread_over_range").head(6)
worst_ratio = bins[bins["days"] > n_days * 0.5].sort_values("spread_over_range", ascending=False).head(6)

def table(frame, cols):
    head = "| New York | " + " | ".join(cols) + " |\n|---|" + "---|" * len(cols) + "\n"
    rows = "".join(f"| {r.ny_time} | " + " | ".join(str(getattr(r, c)) for c in cols) + " |\n" for r in frame.itertuples())
    return head + rows

lines = [
    "# Results: gold by the half hour, New York clock",
    "",
    f"Data: {BARS.replace(os.path.expanduser('~'), '~')}, {len(df):,} minutes from {first:%Y-%m-%d %H:%M} to {last:%Y-%m-%d %H:%M} New York; {n_days} weekdays. "
    f"Generated by profile.py on {pd.Timestamp.now(tz='UTC'):%Y-%m-%d %H:%M} UTC. Definitions in README.md.",
    "",
    "## Verdicts",
    "",
    f"- Claim 1, overlap (08:00–12:00 NY carries > 33% of daily variance and the busiest half hour): share = **{overlap_share:.1%}** "
    f"(pro rata would be 16.7%), busiest half hour = **{bins.loc[busiest, 'ny_time']}** → **{c1_overlap}**.",
    f"- Claim 1, Asian quiet (four quietest bins in 19:00–03:00 NY all below the day's median bin): quietest four = "
    f"{', '.join(bins.loc[b, 'ny_time'] for b in asian_var.head(4).index)} at "
    f"{', '.join(f'{np.sqrt(v):.2f}' for v in asian_var.head(4))} USD vs median bin {np.sqrt(day_median_var):.2f} USD → **{c1_asian}**.",
    f"- Claim 2, spread, rule as written (17:00–18:00 NY median spread ≥ 2× day median: {roll_ratio:.2f}×, on {int(bins.loc[34, 'days'])} days only, because that hour is the broker's break; "
    f"Asian median spread ≥ 1.5× overlap's: {asian_ratio:.2f}×) → **{c2}**.",
    f"- Claim 2, measured at the edges of the break instead (added after seeing that the 17:00 hour is the break): minutes in 16:45–16:59 and 18:00–18:14 whose median spread is ≥ 2× the day's median of {day_median_spread:.2f}: "
    + (", ".join(f"{k} ({v:.2f})" for k, v in wide_minutes.items()) if len(wide_minutes) else "none")
    + f". 90th-percentile spread, Asian window {asian_p90:.2f} vs overlap {overlap_p90:.2f}.",
    "",
    "## Movement",
    "",
    f"Typical whole-day move (root of mean daily realised variance, weekdays): **{np.sqrt(day_var.mean()):.1f} USD/oz**; median day {np.sqrt(day_var.median()):.1f}.",
    "",
    "Busiest half hours (typical move = root mean realised variance in the bin, USD/oz):",
    "",
    table(top, ["typical_move_usd", "variance_share", "largest_day_share", "median_range_usd", "median_ticks_per_min"]),
    "The same ranking on the median day instead of the mean (root of the median daily variance in the bin), which one wild day cannot move:",
    "",
    table(top_robust, ["median_day_move_usd", "typical_move_usd", "median_range_usd"]),
    "Quietest half hours (bins present on more than half the days):",
    "",
    table(low, ["typical_move_usd", "variance_share", "median_range_usd", "median_spread_usd", "median_ticks_per_min"]),
    "Minutes with the largest median 1-minute range (the shape inside the half hours; minutes.csv has all 1,440):",
    "",
    "| minute NY | median range | median abs move | RMS move | days |\n|---|---|---|---|---|\n"
    + "".join(f"| {r.Index} | {r.median_range_usd} | {r.median_abs_move_usd} | {r.rms_move_usd} | {r.days} |\n" for r in top_minutes.itertuples()),
    "## Cost against movement",
    "",
    "Cheapest half hours to trade (median spread as a multiple of the median 1-minute range):",
    "",
    table(best_ratio, ["spread_over_range", "median_spread_usd", "median_range_usd"]),
    "Most expensive:",
    "",
    table(worst_ratio, ["spread_over_range", "median_spread_usd", "median_range_usd"]),
    "## The daily break",
    "",
    f"Gaps of 30–180 minutes occur {len(gaps)} times; the bins they end in: "
    + ", ".join(f"{bins.loc[b, 'ny_time']} ({n})" for b, n in gap_bins.items())
    + f". Absolute jump across them: median {break_jump.median():.2f} USD, 90th percentile {break_jump.quantile(0.9):.2f}, max {break_jump.max():.2f}.",
    "",
    "## Day of week",
    "",
    weekdays.to_markdown(),
    "",
    "Full table: bins.csv (48 rows). Not advice; one broker, one instrument, eight months.",
]
open(os.path.join(HERE, "results.md"), "w").write("\n".join(lines) + "\n")
print("\n".join(lines[:12]))