"""Backtest the live price forecast against what the price actually did.

The forecast being tested is the one used live:
    centre = the current price (a random walk)
    range  = current price x exp(+/- z * sigma * sqrt(h)), where sigma is the
             standard deviation of hourly log returns over the previous 168
             hours

It is issued at every hour of the test period and scored against the real
price h hours later. Nothing about the test period is used to set the range:
sigma at hour t comes only from hours before t.

A probabilistic forecast is right when its stated confidence matches the
outcome: an 80% range should contain the real price 80% of the time, a 50%
range half the time. Too often means the range is wider than it needs to be;
too rarely means it understates the risk.

Usage:
    python scripts/backtest_price_range.py
"""
from __future__ import annotations

import json
import math
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

ROOT = Path(__file__).resolve().parents[1]
HORIZONS = (1, 6, 24, 72)
LEVELS = (0.5, 0.8, 0.95)
LOOKBACK = 168
START = "2023-01-01"
SURFACE, INK, INK_2, INK_3, GRID = "#fcfcfb", "#0b0b0b", "#52514e", "#8b8a85", "#e6e5e1"
SERIES = ["#2a78d6", "#eb6834", "#1baf7a", "#4a3aa7"]       # validated order, slots 1,2,3,7


def z_for(level: float) -> float:
    """Two-sided normal quantile without scipy (blocked on some machines)."""
    p = 0.5 + level / 2
    lo, hi = 0.0, 10.0
    for _ in range(80):
        mid = (lo + hi) / 2
        if 0.5 * math.erfc(-mid / math.sqrt(2)) < p:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


def load() -> pd.Series:
    btc = pd.read_csv(ROOT / "outputs" / "volatility" / "btc_1h.csv",
                      usecols=["timestamp", "close"])
    btc["timestamp"] = pd.to_datetime(btc["timestamp"], utc=True)
    close = btc.set_index("timestamp")["close"].sort_index()
    return close[~close.index.duplicated()].asfreq("1h")


def backtest(close: pd.Series) -> pd.DataFrame:
    logp = np.log(close)
    ret = logp.diff()
    sigma = ret.rolling(LOOKBACK, min_periods=int(LOOKBACK * 0.9)).std()   # known at t
    rows = []
    for h in HORIZONS:
        future = logp.shift(-h)
        frame = pd.DataFrame({"logp": logp, "future": future, "sigma": sigma}).dropna()
        frame = frame[frame.index >= START]
        z_score = (frame["future"] - frame["logp"]) / (frame["sigma"] * math.sqrt(h))
        price, actual = np.exp(frame["logp"]), np.exp(frame["future"])
        for t, zs, p, a, s in zip(frame.index, z_score, price, actual, frame["sigma"]):
            rows.append((h, t, zs, p, a, s))
    return pd.DataFrame(rows, columns=["h", "t", "z", "price", "actual", "sigma"])


def summarise(bt: pd.DataFrame) -> dict:
    out = {}
    for h, g in bt.groupby("h"):
        absz = g["z"].abs()
        err = (g["actual"] - g["price"]).abs()
        z80 = z_for(0.8)
        width = (np.exp(z80 * g["sigma"] * math.sqrt(h)) - np.exp(-z80 * g["sigma"] * math.sqrt(h))) * g["price"]
        by_year = {int(y): float((gy["z"].abs() <= z80).mean())
                   for y, gy in g.groupby(g["t"].dt.year)}
        weekend = g["t"].dt.dayofweek >= 5
        out[f"{h}h"] = {
            "forecasts": int(len(g)),
            "coverage": {f"{int(l * 100)}%": float((absz <= z_for(l)).mean()) for l in LEVELS},
            "coverage_80_by_year": by_year,
            "coverage_80_weekday": float((absz[~weekend] <= z80).mean()),
            "coverage_80_weekend": float((absz[weekend] <= z80).mean()),
            "mean_abs_error_usd": float(err.mean()),
            "median_abs_error_pct": float((err / g["price"]).median() * 100),
            "mean_range_width_usd": float(width.mean()),
            "share_up": float((g["actual"] > g["price"]).mean()),
        }
    return out


