"""Revision analysis (editorial-board review, 2026-10-17) for UQO WP 8.

Re-derives every table and figure of the revised paper from the *local* raw data
(``data/raw/``) without any network call.  Compared with ``run_analysis.py`` the
script

* labels each Polymarket daily observation by the UTC day it closes (the CLOB
  ``prices-history`` end-point stamps the daily point at 00:00 UTC of the *next*
  day), so that Polymarket closes, U.S. market closes and announcement dates share
  a common calendar;
* flags the listing day of the Fed group (10 Dec 2025, seeded quotes) and reports
  the completeness statistics with and without it;
* tests the no-arbitrage restriction on the price of the complete set against
  discount-adjusted bounds built from the effective federal funds rate (FRED DFF);
* replaces the per-contract CAR event study (2–3 events) by an announcement-day
  analysis of the risk-neutral expected policy change (14 scheduled releases);
* re-runs Granger tests on a trading-day grid with information-criterion lag
  selection (the original code reported the minimum p-value over 1–5 lags) and
  adds the 10-year Treasury yield;
* re-estimates GARCH/GJR/HAR models with the listing day excluded and QML errors;
* generates descriptive, order-book and election tables from the raw metadata.

Usage:  python3 scripts/02_revision_analysis.py          (≈ 1 min, offline)
Outputs: results/tables/*.tex, results/revision_stats.json, figures/fig*.pdf|png
"""

from __future__ import annotations

import json
import sys
import warnings
from datetime import date
from pathlib import Path

import numpy as np
import pandas as pd
from scipy import stats

ROOT = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(ROOT))

from src.config import (  # noqa: E402
    CPI_RELEASE_DATES,
    DATA_RAW_FRED,
    DATA_RAW_POLY,
    DATA_RAW_YAHOO,
    ELECTION_2024_EVENT_ID,
    FED_JUNE_2026_EVENT_ID,
    FIGURES,
    FOMC_DATES,
    NFP_RELEASE_DATES,
    TABLES,
    ELECTION_KEY_DATES,
)

warnings.filterwarnings("ignore")

import matplotlib  # noqa: E402

matplotlib.use("Agg")
import matplotlib.dates as mdates  # noqa: E402
import matplotlib.pyplot as plt  # noqa: E402

# Okabe–Ito palette (CVD-safe, validated), fixed order per state
PAL = ["#0072B2", "#D55E00", "#009E73", "#CC79A7", "#56B4E9", "#E69F00", "#000000"]
LS = ["-", "--", "-", "-.", ":", "-", "--"]
plt.rcParams.update({
    "font.family": "serif", "font.serif": ["Times New Roman", "DejaVu Serif"],
    "font.size": 9, "axes.titlesize": 10, "axes.labelsize": 9, "legend.fontsize": 8,
    "xtick.labelsize": 8, "ytick.labelsize": 8, "axes.spines.top": False,
    "axes.spines.right": False, "axes.grid": True, "grid.alpha": 0.25,
    "grid.linestyle": "--", "lines.linewidth": 1.3, "figure.dpi": 150,
})

STATES = ["50+ bps decrease", "25 bps decrease", "No change", "25 bps increase", "50+ bps increase"]
STATE_BPS = np.array([-50.0, -25.0, 0.0, 25.0, 50.0])
RESOLUTION_FED = pd.Timestamp("2026-06-17")
RESOLUTION_ELECTION = pd.Timestamp("2024-11-05")
TICK = 0.001
STATS: dict = {}


# ----------------------------------------------------------------------------
# helpers
# ----------------------------------------------------------------------------
def save_fig(fig, name: str) -> None:
    FIGURES.mkdir(parents=True, exist_ok=True)
    for ext in ("pdf", "png"):
        fig.savefig(FIGURES / f"{name}.{ext}", bbox_inches="tight", dpi=300)
    plt.close(fig)


def save_table(tex: str, name: str) -> None:
    TABLES.mkdir(parents=True, exist_ok=True)
    (TABLES / f"{name}.tex").write_text(tex)


def esc(s: str) -> str:
    return (str(s).replace("&", r"\&").replace("%", r"\%").replace("$", r"\$")
            .replace("_", r"\_").replace("#", r"\#"))


def f(x, d=3, blank="--"):
    if x is None or (isinstance(x, float) and np.isnan(x)):
        return blank
    return f"{x:,.{d}f}"


def stars(p):
    if p is None or np.isnan(p):
        return ""
    return "$^{***}$" if p < 0.01 else "$^{**}$" if p < 0.05 else "$^{*}$" if p < 0.10 else ""


def nw_mean_test(x: pd.Series, lags: int = 5, mu0: float = 0.0):
    """HAC (Newey–West) t-test of H0: E[x] = mu0."""
    import statsmodels.api as sm
    x = x.dropna()
    mod = sm.OLS(x.values - mu0, np.ones((len(x), 1))).fit(cov_type="HAC", cov_kwds={"maxlags": lags})
    return float(mod.params[0] + mu0), float(mod.bse[0]), float(mod.tvalues[0]), float(mod.pvalues[0])


def load_history(token: str) -> pd.Series:
    """Daily close labelled by the UTC day it closes (point stamped 00:00 UTC of day t+1 → day t)."""
    path = DATA_RAW_POLY / f"price_history_{token[:16]}.json"
    raw = json.load(open(path)) if path.exists() else []
    if not raw:
        return pd.Series(dtype=float)
    df = pd.DataFrame(raw)
    ts = pd.to_datetime(df["t"], unit="s", utc=True)
    df["day"] = (ts - pd.Timedelta(minutes=5)).dt.normalize().dt.tz_localize(None)
    df["p"] = df["p"].astype(float)
    return df.drop_duplicates("day", keep="last").set_index("day")["p"].sort_index()


def load_event(event_id: int) -> tuple[dict, list[dict]]:
    ev = json.load(open(DATA_RAW_POLY / f"event_{event_id}.json"))
    out = []
    for m in ev["markets"]:
        toks = json.loads(m.get("clobTokenIds", "[]"))
        out.append({"id": m["id"], "name": m.get("groupItemTitle", m["question"][:40]),
                    "question": m["question"], "yes_token": toks[0] if toks else "",
                    "volume": float(m.get("volumeNum") or 0), "liquidity": float(m.get("liquidityNum") or 0),
                    "closed": m.get("closed", False), "start": m.get("startDate", ""),
                    "end": m.get("endDate", ""), "spread": m.get("spread"),
                    "best_bid": m.get("bestBid"), "best_ask": m.get("bestAsk")})
    return ev, out


def build_panel(markets: list[dict], order: list[str] | None = None) -> pd.DataFrame:
    cols = {m["name"]: load_history(m["yes_token"]) for m in markets if m["yes_token"]}
    cols = {k: v for k, v in cols.items() if len(v)}
    panel = pd.DataFrame(cols).sort_index()
    if order:
        panel = panel[[c for c in order if c in panel.columns]]
    return panel.ffill(limit=3)


# ----------------------------------------------------------------------------
# 1. Data
# ----------------------------------------------------------------------------
fed_ev, fed_mk = load_event(FED_JUNE_2026_EVENT_ID)
fed = build_panel(fed_mk, STATES)
el_ev, el_mk = load_event(ELECTION_2024_EVENT_ID)
el_order = sorted(el_mk, key=lambda m: -m["volume"])
election = build_panel(el_order, [m["name"] for m in el_order])

dff = pd.read_csv(DATA_RAW_FRED / "fed_funds_rate.csv", parse_dates=["date"], index_col="date")["DFF"]
dff = dff.asfreq("D").ffill()
dff = dff.reindex(pd.date_range(dff.index.min(), "2026-06-30")).ffill()


def yahoo(ticker: str) -> pd.Series:
    df = pd.read_csv(DATA_RAW_YAHOO / f"{ticker}.csv", parse_dates=["date"], index_col="date")
    df.index = pd.to_datetime(df.index, utc=True).normalize().tz_localize(None)
    return df["Close"].astype(float)


spy, vix, tnx, btc = yahoo("SPY"), yahoo("VIX"), yahoo("TNX"), yahoo("BTC_USD")

listing_day = fed.index.min()
fed_complete = fed.dropna()
fed_main = fed_complete[fed_complete.index > listing_day]  # listing day excluded
sigma_all = fed_complete.sum(axis=1)
sigma = fed_main.sum(axis=1)
STATS["fed"] = {
    "first_day": str(fed.index.min().date()), "last_day": str(fed.index.max().date()),
    "n_rows_union": int(len(fed)), "n_complete": int(len(fed_complete)), "n_main": int(len(fed_main)),
    "listing_sigma": float(sigma_all.iloc[0]), "listing_prices": fed_complete.iloc[0].round(4).to_dict(),
    "n_weekend_days": int((fed_main.index.dayofweek >= 5).sum()),
    "calendar_span_days": int((fed.index.max() - fed.index.min()).days) + 1,
    "p_nochange_first": float(fed_main["No change"].iloc[0]), "p_nochange_last": float(fed_main["No change"].iloc[-1]),
    "p_cut25_first": float(fed_main["25 bps decrease"].iloc[0]), "p_cut25_last": float(fed_main["25 bps decrease"].iloc[-1]),
    "dff_start": float(dff.loc[fed_main.index.min()]), "dff_end": float(dff.loc[fed_main.index.max()]),
    "dff_before_listing": float(dff.loc[listing_day - pd.Timedelta(days=1)]),
    "min_price_by_state": fed_main.min().round(4).to_dict(),
    "volume_usd": {m["name"]: m["volume"] for m in fed_mk},
    "liquidity_usd": {m["name"]: m["liquidity"] for m in fed_mk},
    "total_volume_usd": float(sum(m["volume"] for m in fed_mk)),
}

# ----------------------------------------------------------------------------
# 2. Table 1 — markets analysed (descriptive)
# ----------------------------------------------------------------------------
rows = []
for m in fed_mk:
    s = fed[m["name"]].dropna().iloc[1:]  # drop seeded listing quote
    rows.append(f"    {esc(m['name'])} & {f(m['volume'],0)} & {f(m['liquidity'],0)} & {len(s)} & "
                f"{f(s.iloc[0])} & {f(s.mean())} & {f(s.min())} & {f(s.max())} & {f(s.iloc[-1])} \\\\")
fed_rows = "\n".join(rows)
rows = []
for m in el_order:
    s = election[m["name"]].dropna().iloc[1:] if m["name"] in election else pd.Series(dtype=float)
    rows.append(f"    {esc(m['name'])} & {f(m['volume']/1e6,1)} & -- & {len(s)} & "
                f"{f(s.iloc[0]) if len(s) else '--'} & {f(s.mean()) if len(s) else '--'} & {f(s.min()) if len(s) else '--'} & "
                f"{f(s.max()) if len(s) else '--'} & {f(s.iloc[-1]) if len(s) else '--'} \\\\")
el_rows = "\n".join(rows)
tex = rf"""\begin{{table}}[htbp]
\centering\footnotesize
\begin{{threeparttable}}
\caption{{Markets analysed: Polymarket metadata and daily closing prices}}
\label{{tab:descriptive_stats}}
\setlength{{\tabcolsep}}{{4pt}}
\begin{{tabular}}{{l r r r r r r r r}}
\toprule
Contract & Volume & Liquidity & $T$ & First & Mean & Min & Max & Last \\
\midrule
\multicolumn{{9}}{{l}}{{\textit{{Panel A: Fed decision, June 2026 FOMC meeting (event 101772; volume and liquidity in USD)}}}} \\
{fed_rows}
\midrule
\multicolumn{{9}}{{l}}{{\textit{{Panel B: Presidential Election Winner 2024 (event 903193; volume, USD million)}}}} \\
{el_rows}
\bottomrule
\end{{tabular}}
\begin{{tablenotes}}
\small
\item \textit{{Notes:}} Volume is cumulative traded notional reported by the Gamma API at the snapshot date (22 May 2026); liquidity is the resting order-book depth reported by the same end-point (not available for closed markets). $T$ is the number of daily closes after exclusion of each contract's first recorded point, which carries the quote seeded at listing (0.50 for most election contracts; the Fed seeds are reported in Section~\ref{{sec:data}}); ``First'' is the first regular close. The maximum of 0.500 for the AOC contract is an isolated print on 9 April 2024. Prices are in USDC per \$1 payoff. Panel B contracts are ordered by volume; the sample ends at the close of 5 November 2024 (00:00 UTC on 6 November), before results were known.
\end{{tablenotes}}
\end{{threeparttable}}
\end{{table}}"""
save_table(tex, "table1_descriptive_stats")

# ----------------------------------------------------------------------------
# 3. Table 2 — summary statistics of state prices and daily changes
# ----------------------------------------------------------------------------
rows = []
for c in STATES:
    s = fed_main[c]
    d = s.diff().dropna()
    jb = stats.jarque_bera(d)
    adf = __import__("statsmodels.tsa.stattools", fromlist=["adfuller"]).adfuller(d, autolag="AIC")
    rows.append(f"    {esc(c)} & {f(s.mean())} & {f(s.std())} & {f(s.min())} & {f(s.max())} & "
                f"{f(d.mean()*100,3)} & {f(d.std()*100,2)} & {f(d.skew(),2)} & {f(d.kurtosis(),1)} & {f(adf[0],2)} \\\\")
sig = sigma
d = sig.diff().dropna()
rows.append(f"    State-price sum $\\Sigma_t$ & {f(sig.mean(),4)} & {f(sig.std(),4)} & {f(sig.min(),4)} & {f(sig.max(),4)} & "
            f"{f(d.mean()*100,3)} & {f(d.std()*100,2)} & {f(d.skew(),2)} & {f(d.kurtosis(),1)} & "
            f"{f(__import__('statsmodels.tsa.stattools', fromlist=['adfuller']).adfuller(d, autolag='AIC')[0],2)} \\\\")
tex = rf"""\begin{{table}}[htbp]
\centering\small
\begin{{threeparttable}}
\caption{{Summary statistics: Fed June 2026 state prices (11 December 2025 -- 22 May 2026)}}
\label{{tab:summary_stats}}
\setlength{{\tabcolsep}}{{4pt}}
\begin{{tabular}}{{l rrrr rrrr r}}
\toprule
 & \multicolumn{{4}}{{c}}{{Level $p_{{j,t}}$}} & \multicolumn{{5}}{{c}}{{Daily change $\Delta p_{{j,t}}$ (percentage points)}} \\
\cmidrule(lr){{2-5}}\cmidrule(lr){{6-10}}
Contract & Mean & SD & Min & Max & Mean & SD & Skew & Ex.\ kurt & ADF \\
\midrule
{chr(10).join(rows)}
\bottomrule
\end{{tabular}}
\begin{{tablenotes}}
\small
\item \textit{{Notes:}} $T={len(fed_main)}$ daily closes (calendar days, Polymarket trades continuously) after exclusion of the listing day. Levels are in USDC; changes are in percentage points of the \$1 payoff. ADF is the augmented Dickey--Fuller statistic (constant, AIC lag) for the daily changes; the 1\% critical value is $-3.48$. Prices are bounded in $[0,1]$ and converge to 0 or 1 at resolution, so level series are not stationary by construction.
\end{{tablenotes}}
\end{{threeparttable}}
\end{{table}}"""
save_table(tex, "table2_summary_stats")

# ----------------------------------------------------------------------------
# 4. Completeness / no-arbitrage on the complete set (Table 3)
# ----------------------------------------------------------------------------
from statsmodels.tsa.stattools import adfuller  # noqa: E402