def draw(close: pd.Series, bt: pd.DataFrame, summary: dict) -> None:
    plt.rcParams.update({"font.family": "DejaVu Sans", "text.color": INK,
                         "axes.labelcolor": INK_2})
    fig = plt.figure(figsize=(13, 9.5), facecolor=SURFACE)
    grid = fig.add_gridspec(2, 2, height_ratios=[1.2, 1], hspace=0.42, wspace=0.25,
                            left=0.07, right=0.98, top=0.86, bottom=0.08)

    def style(ax):
        ax.set_facecolor(SURFACE)
        for side in ("top", "right"):
            ax.spines[side].set_visible(False)
        for side in ("left", "bottom"):
            ax.spines[side].set_color(GRID)
        ax.tick_params(colors=INK_2, labelsize=8.5, length=0)
        ax.grid(color=GRID, linewidth=0.8)
        ax.set_axisbelow(True)

    # --- the last 60 days: forecast range issued 24h earlier vs the real price ----
    ax = fig.add_subplot(grid[0, :]); style(ax)
    g = bt[bt["h"] == 24].set_index("t")
    end = g.index.max()
    g = g[g.index >= end - pd.Timedelta(days=60)]
    z80 = z_for(0.8)
    shift = pd.Timedelta(hours=24)
    when = g.index + shift                         # the hour each forecast is about
    lo = g["price"] * np.exp(-z80 * g["sigma"] * math.sqrt(24))
    hi = g["price"] * np.exp(z80 * g["sigma"] * math.sqrt(24))
    ax.fill_between(when, lo, hi, color=SERIES[0], alpha=0.16, linewidth=0,
                    label="80% range, forecast 24h earlier")
    ax.plot(when, g["actual"], color=INK, linewidth=1.3, label="actual price")
    miss = (g["actual"] < lo) | (g["actual"] > hi)
    ax.scatter(when[miss], g["actual"][miss], color=SERIES[1], s=10, zorder=3,
               label=f"outside the range ({miss.mean():.0%} of hours)")
    ax.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(lambda v, _: f"${v / 1000:,.0f}k"))
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%d %b"))
    ax.legend(loc="upper left", frameon=False, fontsize=8.5)
    ax.set_title("The last 60 days of the test: 24-hour forecast range vs the real price",
                 loc="left", fontsize=11, fontweight="bold")

    # --- calibration: stated confidence vs how often it held ------------------------
    ax = fig.add_subplot(grid[1, 0]); style(ax)
    ax.plot([40, 100], [40, 100], color=INK_3, linewidth=1, linestyle=(0, (4, 2)))
    for i, h in enumerate(HORIZONS):
        cov = summary[f"{h}h"]["coverage"]
        xs = [int(k[:-1]) for k in cov]
        ys = [v * 100 for v in cov.values()]
        ax.plot(xs, ys, color=SERIES[i], linewidth=2, marker="o", markersize=5,
                markeredgecolor=SURFACE, label=f"{h}h ahead")
    ax.set_xlim(45, 100); ax.set_ylim(45, 100)
    ax.set_xticks([50, 80, 95]); ax.set_yticks([50, 60, 70, 80, 90, 95, 100])
    ax.set_xlabel("stated confidence of the range, %", fontsize=8.5)
    ax.set_ylabel("how often the real price landed inside, %", fontsize=8.5)
    ax.text(97, 49, "on the dashed line = honest", ha="right", fontsize=8, color=INK_2)
    ax.legend(loc="upper left", frameon=False, fontsize=8.5)
    ax.set_title("Are the ranges honest?", loc="left", fontsize=11, fontweight="bold")

    # --- 80% coverage by year -----------------------------------------------------------
    ax = fig.add_subplot(grid[1, 1]); style(ax)
    years = sorted(summary["24h"]["coverage_80_by_year"])
    width = 0.2
    for i, h in enumerate(HORIZONS):
        vals = [summary[f"{h}h"]["coverage_80_by_year"].get(y, np.nan) * 100 for y in years]
        ax.bar(np.arange(len(years)) + (i - 1.5) * width, vals, width * 0.9, color=SERIES[i],
               label=f"{h}h")
    ax.axhline(80, color=INK_2, linewidth=1, linestyle=(0, (4, 2)))
    ax.set_xticks(range(len(years)))
    ax.set_xticklabels([str(y) + (" (to Sep)" if y == years[-1] else "") for y in years])
    ax.set_ylim(60, 100)
    ax.set_ylabel("80% range: share of real prices inside", fontsize=8.5)
    ax.legend(loc="upper center", frameon=False, fontsize=8, ncol=4)
    ax.set_title("Does it hold every year?", loc="left", fontsize=11, fontweight="bold")

    n = summary["24h"]["forecasts"]
    fig.suptitle("Backtest: the live price forecast against the real price", x=0.07, ha="left",
                 fontsize=15, fontweight="bold", y=0.975)
    fig.text(0.07, 0.905,
             f"Issued every hour, {START[:4]} – Sep 2026 ({n:,} forecasts per horizon). Centre = "
             "current price; range from the previous 7 days' hourly volatility only, so no "
             "forecast sees its own future.",
             fontsize=9.3, color=INK_2, ha="left")
    fig.savefig(ROOT / "outputs" / "figures" / "fig11_price_backtest.png", dpi=150,
                facecolor=SURFACE)
    plt.close(fig)
    print("written: outputs/figures/fig11_price_backtest.png")


def main() -> None:
    close = load()
    bt = backtest(close)
    summary = summarise(bt)
    draw(close, bt, summary)
    (ROOT / "outputs" / "price_backtest.json").write_text(json.dumps({
        "method": f"centre = current price; range = price*exp(+/-z*sigma*sqrt(h)), sigma = std "
                  f"of the previous {LOOKBACK} hourly log returns; scored every hour from {START}",
        "horizons": summary}, indent=2))
    for h, s in summary.items():
        cov = "  ".join(f"{k}->{v:.1%}" for k, v in s["coverage"].items())
        print(f"{h:>4}: {s['forecasts']:,} forecasts | {cov} | weekday {s['coverage_80_weekday']:.1%} "
              f"weekend {s['coverage_80_weekend']:.1%} | median miss {s['median_abs_error_pct']:.2f}% "
              f"| up {s['share_up']:.1%}")
        print(f"       80% by year: " + "  ".join(f"{y}:{v:.0%}" for y, v in s["coverage_80_by_year"].items()))


if __name__ == "__main__":
    main()