def completeness_block(sig: pd.Series, resolution: pd.Timestamp, n_states: int, tick: float, label: str) -> dict:
    dev = sig - 1.0
    tau = (resolution - sig.index).days.values / 365.0
    r = dff.reindex(sig.index).ffill().bfill().values / 100.0
    lower = np.exp(-r * tau)                      # discounted value of the complete set (buy-and-hold bound)
    dev_disc = sig.values - lower
    mean, se, t, p = nw_mean_test(dev)
    adf = adfuller(dev, autolag="AIC", regression="c")
    phi = np.polyfit(dev.values[:-1], dev.values[1:], 1)[0]
    hl = -np.log(2) / np.log(abs(phi)) if 0 < abs(phi) < 1 else np.nan
    band = 0.005  # half a cent on the complete set
    out = {
        "label": label, "T": int(len(sig)), "mean": float(sig.mean()), "sd": float(sig.std()),
        "median": float(sig.median()), "min": float(sig.min()), "max": float(sig.max()),
        "mean_dev": float(dev.mean()), "mean_abs_dev": float(dev.abs().mean()), "max_abs_dev": float(dev.abs().max()),
        "nw_se": se, "nw_t": t, "nw_p": p,
        "share_above_1": float((sig > 1).mean()), "share_below_1": float((sig < 1).mean()),
        "share_within_tickband": float((dev.abs() <= band).mean()),
        "share_within_1pp": float((dev.abs() <= 0.01).mean()), "share_within_2pp": float((dev.abs() <= 0.02).mean()),
        "share_within_5pp": float((dev.abs() <= 0.05).mean()),
        "mean_lower_bound": float(lower.mean()), "min_lower_bound": float(lower.min()),
        "mean_dev_from_lower": float(dev_disc.mean()), "share_inside_band": float(((sig.values >= lower - band) & (sig.values <= 1 + band)).mean()),
        "share_above_upper_band": float((sig.values > 1 + band).mean()),
        "adf_stat": float(adf[0]), "adf_p": float(adf[1]), "ar1": float(phi), "half_life": float(hl),
        "tick_band": band,
        "first_half_mean": float(sig.iloc[: len(sig) // 2].mean()), "second_half_mean": float(sig.iloc[len(sig) // 2:].mean()),
        "last30_mean": float(sig.iloc[-30:].mean()), "last30_sd": float(sig.iloc[-30:].std()),
        "first30_mean": float(sig.iloc[:30].mean()), "first30_sd": float(sig.iloc[:30].std()),
    }
    # overround vs. time-to-resolution regression
    import statsmodels.api as sm
    X = sm.add_constant(tau)
    m = sm.OLS(dev.values, X).fit(cov_type="HAC", cov_kwds={"maxlags": 5})
    out["slope_tau"] = float(m.params[1]); out["slope_tau_t"] = float(m.tvalues[1]); out["const_tau"] = float(m.params[0])
    return out, lower


el_complete_all = election.dropna()
el_complete = el_complete_all.iloc[1:]  # first complete day carries seeded listing quotes
sigma_el = el_complete.sum(axis=1)
sigma_el_all = el_complete_all.sum(axis=1)
C_fed, lower_fed = completeness_block(sigma, RESOLUTION_FED, 5, TICK, "Fed June 2026")
C_fed_all, _ = completeness_block(sigma_all, RESOLUTION_FED, 5, TICK, "Fed incl. listing day")
C_el, lower_el = completeness_block(sigma_el, RESOLUTION_ELECTION, len(el_complete.columns), 0.001, "Election 2024")
C_el_all, _ = completeness_block(sigma_el_all, RESOLUTION_ELECTION, len(el_complete.columns), 0.001, "Election incl. listing day")
STATS["completeness"] = {"fed": C_fed, "fed_incl_listing": C_fed_all, "election": C_el, "election_incl_listing": C_el_all,
                         "election_n_contracts": int(len(el_complete.columns)), "election_T_union": int(len(election)),
                         "election_first": str(election.index.min().date()), "election_last": str(election.index.max().date())}


def row(name, a, b, d=4, pct=False):
    fa = f(a * 100, 1) + r"\%" if pct else f(a, d)
    fb = f(b * 100, 1) + r"\%" if pct else f(b, d)
    return f"    {name} & {fa} & {fb} \\\\"


tex = rf"""\begin{{table}}[htbp]
\centering\small
\begin{{threeparttable}}
\caption{{Price of the complete set of state claims: law-of-one-price tests}}
\label{{tab:completeness}}
\begin{{tabular}}{{l r r}}
\toprule
 & Fed June 2026 & Election 2024 \\
 & ($N=5$ states) & ($N={len(el_complete.columns)}$ candidates) \\
\midrule
\multicolumn{{3}}{{l}}{{\textit{{Panel A: distribution of $\Sigma_t=\sum_j p_{{j,t}}$}}}} \\
{row('Observations $T$', C_fed['T'], C_el['T'], 0)}
{row('Mean', C_fed['mean'], C_el['mean'])}
{row('Median', C_fed['median'], C_el['median'])}
{row('Standard deviation', C_fed['sd'], C_el['sd'])}
{row('Minimum', C_fed['min'], C_el['min'])}
{row('Maximum', C_fed['max'], C_el['max'])}
\midrule
\multicolumn{{3}}{{l}}{{\textit{{Panel B: deviation $\Sigma_t-1$}}}} \\
{row('Mean deviation (overround)', C_fed['mean_dev'], C_el['mean_dev'])}
{row('HAC standard error (5 lags)', C_fed['nw_se'], C_el['nw_se'])}
{row('$t$-statistic, $H_0$: $\\mathbb{{E}}[\\Sigma_t-1]=0$', C_fed['nw_t'], C_el['nw_t'], 2)}
{row('Mean absolute deviation', C_fed['mean_abs_dev'], C_el['mean_abs_dev'])}
{row('Share of days with $\\Sigma_t>1$', C_fed['share_above_1'], C_el['share_above_1'], pct=True)}
{row('Share within $\\pm 0.5$ pp', C_fed['share_within_tickband'], C_el['share_within_tickband'], pct=True)}
{row('Share within $\\pm 1$ pp', C_fed['share_within_1pp'], C_el['share_within_1pp'], pct=True)}
{row('Share within $\\pm 2$ pp', C_fed['share_within_2pp'], C_el['share_within_2pp'], pct=True)}
{row('Share within $\\pm 5$ pp', C_fed['share_within_5pp'], C_el['share_within_5pp'], pct=True)}
\midrule
\multicolumn{{3}}{{l}}{{\textit{{Panel C: discount-adjusted bounds, $e^{{-r_t\tau_t}}\le\Sigma_t\le 1$}}}} \\
{row('Mean lower bound $e^{{-r_t\\tau_t}}$', C_fed['mean_lower_bound'], C_el['mean_lower_bound'])}
{row('Mean $\\Sigma_t-e^{{-r_t\\tau_t}}$', C_fed['mean_dev_from_lower'], C_el['mean_dev_from_lower'])}
{row('Share inside bounds ($\\pm 0.5$ pp)', C_fed['share_inside_band'], C_el['share_inside_band'], pct=True)}
{row('Share above upper bound', C_fed['share_above_upper_band'], C_el['share_above_upper_band'], pct=True)}
\midrule
\multicolumn{{3}}{{l}}{{\textit{{Panel D: dynamics of $\Sigma_t-1$}}}} \\
{row('AR(1) coefficient', C_fed['ar1'], C_el['ar1'], 3)}
{row('Half-life (days)', C_fed['half_life'], C_el['half_life'], 1)}
{row('ADF statistic (constant, AIC lag)', C_fed['adf_stat'], C_el['adf_stat'], 2)}
{row('ADF $p$-value', C_fed['adf_p'], C_el['adf_p'], 3)}
{row('Slope on time-to-resolution $\\tau_t$ (years)', C_fed['slope_tau'], C_el['slope_tau'], 3)}
{row('\\quad HAC $t$-statistic', C_fed['slope_tau_t'], C_el['slope_tau_t'], 2)}
\bottomrule
\end{{tabular}}
\begin{{tablenotes}}
\small
\item \textit{{Notes:}} $\Sigma_t$ is the sum of the daily closing prices of all contracts in the neg-risk group on days when all contracts are quoted. The Fed sample excludes the listing day (10 December 2025, $\Sigma=1.570$ on seeded quotes; including it gives mean {f(C_fed_all['mean'],4)}, SD {f(C_fed_all['sd'],4)}, maximum absolute deviation {f(C_fed_all['max_abs_dev'],3)}); the election sample likewise excludes its first complete day (4 January 2024, $\Sigma={f(C_el_all['max'],3)}$). The lower bound discounts the \$1 payoff of the complete set at the effective federal funds rate (FRED DFF) over the remaining time to resolution; the upper bound is one because a complete set can be minted and its NO legs converted immediately through the neg-risk adapter. HAC standard errors use the Newey--West estimator with 5 lags.
\end{{tablenotes}}
\end{{threeparttable}}
\end{{table}}"""
save_table(tex, "table3_completeness")

# ----------------------------------------------------------------------------
# 5. Order-book snapshot (Table 9)
# ----------------------------------------------------------------------------
ob_rows, ob = [], {}
for m in fed_mk:
    path = DATA_RAW_POLY / f"order_book_{m['yes_token'][:16]}.json"
    if not path.exists():
        continue
    d = json.load(open(path))
    bids = sorted([(float(b["price"]), float(b["size"])) for b in d["bids"]], reverse=True)
    asks = sorted([(float(a["price"]), float(a["size"])) for a in d["asks"]])
    bb, ba = bids[0][0], asks[0][0]
    ob[m["name"]] = {"bid": bb, "ask": ba, "mid": (bb + ba) / 2, "spread": ba - bb,
                     "rel_spread_pct": 100 * (ba - bb) / ((bb + ba) / 2),
                     "depth_best_bid": bids[0][1], "depth_best_ask": asks[0][1],
                     "depth_bid_usd": sum(p * s for p, s in bids), "depth_ask_usd": sum(p * s for p, s in asks),
                     "shares_bid": sum(s for _, s in bids), "shares_ask": sum(s for _, s in asks),
                     "depth_within_1c_bid_usd": sum(p * s for p, s in bids if p >= bb - 0.01),
                     "depth_within_1c_ask_usd": sum(p * s for p, s in asks if p <= ba + 0.01),
                     "timestamp": int(d["timestamp"])}
STATS["orderbook"] = ob
STATS["orderbook"]["sum_bid"] = float(sum(v["bid"] for k, v in ob.items() if isinstance(v, dict)))
STATS["orderbook"]["sum_ask"] = float(sum(v["ask"] for k, v in ob.items() if isinstance(v, dict)))
STATS["orderbook"]["sum_mid"] = float(sum(v["mid"] for k, v in ob.items() if isinstance(v, dict)))
STATS["orderbook"]["snapshot_utc"] = str(pd.to_datetime(list(ob.values())[0]["timestamp"], unit="ms"))
for c in STATES:
    v = ob[c]
    ob_rows.append(f"    {esc(c)} & {v['bid']:.3f} & {v['ask']:.3f} & {v['spread']:.3f} & {f(v['rel_spread_pct'],1)} & "
                   f"{f(v['depth_best_bid'],0)} & {f(v['depth_best_ask'],0)} & {f(v['depth_within_1c_bid_usd']/1e3,0)} & {f(v['depth_within_1c_ask_usd']/1e3,0)} \\\\")
tex = rf"""\begin{{table}}[htbp]
\centering\small
\begin{{threeparttable}}
\caption{{Order-book snapshot, Fed June 2026 contracts (22 May 2026, {STATS['orderbook']['snapshot_utc'][11:16]} UTC)}}
\label{{tab:orderbook}}
\setlength{{\tabcolsep}}{{4pt}}
\begin{{tabular}}{{l rrrr rr rr}}
\toprule
 & Best & Best & Quoted & Relative & \multicolumn{{2}}{{c}}{{Depth at best (shares)}} & \multicolumn{{2}}{{c}}{{Depth within 1\textcent\ (USD 000)}} \\
\cmidrule(lr){{6-7}}\cmidrule(lr){{8-9}}
Contract & bid & ask & spread & spread (\%) & Bid & Ask & Bid & Ask \\
\midrule
{chr(10).join(ob_rows)}
\midrule
Sum of quotes & {STATS['orderbook']['sum_bid']:.3f} & {STATS['orderbook']['sum_ask']:.3f} & & & & & & \\
\bottomrule
\end{{tabular}}
\begin{{tablenotes}}
\small
\item \textit{{Notes:}} Snapshot of the CLOB order book (YES tokens) retrieved with the price histories. Quoted spread is in USDC; the relative spread divides it by the midpoint. Depth within 1\textcent\ is the notional resting within one cent of the best quote. The sum of best bids is the proceeds from selling one complete set at market; the sum of best asks is the cost of buying one.
\end{{tablenotes}}
\end{{threeparttable}}
\end{{table}}"""
save_table(tex, "table9_orderbook")

# ----------------------------------------------------------------------------
# 6. Risk-neutral density, expected policy change, announcement days (Table 4)
# ----------------------------------------------------------------------------
q = fed_main.div(sigma, axis=0)
eq = (q.values @ STATE_BPS)
eq = pd.Series(eq, index=fed_main.index, name="EQ_bps")
var_q = pd.Series(((STATE_BPS[None, :] - eq.values[:, None]) ** 2 * q.values).sum(axis=1), index=fed_main.index)
d_eq = eq.diff().dropna()
d_pnc = fed_main["No change"].diff().dropna() * 100  # pp

def _map(d0: pd.Timestamp) -> pd.Timestamp | None:
    """First available Polymarket close on or after d0 (gaps of at most 2 days)."""
    pos = d_eq.index.searchsorted(d0)
    if pos < len(d_eq.index) and (d_eq.index[pos] - d0).days <= 2:
        return d_eq.index[pos]
    return None


EVENT_MAP = {}
def in_sample(ds):
    out = []
    for x in ds:
        x = pd.Timestamp(x)
        if fed_main.index.min() < x <= fed_main.index.max():
            m = _map(x)
            if m is not None:
                EVENT_MAP[str(m.date())] = str(x.date()); out.append(m)
    return out
fomc = in_sample(FOMC_DATES)
cpi = in_sample(CPI_RELEASE_DATES)
nfp = in_sample(NFP_RELEASE_DATES)
ann = {"FOMC": fomc, "CPI": cpi, "Employment report": nfp}
all_ann = sorted(set(fomc + cpi + nfp))
STATS["events"] = {"event_map": EVENT_MAP, "fomc": [str(x.date()) for x in fomc], "cpi": [str(x.date()) for x in cpi], "nfp": [str(x.date()) for x in nfp],
                   "eq_first": float(eq.iloc[0]), "eq_last": float(eq.iloc[-1]), "eq_min": float(eq.min()), "eq_max": float(eq.max()),
                   "sd_q_first": float(np.sqrt(var_q.iloc[0])), "sd_q_last": float(np.sqrt(var_q.iloc[-1]))}


def cum(series: pd.Series, d0: pd.Timestamp, a: int, b: int) -> float:
    idx = series.index
    pos = idx.searchsorted(d0)
    if pos >= len(idx) or idx[pos] != d0:
        return np.nan
    lo, hi = pos + a, pos + b
    if lo < 0 or hi >= len(idx):
        return np.nan
    return float(series.iloc[lo:hi + 1].sum())


ev_rows, ev_records = [], []
for etype, ds in ann.items():
    for d0 in ds:
        rec = {"type": etype, "date": str(d0.date()),
               "d_eq_0": cum(d_eq, d0, 0, 0), "d_eq_m1p1": cum(d_eq, d0, -1, 1), "d_eq_m3p3": cum(d_eq, d0, -3, 3),
               "d_pnc_0": cum(d_pnc, d0, 0, 0), "d_pnc_m1p1": cum(d_pnc, d0, -1, 1), "d_pnc_m3p3": cum(d_pnc, d0, -3, 3),
               "eq_level": float(eq.loc[d0]) if d0 in eq.index else np.nan}
        ev_records.append(rec)
        ev_rows.append(f"    {etype} & {d0.strftime('%d %b %Y')} & {f(rec['eq_level'],1)} & {f(rec['d_eq_0'],2)} & {f(rec['d_eq_m1p1'],2)} & "
                       f"{f(rec['d_eq_m3p3'],2)} & {f(rec['d_pnc_0'],2)} & {f(rec['d_pnc_m1p1'],2)} & {f(rec['d_pnc_m3p3'],2)} \\\\")
STATS["events"]["records"] = ev_records

# announcement-day vs other-day comparison
is_ann = d_eq.index.isin(all_ann)
is_wkend = d_eq.index.dayofweek >= 5


def compare(x: pd.Series, mask: np.ndarray, label: str) -> dict:
    a, b = x[mask].abs(), x[~mask].abs()
    u = stats.mannwhitneyu(a, b, alternative="greater")
    lev = stats.levene(x[mask], x[~mask], center="median")
    return {"label": label, "n_a": int(len(a)), "n_b": int(len(b)), "mean_abs_a": float(a.mean()), "mean_abs_b": float(b.mean()),
            "median_abs_a": float(a.median()), "median_abs_b": float(b.median()),
            "ratio_mean": float(a.mean() / b.mean()), "sd_a": float(x[mask].std()), "sd_b": float(x[~mask].std()),
            "var_ratio": float(x[mask].var() / x[~mask].var()), "mw_p": float(u.pvalue), "bf_stat": float(lev.statistic), "bf_p": float(lev.pvalue),
            "mean_a": float(x[mask].mean()), "mean_b": float(x[~mask].mean())}


cmp_rows = []
COMP = {}
for lab, mask in [("All scheduled releases", is_ann), ("FOMC", d_eq.index.isin(fomc)), ("CPI", d_eq.index.isin(cpi)),
                  ("Employment report", d_eq.index.isin(nfp))]:
    c1 = compare(d_eq, mask, lab); c2 = compare(d_pnc, mask, lab)
    COMP[lab] = {"eq": c1, "pnc": c2}
    cmp_rows.append(f"    {lab} & {c1['n_a']} & {f(c1['mean_abs_a'],2)} & {f(c1['mean_abs_b'],2)} & {f(c1['var_ratio'],2)} & {f(c1['bf_p'],3)} & {f(c1['mw_p'],3)} & "
                    f"{f(c2['mean_abs_a'],2)} & {f(c2['mean_abs_b'],2)} & {f(c2['var_ratio'],2)} & {f(c2['bf_p'],3)} & {f(c2['mw_p'],3)} \\\\")
# weekend vs weekday (non-announcement days only)
wk_mask = is_wkend & ~is_ann
c1 = compare(d_eq[~is_ann], wk_mask[~is_ann], "Weekend"); c2 = compare(d_pnc[~is_ann], wk_mask[~is_ann], "Weekend")
COMP["Weekend"] = {"eq": c1, "pnc": c2}
# TNX benchmark: |Δ10y yield| on same announcement days
tnx_d = (tnx.diff() * 100).dropna()
tnx_d = tnx_d[(tnx_d.index > fed_main.index.min()) & (tnx_d.index <= fed_main.index.max())]
ct = compare(tnx_d, tnx_d.index.isin(all_ann), "TNX all releases")
COMP["TNX"] = ct
STATS["announcement"] = COMP

tex = rf"""\begin{{table}}[htbp]
\centering\footnotesize
\begin{{threeparttable}}
\caption{{Scheduled macroeconomic releases and the risk-neutral expected policy change}}
\label{{tab:event_study}}
\setlength{{\tabcolsep}}{{4pt}}
\begin{{tabular}}{{l l r rrr rrr}}
\toprule
\multicolumn{{9}}{{l}}{{\textit{{Panel A: individual releases in the sample}}}} \\
\addlinespace[2pt]
 & & $\mathbb{{E}}^{{\mathbb{{Q}}}}_t[\Delta r]$ & \multicolumn{{3}}{{c}}{{$\Delta \mathbb{{E}}^{{\mathbb{{Q}}}}[\Delta r]$ (bps)}} & \multicolumn{{3}}{{c}}{{$\Delta p_{{\text{{No change}}}}$ (pp)}} \\
\cmidrule(lr){{4-6}}\cmidrule(lr){{7-9}}
Release & Date & level (bps) & $[0]$ & $[-1,+1]$ & $[-3,+3]$ & $[0]$ & $[-1,+1]$ & $[-3,+3]$ \\
\midrule
{chr(10).join(ev_rows)}
\bottomrule
\end{{tabular}}\par\bigskip
\begin{{tabular}}{{l r rrrrr rrrrr}}
\toprule
\multicolumn{{12}}{{l}}{{\textit{{Panel B: absolute daily changes on release days versus all other days}}}} \\
\addlinespace[2pt]
 & & \multicolumn{{5}}{{c}}{{$|\Delta \mathbb{{E}}^{{\mathbb{{Q}}}}[\Delta r]|$ (bps)}} & \multicolumn{{5}}{{c}}{{$|\Delta p_{{\text{{No change}}}}|$ (pp)}} \\
\cmidrule(lr){{3-7}}\cmidrule(lr){{8-12}}
Days & $K$ & Ann. & Other & VR & $p_{{BF}}$ & $p_{{MW}}$ & Ann. & Other & VR & $p_{{BF}}$ & $p_{{MW}}$ \\
\midrule
"""
for lab in ["All scheduled releases", "FOMC", "CPI", "Employment report"]:
    c1, c2 = COMP[lab]["eq"], COMP[lab]["pnc"]
    tex += (f"    {lab} & {c1['n_a']} & {f(c1['mean_abs_a'],2)} & {f(c1['mean_abs_b'],2)} & {f(c1['var_ratio'],2)} & {f(c1['bf_p'],3)} & {f(c1['mw_p'],3)} & "
            f"{f(c2['mean_abs_a'],2)} & {f(c2['mean_abs_b'],2)} & {f(c2['var_ratio'],2)} & {f(c2['bf_p'],3)} & {f(c2['mw_p'],3)} \\\\\n")
c1, c2 = COMP["Weekend"]["eq"], COMP["Weekend"]["pnc"]
tex += (f"    Weekend vs.\\ weekday (no release) & {c1['n_a']} & {f(c1['mean_abs_a'],2)} & {f(c1['mean_abs_b'],2)} & {f(c1['var_ratio'],2)} & {f(c1['bf_p'],3)} & {f(c1['mw_p'],3)} & "
        f"{f(c2['mean_abs_a'],2)} & {f(c2['mean_abs_b'],2)} & {f(c2['var_ratio'],2)} & {f(c2['bf_p'],3)} & {f(c2['mw_p'],3)} \\\\\n")
tex += rf"""\bottomrule
\end{{tabular}}
\begin{{tablenotes}}
\small
\item \textit{{Notes:}} $\mathbb{{E}}^{{\mathbb{{Q}}}}_t[\Delta r]=\sum_j q_{{j,t}}\,x_j$ with $x=(-50,-25,0,25,50)$ basis points and $q_{{j,t}}=p_{{j,t}}/\Sigma_t$; the open-ended states are valued at $\pm 50$ bps. Day 0 is the UTC calendar day of the release (FOMC statements at 14:00 ET, CPI and employment reports at 08:30 ET), which is contained in the Polymarket close of that day (00:00 UTC). Windows cumulate daily changes; the 3 April 2026 employment report (Good Friday) is assigned the first subsequent close, and the 14 April 2026 CPI release has no close within two days and is dropped. Panel B compares absolute daily changes on the $K$ release days with all other days ($T={len(d_eq)}$ daily changes, listing day excluded): ``Ann.'' and ``Other'' are means, VR is the ratio of variances, $p_{{BF}}$ the Brown--Forsythe test of equal variances and $p_{{MW}}$ the one-sided Mann--Whitney test that absolute changes are larger on release days. The last row compares weekend days with weekdays among non-release days. For comparison, the absolute daily change in the 10-year Treasury yield (CBOE TNX) is {f(ct['mean_abs_a'],1)} bps on release days versus {f(ct['mean_abs_b'],1)} bps otherwise (variance ratio {f(ct['var_ratio'],2)}, $p_{{BF}}={f(ct['bf_p'],3)}$).
\end{{tablenotes}}
\end{{threeparttable}}
\end{{table}}"""
save_table(tex, "table4_event_study")

# ----------------------------------------------------------------------------
# 7. Granger causality / VAR on trading-day grid (Tables 5, 6)
# ----------------------------------------------------------------------------
from statsmodels.tsa.api import VAR  # noqa: E402
import statsmodels.api as sm  # noqa: E402

grid = spy.index[(spy.index >= fed_main.index.min()) & (spy.index <= fed_main.index.max())]
eq_grid = eq.reindex(grid).ffill()
pnc_grid = fed_main["No change"].reindex(grid).ffill()
panel = pd.DataFrame({
    "Fed": eq_grid.diff(),                       # Δ E_Q[Δr], bps
    "TNX": (tnx.reindex(grid).diff() * 100),      # Δ 10-year yield, bps
    "SPY": np.log(spy.reindex(grid)).diff() * 100,
    "VIX": np.log(vix.reindex(grid)).diff() * 100,
    "BTC": np.log(btc.reindex(grid)).diff() * 100,
}).dropna()
STATS["granger"] = {"T": int(len(panel)), "first": str(panel.index.min().date()), "last": str(panel.index.max().date())}
sel = VAR(panel).select_order(maxlags=5)
STATS["granger"]["aic_lag"] = int(sel.aic); STATS["granger"]["bic_lag"] = int(sel.bic); STATS["granger"]["hqic_lag"] = int(sel.hqic)
p_main = max(int(sel.aic), 1)


def gc_pair(y: str, x: str, p: int):
    """Bivariate VAR(p) Granger F-test of x → y with HAC-robust Wald as well."""
    df = panel[[y, x]]
    r = VAR(df).fit(p)
    t = r.test_causality(y, [x], kind="f")
    # HAC-robust Wald on single equation
    Y = df[y].iloc[p:]
    X = pd.concat([df.shift(k).add_suffix(f"_L{k}") for k in range(1, p + 1)], axis=1).iloc[p:]
    X = sm.add_constant(X)
    m = sm.OLS(Y, X).fit(cov_type="HAC", cov_kwds={"maxlags": 5})
    R = np.zeros((p, X.shape[1]))
    for k in range(1, p + 1):
        R[k - 1, list(X.columns).index(f"{x}_L{k}")] = 1
    w = m.wald_test(R, use_f=True, scalar=True)
    return float(t.test_statistic), float(t.pvalue), float(w.statistic), float(w.pvalue)


names = {"Fed": r"$\Delta\mathbb{E}^{\mathbb{Q}}[\Delta r]$", "TNX": "10-year yield", "SPY": "SPY", "VIX": "VIX", "BTC": "BTC"}
g_rows, G = [], {}
pairs = [("Fed", v) for v in ["TNX", "SPY", "VIX", "BTC"]] + [(v, "Fed") for v in ["TNX", "SPY", "VIX", "BTC"]]
for x, y in pairs:
    res = {}
    for p in [p_main, 1, 2, 5]:
        res[p] = gc_pair(y, x, p)
    G[f"{x}->{y}"] = {str(k): v for k, v in res.items()}
    Fm, pm, Wm, pwm = res[p_main]
    g_rows.append(f"    {names[x]} $\\to$ {names[y]} & {f(Fm,2)}{stars(pm)} & {f(pm,3)} & {f(pwm,3)} & "
                  f"{f(res[1][1],3)} & {f(res[2][1],3)} & {f(res[5][1],3)} \\\\")
STATS["granger"]["pairs"] = G
# old-style (min-p over lags 1..5, misaligned returns) for the changelog
from statsmodels.tsa.stattools import grangercausalitytests  # noqa: E402
old = {}
for x, y in [("Fed", "VIX"), ("VIX", "Fed"), ("Fed", "SPY")]:
    r = grangercausalitytests(panel[[y, x]].values, maxlag=5, verbose=False)
    best = min(r, key=lambda k: r[k][0]["ssr_ftest"][1])
    old[f"{x}->{y}"] = {"minp_lag": int(best), "F": float(r[best][0]["ssr_ftest"][0]), "p": float(r[best][0]["ssr_ftest"][1])}
STATS["granger"]["min_p_over_lags_aligned"] = old

tex = rf"""\begin{{table}}[htbp]
\centering\small
\begin{{threeparttable}}
\caption{{Granger causality between the Fed state-price panel and U.S. asset returns}}
\label{{tab:granger}}
\setlength{{\tabcolsep}}{{4pt}}
\begin{{tabular}}{{l rrr rrr}}
\toprule
 & \multicolumn{{3}}{{c}}{{VAR($p^*$), $p^*={p_main}$ (AIC)}} & \multicolumn{{3}}{{c}}{{$p$-value at fixed lag}} \\
\cmidrule(lr){{2-4}}\cmidrule(lr){{5-7}}
Direction & $F$ & $p$ & $p$ (HAC) & $p=1$ & $p=2$ & $p=5$ \\
\midrule
{chr(10).join(g_rows)}
\bottomrule
\end{{tabular}}
\begin{{tablenotes}}
\small
\item \textit{{Notes:}} Bivariate VAR Granger tests on $T={len(panel)}$ U.S. trading days ({panel.index.min().strftime('%d %b %Y')} -- {panel.index.max().strftime('%d %b %Y')}). The Fed variable is the daily change in the risk-neutral expected June 2026 policy change (basis points), measured at the Polymarket close of the same calendar day (00:00 UTC, about four hours after the U.S. cash close); changes over weekends are assigned to Monday, as for the financial series. TNX is the daily change in the CBOE 10-year Treasury yield index (bps); SPY, VIX and BTC are daily log changes in percent. The lag order is selected by the AIC on the five-variable system (BIC selects {int(sel.bic)}, HQIC {int(sel.hqic)}). $F$ and $p$ are the conventional Granger $F$-test; $p$ (HAC) is the Wald test with Newey--West (5 lags) covariance. $^{{***}}$, $^{{**}}$, $^{{*}}$: 1\%, 5\%, 10\%.
\end{{tablenotes}}
\end{{threeparttable}}
\end{{table}}"""
save_table(tex, "table5_granger")

# VAR(1) coefficient table with HAC t-stats (equation by equation)
var_rows = []
Xv = sm.add_constant(panel.shift(1).iloc[1:])
coef, tval = {}, {}
for y in panel.columns:
    m = sm.OLS(panel[y].iloc[1:], Xv).fit(cov_type="HAC", cov_kwds={"maxlags": 5})
    coef[y] = m.params; tval[y] = m.tvalues
    STATS["granger"].setdefault("var1_r2", {})[y] = float(m.rsquared)
for reg in ["const"] + list(panel.columns):
    lab = "Constant" if reg == "const" else f"{names[reg]}$_{{t-1}}$"
    cells, tcells = [], []
    for y in panel.columns:
        c, t = coef[y][reg], tval[y][reg]
        p = 2 * (1 - stats.norm.cdf(abs(t)))
        cells.append(f"{f(c,3)}{stars(p)}")
        tcells.append(f"({f(t,2)})")
    var_rows.append(f"    {lab} & " + " & ".join(cells) + " \\\\")
    var_rows.append("     & " + " & ".join(tcells) + " \\\\")
r2 = " & ".join(f(STATS["granger"]["var1_r2"][y], 3) for y in panel.columns)
tex = rf"""\begin{{table}}[htbp]
\centering\small
\begin{{threeparttable}}
\caption{{VAR(1) estimates on the trading-day grid}}
\label{{tab:var}}
\setlength{{\tabcolsep}}{{4pt}}
\begin{{tabular}}{{l rrrrr}}
\toprule
 & \multicolumn{{5}}{{c}}{{Dependent variable (day $t$)}} \\
\cmidrule(lr){{2-6}}
Regressor & {names['Fed']} & 10-year yield & SPY & VIX & BTC \\
\midrule
{chr(10).join(var_rows)}
\midrule
$R^2$ & {r2} \\
Observations & \multicolumn{{5}}{{c}}{{{len(panel)-1}}} \\
\bottomrule
\end{{tabular}}
\begin{{tablenotes}}
\small
\item \textit{{Notes:}} Equation-by-equation OLS of each variable on one lag of all five variables (the AIC-selected order; see Table~\ref{{tab:granger}}). Newey--West (5 lags) $t$-statistics in parentheses. Variables are defined in Table~\ref{{tab:granger}}. $^{{***}}$, $^{{**}}$, $^{{*}}$: 1\%, 5\%, 10\%.
\end{{tablenotes}}
\end{{threeparttable}}
\end{{table}}"""
save_table(tex, "table6_var")

# contemporaneous correlations (for text)
STATS["granger"]["contemp_corr"] = panel.corr().round(3).to_dict()

# ----------------------------------------------------------------------------
# 8. Volatility: GARCH / GJR (QML) and HAR (appendix tables 7, 8)
# ----------------------------------------------------------------------------
from arch import arch_model  # noqa: E402

r_nc = np.log(fed_main["No change"]).diff().dropna() * 100  # percent log returns, calendar days
STATS["garch"] = {"T": int(len(r_nc)), "mean": float(r_nc.mean()), "sd": float(r_nc.std()), "min": float(r_nc.min()), "max": float(r_nc.max())}
# ARCH-LM
from statsmodels.stats.diagnostic import het_arch  # noqa: E402
lm = het_arch(r_nc - r_nc.mean(), nlags=5)
STATS["garch"]["arch_lm_stat"] = float(lm[0]); STATS["garch"]["arch_lm_p"] = float(lm[1])
lb = __import__("statsmodels.stats.diagnostic", fromlist=["acorr_ljungbox"]).acorr_ljungbox(r_nc ** 2, lags=[10], return_df=True)
STATS["garch"]["lb_sq_10_stat"] = float(lb["lb_stat"].iloc[0]); STATS["garch"]["lb_sq_10_p"] = float(lb["lb_pvalue"].iloc[0])
lb1 = __import__("statsmodels.stats.diagnostic", fromlist=["acorr_ljungbox"]).acorr_ljungbox(r_nc, lags=[10], return_df=True)
STATS["garch"]["lb_ret_10_stat"] = float(lb1["lb_stat"].iloc[0]); STATS["garch"]["lb_ret_10_p"] = float(lb1["lb_pvalue"].iloc[0])

garch_fits = {}
for label, o, dist in [("GARCH(1,1), normal", 0, "normal"), ("GJR-GARCH(1,1), normal", 1, "normal"),
                       ("GARCH(1,1), Student-$t$", 0, "t"), ("GJR-GARCH(1,1), Student-$t$", 1, "t")]:
    am = arch_model(r_nc, vol="GARCH", p=1, o=o, q=1, dist=dist, mean="Constant")
    res = am.fit(disp="off", cov_type="robust")
    garch_fits[label] = res
    STATS["garch"][label] = {"params": {k: float(v) for k, v in res.params.items()}, "se": {k: float(v) for k, v in res.std_err.items()},
                             "pvalues": {k: float(v) for k, v in res.pvalues.items()},
                             "loglik": float(res.loglikelihood), "aic": float(res.aic), "bic": float(res.bic),
                             "persistence": float(res.params.get("alpha[1]", 0) + res.params.get("beta[1]", 0) + 0.5 * res.params.get("gamma[1]", 0))}
    # Ljung-Box on squared standardised residuals
    z = res.std_resid.dropna()
    lbz = __import__("statsmodels.stats.diagnostic", fromlist=["acorr_ljungbox"]).acorr_ljungbox(z ** 2, lags=[10], return_df=True)
    STATS["garch"][label]["lb_z2_10_p"] = float(lbz["lb_pvalue"].iloc[0])
    STATS["garch"][label]["kurt_z"] = float(stats.kurtosis(z, fisher=False))
# LR test GARCH vs GJR (normal)
lr = 2 * (garch_fits["GJR-GARCH(1,1), normal"].loglikelihood - garch_fits["GARCH(1,1), normal"].loglikelihood)
STATS["garch"]["lr_gjr_vs_garch"] = float(lr); STATS["garch"]["lr_p"] = float(1 - stats.chi2.cdf(lr, 1))

param_order = ["mu", "omega", "alpha[1]", "gamma[1]", "beta[1]", "nu"]
pnames = {"mu": r"$\mu$", "omega": r"$\omega$", "alpha[1]": r"$\alpha_1$", "gamma[1]": r"$\gamma_1$", "beta[1]": r"$\beta_1$", "nu": r"$\nu$"}
g_rows = []
labels = list(garch_fits.keys())
for pn in param_order:
    cells, secells = [], []
    for lab in labels:
        res = garch_fits[lab]
        if pn in res.params:
            cells.append(f"{f(res.params[pn],3)}{stars(res.pvalues[pn])}"); secells.append(f"({f(res.std_err[pn],3)})")
        else:
            cells.append("--"); secells.append("")
    g_rows.append(f"    {pnames[pn]} & " + " & ".join(cells) + " \\\\")
    g_rows.append("     & " + " & ".join(secells) + " \\\\")
ll = " & ".join(f(garch_fits[l].loglikelihood, 1) for l in labels)
aic = " & ".join(f(garch_fits[l].aic, 1) for l in labels)
bic = " & ".join(f(garch_fits[l].bic, 1) for l in labels)
pers = " & ".join(f(STATS["garch"][l]["persistence"], 3) for l in labels)
lbz = " & ".join(f(STATS["garch"][l]["lb_z2_10_p"], 3) for l in labels)
tex = rf"""\begin{{table}}[htbp]
\centering\small
\begin{{threeparttable}}
\caption{{GARCH-family estimates for the ``No change'' contract (daily log returns, percent)}}
\label{{tab:garch}}
\setlength{{\tabcolsep}}{{4pt}}
\begin{{tabular}}{{l rrrr}}
\toprule
 & GARCH & GJR-GARCH & GARCH & GJR-GARCH \\
Parameter & normal & normal & Student-$t$ & Student-$t$ \\
\midrule
{chr(10).join(g_rows)}
\midrule
Log-likelihood & {ll} \\
AIC & {aic} \\
BIC & {bic} \\
Persistence $\alpha_1+\beta_1+\gamma_1/2$ & {pers} \\
Ljung--Box $p$, $z_t^2$ (10 lags) & {lbz} \\
Observations & \multicolumn{{4}}{{c}}{{{len(r_nc)}}} \\
\bottomrule
\end{{tabular}}
\begin{{tablenotes}}
\small
\item \textit{{Notes:}} Constant-mean models estimated by (quasi-)maximum likelihood on calendar-day log returns of the ``No change'' contract, listing day excluded; Bollerslev--Wooldridge robust standard errors in parentheses. The GJR variance equation is $\sigma_t^2=\omega+\alpha_1\varepsilon_{{t-1}}^2+\gamma_1\varepsilon_{{t-1}}^2\mathbf{{1}}\{{\varepsilon_{{t-1}}<0\}}+\beta_1\sigma_{{t-1}}^2$. The likelihood-ratio statistic for $\gamma_1=0$ (normal innovations) is {f(lr,2)} ($p={f(STATS['garch']['lr_p'],3)}$). Engle's ARCH-LM test (5 lags) on the returns: {f(STATS['garch']['arch_lm_stat'],2)} ($p={f(STATS['garch']['arch_lm_p'],3)}$). $^{{***}}$, $^{{**}}$, $^{{*}}$: 1\%, 5\%, 10\%.
\end{{tablenotes}}
\end{{threeparttable}}
\end{{table}}"""
save_table(tex, "table7_garch")

# HAR with standard component definitions
rv_d = r_nc.abs()
har = pd.DataFrame({"rv_d": rv_d, "rv_w": rv_d.rolling(5).mean(), "rv_m": rv_d.rolling(22).mean()})
y = rv_d.shift(-1)
dat = pd.concat([y.rename("y"), har], axis=1).dropna()
hm = sm.OLS(dat["y"], sm.add_constant(dat[["rv_d", "rv_w", "rv_m"]])).fit(cov_type="HAC", cov_kwds={"maxlags": 5})
# squared-return version
rv2 = r_nc ** 2
har2 = pd.DataFrame({"rv_d": rv2, "rv_w": rv2.rolling(5).mean(), "rv_m": rv2.rolling(22).mean()})
dat2 = pd.concat([rv2.shift(-1).rename("y"), har2], axis=1).dropna()
hm2 = sm.OLS(dat2["y"], sm.add_constant(dat2[["rv_d", "rv_w", "rv_m"]])).fit(cov_type="HAC", cov_kwds={"maxlags": 5})
STATS["har"] = {"abs": {"params": hm.params.to_dict(), "t": hm.tvalues.to_dict(), "r2": float(hm.rsquared), "adj_r2": float(hm.rsquared_adj), "n": int(hm.nobs),
                        "F_p": float(hm.f_pvalue)},
                "sq": {"params": hm2.params.to_dict(), "t": hm2.tvalues.to_dict(), "r2": float(hm2.rsquared), "adj_r2": float(hm2.rsquared_adj), "n": int(hm2.nobs),
                       "F_p": float(hm2.f_pvalue)}}
hn = {"const": "Constant", "rv_d": r"$RV^{(d)}_t$", "rv_w": r"$RV^{(w)}_t$", "rv_m": r"$RV^{(m)}_t$"}
h_rows = []
for k in ["const", "rv_d", "rv_w", "rv_m"]:
    h_rows.append(f"    {hn[k]} & {f(hm.params[k],3)}{stars(hm.pvalues[k])} & ({f(hm.bse[k],3)}) & {f(hm2.params[k],3)}{stars(hm2.pvalues[k])} & ({f(hm2.bse[k],3)}) \\\\")
tex = rf"""\begin{{table}}[htbp]
\centering\small
\begin{{threeparttable}}
\caption{{HAR regressions for the ``No change'' contract}}
\label{{tab:har_rv}}
\begin{{tabular}}{{l rr rr}}
\toprule
 & \multicolumn{{2}}{{c}}{{$RV_t=|r_t|$}} & \multicolumn{{2}}{{c}}{{$RV_t=r_t^2$}} \\
\cmidrule(lr){{2-3}}\cmidrule(lr){{4-5}}
Regressor & Coef. & (SE) & Coef. & (SE) \\
\midrule
{chr(10).join(h_rows)}
\midrule
$R^2$ & \multicolumn{{2}}{{c}}{{{f(hm.rsquared,3)}}} & \multicolumn{{2}}{{c}}{{{f(hm2.rsquared,3)}}} \\
Adjusted $R^2$ & \multicolumn{{2}}{{c}}{{{f(hm.rsquared_adj,3)}}} & \multicolumn{{2}}{{c}}{{{f(hm2.rsquared_adj,3)}}} \\
Joint $F$-test $p$-value & \multicolumn{{2}}{{c}}{{{f(hm.f_pvalue,3)}}} & \multicolumn{{2}}{{c}}{{{f(hm2.f_pvalue,3)}}} \\
Observations & \multicolumn{{2}}{{c}}{{{int(hm.nobs)}}} & \multicolumn{{2}}{{c}}{{{int(hm2.nobs)}}} \\
\bottomrule
\end{{tabular}}
\begin{{tablenotes}}
\small
\item \textit{{Notes:}} $RV^{{(d)}}_{{t+1}}=c+\beta_d RV^{{(d)}}_t+\beta_w RV^{{(w)}}_t+\beta_m RV^{{(m)}}_t+\varepsilon_{{t+1}}$ with $RV^{{(w)}}_t$ and $RV^{{(m)}}_t$ the averages of the daily measure over the past 5 and 22 calendar days \citep{{corsi2009simple}}. Daily log returns in percent; listing day excluded. Newey--West (5 lags) standard errors. $^{{***}}$, $^{{**}}$, $^{{*}}$: 1\%, 5\%, 10\%.
\end{{tablenotes}}
\end{{threeparttable}}
\end{{table}}"""
save_table(tex, "table8_har_rv")

# ----------------------------------------------------------------------------
# 9. Roll spread
# ----------------------------------------------------------------------------
def roll_spread(p: pd.Series, window: int = 20) -> pd.Series:
    dp = p.diff()
    cov = dp.rolling(window).apply(lambda x: np.cov(x[:-1], x[1:])[0, 1] if len(x) > 2 else np.nan, raw=True)
    return 2 * np.sqrt(np.maximum(-cov, 0))


roll = {c: roll_spread(fed_main[c]) for c in STATES}
rs = roll["No change"].dropna()
STATS["roll"] = {"nochange_mean": float(rs.mean()), "nochange_median": float(rs.median()), "nochange_first30_mean": float(rs.iloc[:30].mean()),
                 "nochange_last30_mean": float(rs.iloc[-30:].mean()), "share_zero": float((rs == 0).mean()),
                 "autocov_sign_share_negative": float((fed_main["No change"].diff().rolling(20).apply(lambda x: np.cov(x[:-1], x[1:])[0, 1], raw=True).dropna() < 0).mean())}

# ----------------------------------------------------------------------------
# 10. Election panel statistics
# ----------------------------------------------------------------------------
tr, ha, bi = election["Donald Trump"], election.get("Kamala Harris"), election.get("Joe Biden")
STATS["election"] = {"trump_last": float(tr.dropna().iloc[-1]), "trump_max": float(tr.max()), "trump_min": float(tr.min()),
                     "trump_mean": float(tr.mean()), "harris_last": float(ha.dropna().iloc[-1]), "harris_max": float(ha.max()),
                     "biden_before_dropout": float(bi.loc[:pd.Timestamp("2024-07-20")].iloc[-1]), "biden_after_dropout": float(bi.loc[pd.Timestamp("2024-07-22"):].iloc[0]),
                     "harris_before_dropout": float(ha.loc[:pd.Timestamp("2024-07-20")].iloc[-1]), "harris_after_dropout": float(ha.loc[pd.Timestamp("2024-07-22"):].iloc[0]),
                     "trump_jul12": float(tr.loc[pd.Timestamp("2024-07-12")]), "trump_jul14": float(tr.loc[pd.Timestamp("2024-07-14")]),
                     "trump_sep9": float(tr.loc[pd.Timestamp("2024-09-09")]), "trump_sep11": float(tr.loc[pd.Timestamp("2024-09-11")]),
                     "trump_oct_peak": float(tr.loc["2024-10"].max()), "trump_oct_peak_date": str(tr.loc["2024-10"].idxmax().date()),
                     "trump_nov5": float(tr.loc[pd.Timestamp("2024-11-05")]), "harris_nov5": float(ha.loc[pd.Timestamp("2024-11-05")]),
                     "others_sum_nov5": float(el_complete.loc[pd.Timestamp("2024-11-05")].drop(["Donald Trump", "Kamala Harris"]).sum()),
                     "n_contracts": int(election.shape[1]), "T": int(len(election)), "T_complete": int(len(el_complete))}

# ----------------------------------------------------------------------------
# 11. Figures
# ----------------------------------------------------------------------------
# fig1 state prices
fig, ax = plt.subplots(figsize=(6.5, 3.6))
for i, c in enumerate(STATES):
    ax.plot(fed_main.index, fed_main[c], color=PAL[i], ls=LS[i], label=c)
for d0 in fomc:
    ax.axvline(d0, color="0.6", lw=0.6, ls=":")
ax.set_ylabel("State price (USDC per \\$1 payoff)"); ax.set_ylim(-0.02, 1.02)
ax.legend(ncol=3, loc="center left", frameon=False)
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
save_fig(fig, "fig1_state_prices")

# fig2 Σ with bounds
fig, ax = plt.subplots(figsize=(6.5, 3.4))
ax.plot(sigma.index, sigma.values, color=PAL[0], lw=1.2, label=r"$\Sigma_t$ (sum of closing prices)")
ax.plot(sigma.index, lower_fed, color=PAL[1], lw=1.0, ls="--", label=r"Lower bound $e^{-r_t\tau_t}$ (DFF)")
ax.axhline(1.0, color="black", lw=0.8, label="Upper bound 1")
ax.fill_between(sigma.index, lower_fed - 0.005, 1.005, color="0.85", alpha=0.5, label=r"Bounds $\pm 0.5$ pp")
ax.set_ylim(0.95, 1.12); ax.set_ylabel(r"$\Sigma_t$")
ax.legend(loc="upper right", frameon=False, ncol=2)
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
save_fig(fig, "fig2_completeness_test")

# fig3 RN density snapshots
fig, ax = plt.subplots(figsize=(6.5, 3.4))
idx = q.index
snap = [idx[0], idx[len(idx) // 4], idx[len(idx) // 2], idx[3 * len(idx) // 4], idx[-1]]
w = 3.5
for i, d0 in enumerate(snap):
    ax.bar(STATE_BPS + (i - 2) * w, q.loc[d0].values, width=w, color=PAL[i], label=d0.strftime("%d %b %Y"), edgecolor="white", lw=0.5)
ax.set_xticks(STATE_BPS); ax.set_xticklabels(["$\\leq-50$", "$-25$", "0", "$+25$", "$\\geq+50$"])
ax.set_xlabel("Change in the target range upper bound, June 2026 meeting (bps)"); ax.set_ylabel("Risk-neutral probability $q_{j,t}$")
ax.legend(frameon=False)
save_fig(fig, "fig3_risk_neutral_density")

# fig4 E_Q path with release markers
fig, axes = plt.subplots(2, 1, figsize=(6.5, 5.2), sharex=True, gridspec_kw={"height_ratios": [2, 1.3]})
axes[0].plot(eq.index, eq.values, color=PAL[0], lw=1.2)
for lab, ds, mk, col in [("FOMC", fomc, "v", PAL[1]), ("CPI", cpi, "o", PAL[2]), ("Employment report", nfp, "s", PAL[3])]:
    axes[0].scatter(ds, eq.reindex(ds).values, marker=mk, s=28, color=col, zorder=3, label=lab)
axes[0].set_ylabel(r"$\mathbb{E}^{\mathbb{Q}}_t[\Delta r]$ (bps)"); axes[0].legend(frameon=False, loc="lower left")
axes[1].bar(d_eq.index, d_eq.values, width=1, color=np.where(d_eq.index.isin(all_ann), PAL[1], "0.6"))
axes[1].set_ylabel(r"$\Delta\mathbb{E}^{\mathbb{Q}}_t[\Delta r]$ (bps)")
axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
save_fig(fig, "fig4_expected_policy_change")

# fig5 announcement-day distributions
fig, axes = plt.subplots(1, 2, figsize=(6.5, 3.0))
groups = [("Other", d_eq[~is_ann]), ("Employment", d_eq[d_eq.index.isin(nfp)]), ("CPI", d_eq[d_eq.index.isin(cpi)]), ("FOMC", d_eq[d_eq.index.isin(fomc)])]
for k, (lab, s) in enumerate(groups):
    jitter = (np.random.RandomState(k).rand(len(s)) - 0.5) * 0.3
    axes[0].scatter(np.full(len(s), k) + jitter, s.abs().values, s=10, color=PAL[k], alpha=0.8)
    axes[0].hlines(s.abs().mean(), k - 0.3, k + 0.3, color="black", lw=1.2)
axes[0].set_xticks(range(4)); axes[0].set_xticklabels([g[0] for g in groups]); axes[0].set_ylabel(r"$|\Delta\mathbb{E}^{\mathbb{Q}}_t[\Delta r]|$ (bps)")
axes[0].set_title("A. Expected policy change", loc="left")
groups2 = [("Other", d_pnc[~is_ann]), ("Employment", d_pnc[d_pnc.index.isin(nfp)]), ("CPI", d_pnc[d_pnc.index.isin(cpi)]), ("FOMC", d_pnc[d_pnc.index.isin(fomc)])]
for k, (lab, s) in enumerate(groups2):
    jitter = (np.random.RandomState(k).rand(len(s)) - 0.5) * 0.3
    axes[1].scatter(np.full(len(s), k) + jitter, s.abs().values, s=10, color=PAL[k], alpha=0.8)
    axes[1].hlines(s.abs().mean(), k - 0.3, k + 0.3, color="black", lw=1.2)
axes[1].set_xticks(range(4)); axes[1].set_xticklabels([g[0] for g in groups2]); axes[1].set_ylabel(r"$|\Delta p_{\mathrm{No\ change}}|$ (pp)")
axes[1].set_title("B. ``No change'' state price", loc="left")
fig.tight_layout()
save_fig(fig, "fig5_event_study_cars")

# fig7 Granger p-value heatmap (AIC lag)
mat = pd.DataFrame(np.nan, index=panel.columns, columns=panel.columns)
for x in panel.columns:
    for y in panel.columns:
        if x != y:
            mat.loc[y, x] = gc_pair(y, x, p_main)[1]
fig, ax = plt.subplots(figsize=(4.6, 3.8))
im = ax.imshow(mat.values.astype(float), cmap="Blues_r", vmin=0, vmax=0.5)
ax.set_xticks(range(5)); ax.set_xticklabels([names[c].replace("$\\Delta\\mathbb{E}^{\\mathbb{Q}}[\\Delta r]$", "Fed") for c in panel.columns], rotation=30, ha="right")
ax.set_yticks(range(5)); ax.set_yticklabels([names[c].replace("$\\Delta\\mathbb{E}^{\\mathbb{Q}}[\\Delta r]$", "Fed") for c in panel.columns])
ax.set_xlabel("Causing variable (lagged)"); ax.set_ylabel("Caused variable"); ax.grid(False)
for i in range(5):
    for j in range(5):
        v = mat.values[i, j]
        if not np.isnan(v):
            ax.text(j, i, f"{v:.2f}", ha="center", va="center", fontsize=8, color="white" if v < 0.2 else "black")
fig.colorbar(im, ax=ax, label=f"Granger $p$-value, VAR({p_main})", shrink=0.8)
save_fig(fig, "fig7_granger_heatmap")
STATS["granger"]["full_matrix_p"] = mat.astype(float).round(4).to_dict()

# fig9 volatility
fig, axes = plt.subplots(2, 1, figsize=(6.5, 4.6), sharex=True)
axes[0].bar(r_nc.index, r_nc.values, width=1, color=PAL[0]); axes[0].set_ylabel("Log return (%)"); axes[0].set_title("A. Daily log return, ``No change''", loc="left")
for lab, col, ls_ in [("GARCH(1,1), normal", PAL[0], "-"), ("GJR-GARCH(1,1), Student-$t$", PAL[1], "--")]:
    cv = garch_fits[lab].conditional_volatility
    axes[1].plot(cv.index, cv.values, color=col, ls=ls_, label=lab)
axes[1].set_ylabel("Conditional SD (%)"); axes[1].legend(frameon=False); axes[1].set_title("B. Conditional volatility", loc="left")
axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
save_fig(fig, "fig9_volatility")

# fig10 liquidity: Roll spread + order-book depth
fig, axes = plt.subplots(1, 2, figsize=(6.5, 3.2))
axes[0].plot(roll["No change"].index, roll["No change"].values * 100, color=PAL[0], label="No change")
axes[0].plot(roll["25 bps decrease"].index, roll["25 bps decrease"].values * 100, color=PAL[1], ls="--", label="25 bps decrease")
axes[0].set_ylabel("Roll spread estimate (cents)"); axes[0].legend(frameon=False); axes[0].set_title("A. Roll (1984) estimator, 20-day window", loc="left")
axes[0].xaxis.set_major_formatter(mdates.DateFormatter("%b"))
for i, c in enumerate(STATES):
    path = DATA_RAW_POLY / f"order_book_{[m for m in fed_mk if m['name']==c][0]['yes_token'][:16]}.json"
    d = json.load(open(path))
    bids = sorted([(float(b["price"]), float(b["size"])) for b in d["bids"]], reverse=True)
    asks = sorted([(float(a["price"]), float(a["size"])) for a in d["asks"]])
    mid = ob[c]["mid"]
    xb = [(mid - p) * 100 for p, _ in bids]; yb = np.cumsum([p * s for p, s in bids]) / 1e3
    xa = [(p - mid) * 100 for p, _ in asks]; ya = np.cumsum([p * s for p, s in asks]) / 1e3
    axes[1].step(xb, yb, where="post", color=PAL[i], ls="--", lw=1)
    axes[1].step(xa, ya, where="post", color=PAL[i], ls="-", lw=1, label=c)
axes[1].set_xscale("log"); axes[1].set_xlim(0.05, 100); axes[1].set_yscale("log")
axes[1].set_xlabel("Distance from midpoint (cents)"); axes[1].set_ylabel("Cumulative depth (USD 000)")
axes[1].legend(frameon=False, fontsize=6, loc="lower right"); axes[1].set_title("B. Order book, 22 May 2026 (ask solid, bid dashed)", loc="left")
fig.tight_layout()
save_fig(fig, "fig10_liquidity")

# fig11 election
fig, axes = plt.subplots(2, 1, figsize=(6.5, 5.0), sharex=True, gridspec_kw={"height_ratios": [2, 1]})
for i, c in enumerate(["Donald Trump", "Kamala Harris", "Joe Biden", "Michelle Obama", "Robert F. Kennedy Jr."]):
    axes[0].plot(election.index, election[c], color=PAL[i], ls=LS[i], label=c)
for lab, d0 in ELECTION_KEY_DATES.items():
    axes[0].axvline(pd.Timestamp(d0), color="0.6", lw=0.6, ls=":")
    axes[0].text(pd.Timestamp(d0), 0.98, lab, rotation=90, fontsize=6, va="top", ha="right", color="0.3")
axes[0].set_ylabel("Price (USDC)"); axes[0].legend(frameon=False, ncol=2, loc="upper left", fontsize=7); axes[0].set_ylim(0, 1.0)
axes[1].plot(sigma_el.index, sigma_el.values, color=PAL[0]); axes[1].axhline(1, color="black", lw=0.8)
axes[1].plot(sigma_el.index, lower_el, color=PAL[1], ls="--", lw=0.9)
axes[1].set_ylabel(r"$\Sigma_t$, 17 contracts"); axes[1].set_ylim(0.9, 1.1)
axes[1].xaxis.set_major_locator(mdates.MonthLocator()); axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%b"))
axes[1].set_xlabel("2024")
save_fig(fig, "fig11_election_case_study")

# fig12 deviation distribution
fig, axes = plt.subplots(1, 2, figsize=(6.5, 2.9))
dev = sigma - 1
axes[0].hist(dev.values * 100, bins=30, color=PAL[0], edgecolor="white"); axes[0].axvline(0, color="black", lw=0.8)
axes[0].axvline(0.5, color=PAL[1], ls="--", lw=0.8); axes[0].set_xlabel(r"$\Sigma_t-1$ (pp), Fed group"); axes[0].set_ylabel("Days")
axes[0].set_title("A. Fed June 2026", loc="left")
dev_el = sigma_el - 1
axes[1].hist(dev_el.values * 100, bins=30, color=PAL[0], edgecolor="white"); axes[1].axvline(0, color="black", lw=0.8)
axes[1].set_xlabel(r"$\Sigma_t-1$ (pp), election group"); axes[1].set_title("B. Election 2024", loc="left")
fig.tight_layout()
save_fig(fig, "fig12_arbitrage_violations")


# ----------------------------------------------------------------------------
# 11b. Sensitivity checks reported in Section 7
# ----------------------------------------------------------------------------
SENS = {}
# (a) announcement tests excluding the two largest FOMC days
big2 = d_pnc.abs().nlargest(2).index
mask_ex = ~d_pnc.index.isin(big2)
c_ex = compare(d_pnc[mask_ex], is_ann[mask_ex], "All releases ex 2 largest FOMC days")
c_ex_eq = compare(d_eq[mask_ex], is_ann[mask_ex], "All releases ex 2 largest FOMC days")
SENS["ex_two_largest_days"] = {"days": [str(x.date()) for x in big2], "pnc": c_ex, "eq": c_ex_eq}
# (b) cost allowance of 1 cent on the complete set
tau_f = (RESOLUTION_FED - sigma.index).days.values / 365.0
r_f = dff.reindex(sigma.index).ffill().bfill().values / 100.0
low_f = np.exp(-r_f * tau_f)
SENS["share_inside_band_1c"] = float(((sigma.values >= low_f - 0.01) & (sigma.values <= 1.01)).mean())
SENS["share_inside_band_roll"] = float(((sigma.values >= low_f - roll["No change"].reindex(sigma.index).fillna(0).values) & (sigma.values <= 1 + roll["No change"].reindex(sigma.index).fillna(0).values)).mean())
# (f) lower bound with the target-range upper bound (3.75%) instead of DFF
low_ub = np.exp(-0.0375 * tau_f)
SENS["max_abs_diff_lower_bound_ub_vs_dff"] = float(np.abs(low_ub - low_f).max())
# (e) ±75 valuation of open-ended states
eq75 = pd.Series(q.values @ np.array([-75.0, -25.0, 0.0, 25.0, 75.0]), index=q.index)
d_eq75 = eq75.diff().dropna()
SENS["eq75_corr_with_eq50"] = float(np.corrcoef(d_eq75.values, d_eq.values)[0, 1])
SENS["eq75_sd_ratio"] = float(d_eq75.std() / d_eq.std())
SENS["eq75_all_releases"] = compare(d_eq75, is_ann, "All releases, ±75")
# (c) Granger with Δp(No change) as the Polymarket variable
panel_pnc = panel.copy(); panel_pnc["Fed"] = (pnc_grid.diff() * 100).reindex(panel.index)
panel_pnc = panel_pnc.dropna()
def gc_pair_df(df, y, x, p):
    r = VAR(df[[y, x]]).fit(p); t = r.test_causality(y, [x], kind="f"); return float(t.test_statistic), float(t.pvalue)
SENS["granger_pnc"] = {f"{x}->{y}": {str(p): gc_pair_df(panel_pnc, y, x, p) for p in [1, 2]} for x, y in pairs}
# (d) half-sample Granger at p=1
h = len(panel) // 2
SENS["granger_halves_p1"] = {}
for lab, df in [("first_half", panel.iloc[:h]), ("second_half", panel.iloc[h:])]:
    SENS["granger_halves_p1"][lab] = {f"{x}->{y}": gc_pair_df(df, y, x, 1) for x, y in pairs}
    SENS["granger_halves_p1"][lab]["T"] = int(len(df))
# (g) election facts
SENS["trump_max_date"] = str(tr.idxmax().date()); SENS["trump_max"] = float(tr.max())
SENS["nochange_mar17"] = float(fed_main["No change"].loc["2026-03-17"]); SENS["nochange_mar18"] = float(fed_main["No change"].loc["2026-03-18"])
SENS["cut25_mar17"] = float(fed_main["25 bps decrease"].loc["2026-03-17"]); SENS["cut25_mar18"] = float(fed_main["25 bps decrease"].loc["2026-03-18"])
SENS["n_release_days"] = int(is_ann.sum()); SENS["n_other_days"] = int((~is_ann).sum())
STATS["sensitivity"] = SENS

# ----------------------------------------------------------------------------
# 12. Save stats
# ----------------------------------------------------------------------------
def _clean(o):
    if isinstance(o, dict):
        return {str(k): _clean(v) for k, v in o.items()}
    if isinstance(o, (list, tuple)):
        return [_clean(v) for v in o]
    if isinstance(o, (np.floating, float)):
        return None if np.isnan(o) else float(o)
    if isinstance(o, (np.integer,)):
        return int(o)
    if isinstance(o, (pd.Timestamp,)):
        return str(o)
    return o


(ROOT / "results" / "revision_stats.json").write_text(json.dumps(_clean(STATS), indent=2))
fed_main.assign(state_price_sum=sigma, EQ_bps=eq).to_parquet(ROOT / "data" / "processed" / "fed_june2026_state_prices_aligned.parquet")
election.to_parquet(ROOT / "data" / "processed" / "election_2024_panel_aligned.parquet")
panel.to_parquet(ROOT / "data" / "processed" / "trading_day_panel.parquet")
print(json.dumps(_clean({k: STATS[k] for k in ["fed", "completeness", "orderbook", "announcement", "granger", "garch", "har", "roll", "election", "events"]}), indent=1)[:20000])