Python 56.2%
TeX 43.7%
1"""Revision analysis (editorial-board review, 2026-10-17) for UQO WP 8.23Re-derives every table and figure of the revised paper from the *local* raw data4(``data/raw/``) without any network call. Compared with ``run_analysis.py`` the5script67* labels each Polymarket daily observation by the UTC day it closes (the CLOB8 ``prices-history`` end-point stamps the daily point at 00:00 UTC of the *next*9 day), so that Polymarket closes, U.S. market closes and announcement dates share10 a common calendar;11* flags the listing day of the Fed group (10 Dec 2025, seeded quotes) and reports12 the completeness statistics with and without it;13* tests the no-arbitrage restriction on the price of the complete set against14 discount-adjusted bounds built from the effective federal funds rate (FRED DFF);15* replaces the per-contract CAR event study (2–3 events) by an announcement-day16 analysis of the risk-neutral expected policy change (14 scheduled releases);17* re-runs Granger tests on a trading-day grid with information-criterion lag18 selection (the original code reported the minimum p-value over 1–5 lags) and19 adds the 10-year Treasury yield;20* re-estimates GARCH/GJR/HAR models with the listing day excluded and QML errors;21* generates descriptive, order-book and election tables from the raw metadata.2223Usage: python3 scripts/02_revision_analysis.py (≈ 1 min, offline)24Outputs: results/tables/*.tex, results/revision_stats.json, figures/fig*.pdf|png25"""2627from __future__ import annotations2829import json30import sys31import warnings32from datetime import date33from pathlib import Path3435import numpy as np36import pandas as pd37from scipy import stats3839ROOT = Path(__file__).resolve().parents[1]40sys.path.insert(0, str(ROOT))4142from src.config import ( # noqa: E40243 CPI_RELEASE_DATES,44 DATA_RAW_FRED,45 DATA_RAW_POLY,46 DATA_RAW_YAHOO,47 ELECTION_2024_EVENT_ID,48 FED_JUNE_2026_EVENT_ID,49 FIGURES,50 FOMC_DATES,51 NFP_RELEASE_DATES,52 TABLES,53 ELECTION_KEY_DATES,54)5556warnings.filterwarnings("ignore")5758import matplotlib # noqa: E4025960matplotlib.use("Agg")61import matplotlib.dates as mdates # noqa: E40262import matplotlib.pyplot as plt # noqa: E4026364# Okabe–Ito palette (CVD-safe, validated), fixed order per state65PAL = ["#0072B2", "#D55E00", "#009E73", "#CC79A7", "#56B4E9", "#E69F00", "#000000"]66LS = ["-", "--", "-", "-.", ":", "-", "--"]67plt.rcParams.update({68 "font.family": "serif", "font.serif": ["Times New Roman", "DejaVu Serif"],69 "font.size": 9, "axes.titlesize": 10, "axes.labelsize": 9, "legend.fontsize": 8,70 "xtick.labelsize": 8, "ytick.labelsize": 8, "axes.spines.top": False,71 "axes.spines.right": False, "axes.grid": True, "grid.alpha": 0.25,72 "grid.linestyle": "--", "lines.linewidth": 1.3, "figure.dpi": 150,73})7475STATES = ["50+ bps decrease", "25 bps decrease", "No change", "25 bps increase", "50+ bps increase"]76STATE_BPS = np.array([-50.0, -25.0, 0.0, 25.0, 50.0])77RESOLUTION_FED = pd.Timestamp("2026-06-17")78RESOLUTION_ELECTION = pd.Timestamp("2024-11-05")79TICK = 0.00180STATS: dict = {}818283# ----------------------------------------------------------------------------84# helpers85# ----------------------------------------------------------------------------86def save_fig(fig, name: str) -> None:87 FIGURES.mkdir(parents=True, exist_ok=True)88 for ext in ("pdf", "png"):89 fig.savefig(FIGURES / f"{name}.{ext}", bbox_inches="tight", dpi=300)90 plt.close(fig)919293def save_table(tex: str, name: str) -> None:94 TABLES.mkdir(parents=True, exist_ok=True)95 (TABLES / f"{name}.tex").write_text(tex)969798def esc(s: str) -> str:99 return (str(s).replace("&", r"\&").replace("%", r"\%").replace("$", r"\$")100 .replace("_", r"\_").replace("#", r"\#"))101102103def f(x, d=3, blank="--"):104 if x is None or (isinstance(x, float) and np.isnan(x)):105 return blank106 return f"{x:,.{d}f}"107108109def stars(p):110 if p is None or np.isnan(p):111 return ""112 return "$^{***}$" if p < 0.01 else "$^{**}$" if p < 0.05 else "$^{*}$" if p < 0.10 else ""113114115def nw_mean_test(x: pd.Series, lags: int = 5, mu0: float = 0.0):116 """HAC (Newey–West) t-test of H0: E[x] = mu0."""117 import statsmodels.api as sm118 x = x.dropna()119 mod = sm.OLS(x.values - mu0, np.ones((len(x), 1))).fit(cov_type="HAC", cov_kwds={"maxlags": lags})120 return float(mod.params[0] + mu0), float(mod.bse[0]), float(mod.tvalues[0]), float(mod.pvalues[0])121122123def load_history(token: str) -> pd.Series:124 """Daily close labelled by the UTC day it closes (point stamped 00:00 UTC of day t+1 → day t)."""125 path = DATA_RAW_POLY / f"price_history_{token[:16]}.json"126 raw = json.load(open(path)) if path.exists() else []127 if not raw:128 return pd.Series(dtype=float)129 df = pd.DataFrame(raw)130 ts = pd.to_datetime(df["t"], unit="s", utc=True)131 df["day"] = (ts - pd.Timedelta(minutes=5)).dt.normalize().dt.tz_localize(None)132 df["p"] = df["p"].astype(float)133 return df.drop_duplicates("day", keep="last").set_index("day")["p"].sort_index()134135136def load_event(event_id: int) -> tuple[dict, list[dict]]:137 ev = json.load(open(DATA_RAW_POLY / f"event_{event_id}.json"))138 out = []139 for m in ev["markets"]:140 toks = json.loads(m.get("clobTokenIds", "[]"))141 out.append({"id": m["id"], "name": m.get("groupItemTitle", m["question"][:40]),142 "question": m["question"], "yes_token": toks[0] if toks else "",143 "volume": float(m.get("volumeNum") or 0), "liquidity": float(m.get("liquidityNum") or 0),144 "closed": m.get("closed", False), "start": m.get("startDate", ""),145 "end": m.get("endDate", ""), "spread": m.get("spread"),146 "best_bid": m.get("bestBid"), "best_ask": m.get("bestAsk")})147 return ev, out148149150def build_panel(markets: list[dict], order: list[str] | None = None) -> pd.DataFrame:151 cols = {m["name"]: load_history(m["yes_token"]) for m in markets if m["yes_token"]}152 cols = {k: v for k, v in cols.items() if len(v)}153 panel = pd.DataFrame(cols).sort_index()154 if order:155 panel = panel[[c for c in order if c in panel.columns]]156 return panel.ffill(limit=3)157158159# ----------------------------------------------------------------------------160# 1. Data161# ----------------------------------------------------------------------------162fed_ev, fed_mk = load_event(FED_JUNE_2026_EVENT_ID)163fed = build_panel(fed_mk, STATES)164el_ev, el_mk = load_event(ELECTION_2024_EVENT_ID)165el_order = sorted(el_mk, key=lambda m: -m["volume"])166election = build_panel(el_order, [m["name"] for m in el_order])167168dff = pd.read_csv(DATA_RAW_FRED / "fed_funds_rate.csv", parse_dates=["date"], index_col="date")["DFF"]169dff = dff.asfreq("D").ffill()170dff = dff.reindex(pd.date_range(dff.index.min(), "2026-06-30")).ffill()171172173def yahoo(ticker: str) -> pd.Series:174 df = pd.read_csv(DATA_RAW_YAHOO / f"{ticker}.csv", parse_dates=["date"], index_col="date")175 df.index = pd.to_datetime(df.index, utc=True).normalize().tz_localize(None)176 return df["Close"].astype(float)177178179spy, vix, tnx, btc = yahoo("SPY"), yahoo("VIX"), yahoo("TNX"), yahoo("BTC_USD")180181listing_day = fed.index.min()182fed_complete = fed.dropna()183fed_main = fed_complete[fed_complete.index > listing_day] # listing day excluded184sigma_all = fed_complete.sum(axis=1)185sigma = fed_main.sum(axis=1)186STATS["fed"] = {187 "first_day": str(fed.index.min().date()), "last_day": str(fed.index.max().date()),188 "n_rows_union": int(len(fed)), "n_complete": int(len(fed_complete)), "n_main": int(len(fed_main)),189 "listing_sigma": float(sigma_all.iloc[0]), "listing_prices": fed_complete.iloc[0].round(4).to_dict(),190 "n_weekend_days": int((fed_main.index.dayofweek >= 5).sum()),191 "calendar_span_days": int((fed.index.max() - fed.index.min()).days) + 1,192 "p_nochange_first": float(fed_main["No change"].iloc[0]), "p_nochange_last": float(fed_main["No change"].iloc[-1]),193 "p_cut25_first": float(fed_main["25 bps decrease"].iloc[0]), "p_cut25_last": float(fed_main["25 bps decrease"].iloc[-1]),194 "dff_start": float(dff.loc[fed_main.index.min()]), "dff_end": float(dff.loc[fed_main.index.max()]),195 "dff_before_listing": float(dff.loc[listing_day - pd.Timedelta(days=1)]),196 "min_price_by_state": fed_main.min().round(4).to_dict(),197 "volume_usd": {m["name"]: m["volume"] for m in fed_mk},198 "liquidity_usd": {m["name"]: m["liquidity"] for m in fed_mk},199 "total_volume_usd": float(sum(m["volume"] for m in fed_mk)),200}201202# ----------------------------------------------------------------------------203# 2. Table 1 — markets analysed (descriptive)204# ----------------------------------------------------------------------------205rows = []206for m in fed_mk:207 s = fed[m["name"]].dropna().iloc[1:] # drop seeded listing quote208 rows.append(f" {esc(m['name'])} & {f(m['volume'],0)} & {f(m['liquidity'],0)} & {len(s)} & "209 f"{f(s.iloc[0])} & {f(s.mean())} & {f(s.min())} & {f(s.max())} & {f(s.iloc[-1])} \\\\")210fed_rows = "\n".join(rows)211rows = []212for m in el_order:213 s = election[m["name"]].dropna().iloc[1:] if m["name"] in election else pd.Series(dtype=float)214 rows.append(f" {esc(m['name'])} & {f(m['volume']/1e6,1)} & -- & {len(s)} & "215 f"{f(s.iloc[0]) if len(s) else '--'} & {f(s.mean()) if len(s) else '--'} & {f(s.min()) if len(s) else '--'} & "216 f"{f(s.max()) if len(s) else '--'} & {f(s.iloc[-1]) if len(s) else '--'} \\\\")217el_rows = "\n".join(rows)218tex = rf"""\begin{{table}}[htbp]219\centering\footnotesize220\begin{{threeparttable}}221\caption{{Markets analysed: Polymarket metadata and daily closing prices}}222\label{{tab:descriptive_stats}}223\setlength{{\tabcolsep}}{{4pt}}224\begin{{tabular}}{{l r r r r r r r r}}225\toprule226Contract & Volume & Liquidity & $T$ & First & Mean & Min & Max & Last \\227\midrule228\multicolumn{{9}}{{l}}{{\textit{{Panel A: Fed decision, June 2026 FOMC meeting (event 101772; volume and liquidity in USD)}}}} \\229{fed_rows}230\midrule231\multicolumn{{9}}{{l}}{{\textit{{Panel B: Presidential Election Winner 2024 (event 903193; volume, USD million)}}}} \\232{el_rows}233\bottomrule234\end{{tabular}}235\begin{{tablenotes}}236\small237\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.238\end{{tablenotes}}239\end{{threeparttable}}240\end{{table}}"""241save_table(tex, "table1_descriptive_stats")242243# ----------------------------------------------------------------------------244# 3. Table 2 — summary statistics of state prices and daily changes245# ----------------------------------------------------------------------------246rows = []247for c in STATES:248 s = fed_main[c]249 d = s.diff().dropna()250 jb = stats.jarque_bera(d)251 adf = __import__("statsmodels.tsa.stattools", fromlist=["adfuller"]).adfuller(d, autolag="AIC")252 rows.append(f" {esc(c)} & {f(s.mean())} & {f(s.std())} & {f(s.min())} & {f(s.max())} & "253 f"{f(d.mean()*100,3)} & {f(d.std()*100,2)} & {f(d.skew(),2)} & {f(d.kurtosis(),1)} & {f(adf[0],2)} \\\\")254sig = sigma255d = sig.diff().dropna()256rows.append(f" State-price sum $\\Sigma_t$ & {f(sig.mean(),4)} & {f(sig.std(),4)} & {f(sig.min(),4)} & {f(sig.max(),4)} & "257 f"{f(d.mean()*100,3)} & {f(d.std()*100,2)} & {f(d.skew(),2)} & {f(d.kurtosis(),1)} & "258 f"{f(__import__('statsmodels.tsa.stattools', fromlist=['adfuller']).adfuller(d, autolag='AIC')[0],2)} \\\\")259tex = rf"""\begin{{table}}[htbp]260\centering\small261\begin{{threeparttable}}262\caption{{Summary statistics: Fed June 2026 state prices (11 December 2025 -- 22 May 2026)}}263\label{{tab:summary_stats}}264\setlength{{\tabcolsep}}{{4pt}}265\begin{{tabular}}{{l rrrr rrrr r}}266\toprule267 & \multicolumn{{4}}{{c}}{{Level $p_{{j,t}}$}} & \multicolumn{{5}}{{c}}{{Daily change $\Delta p_{{j,t}}$ (percentage points)}} \\268\cmidrule(lr){{2-5}}\cmidrule(lr){{6-10}}269Contract & Mean & SD & Min & Max & Mean & SD & Skew & Ex.\ kurt & ADF \\270\midrule271{chr(10).join(rows)}272\bottomrule273\end{{tabular}}274\begin{{tablenotes}}275\small276\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.277\end{{tablenotes}}278\end{{threeparttable}}279\end{{table}}"""280save_table(tex, "table2_summary_stats")281282# ----------------------------------------------------------------------------283# 4. Completeness / no-arbitrage on the complete set (Table 3)284# ----------------------------------------------------------------------------285from statsmodels.tsa.stattools import adfuller # noqa: E402286287288def completeness_block(sig: pd.Series, resolution: pd.Timestamp, n_states: int, tick: float, label: str) -> dict:289 dev = sig - 1.0290 tau = (resolution - sig.index).days.values / 365.0291 r = dff.reindex(sig.index).ffill().bfill().values / 100.0292 lower = np.exp(-r * tau) # discounted value of the complete set (buy-and-hold bound)293 dev_disc = sig.values - lower294 mean, se, t, p = nw_mean_test(dev)295 adf = adfuller(dev, autolag="AIC", regression="c")296 phi = np.polyfit(dev.values[:-1], dev.values[1:], 1)[0]297 hl = -np.log(2) / np.log(abs(phi)) if 0 < abs(phi) < 1 else np.nan298 band = 0.005 # half a cent on the complete set299 out = {300 "label": label, "T": int(len(sig)), "mean": float(sig.mean()), "sd": float(sig.std()),301 "median": float(sig.median()), "min": float(sig.min()), "max": float(sig.max()),302 "mean_dev": float(dev.mean()), "mean_abs_dev": float(dev.abs().mean()), "max_abs_dev": float(dev.abs().max()),303 "nw_se": se, "nw_t": t, "nw_p": p,304 "share_above_1": float((sig > 1).mean()), "share_below_1": float((sig < 1).mean()),305 "share_within_tickband": float((dev.abs() <= band).mean()),306 "share_within_1pp": float((dev.abs() <= 0.01).mean()), "share_within_2pp": float((dev.abs() <= 0.02).mean()),307 "share_within_5pp": float((dev.abs() <= 0.05).mean()),308 "mean_lower_bound": float(lower.mean()), "min_lower_bound": float(lower.min()),309 "mean_dev_from_lower": float(dev_disc.mean()), "share_inside_band": float(((sig.values >= lower - band) & (sig.values <= 1 + band)).mean()),310 "share_above_upper_band": float((sig.values > 1 + band).mean()),311 "adf_stat": float(adf[0]), "adf_p": float(adf[1]), "ar1": float(phi), "half_life": float(hl),312 "tick_band": band,313 "first_half_mean": float(sig.iloc[: len(sig) // 2].mean()), "second_half_mean": float(sig.iloc[len(sig) // 2:].mean()),314 "last30_mean": float(sig.iloc[-30:].mean()), "last30_sd": float(sig.iloc[-30:].std()),315 "first30_mean": float(sig.iloc[:30].mean()), "first30_sd": float(sig.iloc[:30].std()),316 }317 # overround vs. time-to-resolution regression318 import statsmodels.api as sm319 X = sm.add_constant(tau)320 m = sm.OLS(dev.values, X).fit(cov_type="HAC", cov_kwds={"maxlags": 5})321 out["slope_tau"] = float(m.params[1]); out["slope_tau_t"] = float(m.tvalues[1]); out["const_tau"] = float(m.params[0])322 return out, lower323324325el_complete_all = election.dropna()326el_complete = el_complete_all.iloc[1:] # first complete day carries seeded listing quotes327sigma_el = el_complete.sum(axis=1)328sigma_el_all = el_complete_all.sum(axis=1)329C_fed, lower_fed = completeness_block(sigma, RESOLUTION_FED, 5, TICK, "Fed June 2026")330C_fed_all, _ = completeness_block(sigma_all, RESOLUTION_FED, 5, TICK, "Fed incl. listing day")331C_el, lower_el = completeness_block(sigma_el, RESOLUTION_ELECTION, len(el_complete.columns), 0.001, "Election 2024")332C_el_all, _ = completeness_block(sigma_el_all, RESOLUTION_ELECTION, len(el_complete.columns), 0.001, "Election incl. listing day")333STATS["completeness"] = {"fed": C_fed, "fed_incl_listing": C_fed_all, "election": C_el, "election_incl_listing": C_el_all,334 "election_n_contracts": int(len(el_complete.columns)), "election_T_union": int(len(election)),335 "election_first": str(election.index.min().date()), "election_last": str(election.index.max().date())}336337338def row(name, a, b, d=4, pct=False):339 fa = f(a * 100, 1) + r"\%" if pct else f(a, d)340 fb = f(b * 100, 1) + r"\%" if pct else f(b, d)341 return f" {name} & {fa} & {fb} \\\\"342343344tex = rf"""\begin{{table}}[htbp]345\centering\small346\begin{{threeparttable}}347\caption{{Price of the complete set of state claims: law-of-one-price tests}}348\label{{tab:completeness}}349\begin{{tabular}}{{l r r}}350\toprule351 & Fed June 2026 & Election 2024 \\352 & ($N=5$ states) & ($N={len(el_complete.columns)}$ candidates) \\353\midrule354\multicolumn{{3}}{{l}}{{\textit{{Panel A: distribution of $\Sigma_t=\sum_j p_{{j,t}}$}}}} \\355{row('Observations $T$', C_fed['T'], C_el['T'], 0)}356{row('Mean', C_fed['mean'], C_el['mean'])}357{row('Median', C_fed['median'], C_el['median'])}358{row('Standard deviation', C_fed['sd'], C_el['sd'])}359{row('Minimum', C_fed['min'], C_el['min'])}360{row('Maximum', C_fed['max'], C_el['max'])}361\midrule362\multicolumn{{3}}{{l}}{{\textit{{Panel B: deviation $\Sigma_t-1$}}}} \\363{row('Mean deviation (overround)', C_fed['mean_dev'], C_el['mean_dev'])}364{row('HAC standard error (5 lags)', C_fed['nw_se'], C_el['nw_se'])}365{row('$t$-statistic, $H_0$: $\\mathbb{{E}}[\\Sigma_t-1]=0$', C_fed['nw_t'], C_el['nw_t'], 2)}366{row('Mean absolute deviation', C_fed['mean_abs_dev'], C_el['mean_abs_dev'])}367{row('Share of days with $\\Sigma_t>1$', C_fed['share_above_1'], C_el['share_above_1'], pct=True)}368{row('Share within $\\pm 0.5$ pp', C_fed['share_within_tickband'], C_el['share_within_tickband'], pct=True)}369{row('Share within $\\pm 1$ pp', C_fed['share_within_1pp'], C_el['share_within_1pp'], pct=True)}370{row('Share within $\\pm 2$ pp', C_fed['share_within_2pp'], C_el['share_within_2pp'], pct=True)}371{row('Share within $\\pm 5$ pp', C_fed['share_within_5pp'], C_el['share_within_5pp'], pct=True)}372\midrule373\multicolumn{{3}}{{l}}{{\textit{{Panel C: discount-adjusted bounds, $e^{{-r_t\tau_t}}\le\Sigma_t\le 1$}}}} \\374{row('Mean lower bound $e^{{-r_t\\tau_t}}$', C_fed['mean_lower_bound'], C_el['mean_lower_bound'])}375{row('Mean $\\Sigma_t-e^{{-r_t\\tau_t}}$', C_fed['mean_dev_from_lower'], C_el['mean_dev_from_lower'])}376{row('Share inside bounds ($\\pm 0.5$ pp)', C_fed['share_inside_band'], C_el['share_inside_band'], pct=True)}377{row('Share above upper bound', C_fed['share_above_upper_band'], C_el['share_above_upper_band'], pct=True)}378\midrule379\multicolumn{{3}}{{l}}{{\textit{{Panel D: dynamics of $\Sigma_t-1$}}}} \\380{row('AR(1) coefficient', C_fed['ar1'], C_el['ar1'], 3)}381{row('Half-life (days)', C_fed['half_life'], C_el['half_life'], 1)}382{row('ADF statistic (constant, AIC lag)', C_fed['adf_stat'], C_el['adf_stat'], 2)}383{row('ADF $p$-value', C_fed['adf_p'], C_el['adf_p'], 3)}384{row('Slope on time-to-resolution $\\tau_t$ (years)', C_fed['slope_tau'], C_el['slope_tau'], 3)}385{row('\\quad HAC $t$-statistic', C_fed['slope_tau_t'], C_el['slope_tau_t'], 2)}386\bottomrule387\end{{tabular}}388\begin{{tablenotes}}389\small390\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.391\end{{tablenotes}}392\end{{threeparttable}}393\end{{table}}"""394save_table(tex, "table3_completeness")395396# ----------------------------------------------------------------------------397# 5. Order-book snapshot (Table 9)398# ----------------------------------------------------------------------------399ob_rows, ob = [], {}400for m in fed_mk:401 path = DATA_RAW_POLY / f"order_book_{m['yes_token'][:16]}.json"402 if not path.exists():403 continue404 d = json.load(open(path))405 bids = sorted([(float(b["price"]), float(b["size"])) for b in d["bids"]], reverse=True)406 asks = sorted([(float(a["price"]), float(a["size"])) for a in d["asks"]])407 bb, ba = bids[0][0], asks[0][0]408 ob[m["name"]] = {"bid": bb, "ask": ba, "mid": (bb + ba) / 2, "spread": ba - bb,409 "rel_spread_pct": 100 * (ba - bb) / ((bb + ba) / 2),410 "depth_best_bid": bids[0][1], "depth_best_ask": asks[0][1],411 "depth_bid_usd": sum(p * s for p, s in bids), "depth_ask_usd": sum(p * s for p, s in asks),412 "shares_bid": sum(s for _, s in bids), "shares_ask": sum(s for _, s in asks),413 "depth_within_1c_bid_usd": sum(p * s for p, s in bids if p >= bb - 0.01),414 "depth_within_1c_ask_usd": sum(p * s for p, s in asks if p <= ba + 0.01),415 "timestamp": int(d["timestamp"])}416STATS["orderbook"] = ob417STATS["orderbook"]["sum_bid"] = float(sum(v["bid"] for k, v in ob.items() if isinstance(v, dict)))418STATS["orderbook"]["sum_ask"] = float(sum(v["ask"] for k, v in ob.items() if isinstance(v, dict)))419STATS["orderbook"]["sum_mid"] = float(sum(v["mid"] for k, v in ob.items() if isinstance(v, dict)))420STATS["orderbook"]["snapshot_utc"] = str(pd.to_datetime(list(ob.values())[0]["timestamp"], unit="ms"))421for c in STATES:422 v = ob[c]423 ob_rows.append(f" {esc(c)} & {v['bid']:.3f} & {v['ask']:.3f} & {v['spread']:.3f} & {f(v['rel_spread_pct'],1)} & "424 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)} \\\\")425tex = rf"""\begin{{table}}[htbp]426\centering\small427\begin{{threeparttable}}428\caption{{Order-book snapshot, Fed June 2026 contracts (22 May 2026, {STATS['orderbook']['snapshot_utc'][11:16]} UTC)}}429\label{{tab:orderbook}}430\setlength{{\tabcolsep}}{{4pt}}431\begin{{tabular}}{{l rrrr rr rr}}432\toprule433 & Best & Best & Quoted & Relative & \multicolumn{{2}}{{c}}{{Depth at best (shares)}} & \multicolumn{{2}}{{c}}{{Depth within 1\textcent\ (USD 000)}} \\434\cmidrule(lr){{6-7}}\cmidrule(lr){{8-9}}435Contract & bid & ask & spread & spread (\%) & Bid & Ask & Bid & Ask \\436\midrule437{chr(10).join(ob_rows)}438\midrule439Sum of quotes & {STATS['orderbook']['sum_bid']:.3f} & {STATS['orderbook']['sum_ask']:.3f} & & & & & & \\440\bottomrule441\end{{tabular}}442\begin{{tablenotes}}443\small444\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.445\end{{tablenotes}}446\end{{threeparttable}}447\end{{table}}"""448save_table(tex, "table9_orderbook")449450# ----------------------------------------------------------------------------451# 6. Risk-neutral density, expected policy change, announcement days (Table 4)452# ----------------------------------------------------------------------------453q = fed_main.div(sigma, axis=0)454eq = (q.values @ STATE_BPS)455eq = pd.Series(eq, index=fed_main.index, name="EQ_bps")456var_q = pd.Series(((STATE_BPS[None, :] - eq.values[:, None]) ** 2 * q.values).sum(axis=1), index=fed_main.index)457d_eq = eq.diff().dropna()458d_pnc = fed_main["No change"].diff().dropna() * 100 # pp459460def _map(d0: pd.Timestamp) -> pd.Timestamp | None:461 """First available Polymarket close on or after d0 (gaps of at most 2 days)."""462 pos = d_eq.index.searchsorted(d0)463 if pos < len(d_eq.index) and (d_eq.index[pos] - d0).days <= 2:464 return d_eq.index[pos]465 return None466467468EVENT_MAP = {}469def in_sample(ds):470 out = []471 for x in ds:472 x = pd.Timestamp(x)473 if fed_main.index.min() < x <= fed_main.index.max():474 m = _map(x)475 if m is not None:476 EVENT_MAP[str(m.date())] = str(x.date()); out.append(m)477 return out478fomc = in_sample(FOMC_DATES)479cpi = in_sample(CPI_RELEASE_DATES)480nfp = in_sample(NFP_RELEASE_DATES)481ann = {"FOMC": fomc, "CPI": cpi, "Employment report": nfp}482all_ann = sorted(set(fomc + cpi + nfp))483STATS["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],484 "eq_first": float(eq.iloc[0]), "eq_last": float(eq.iloc[-1]), "eq_min": float(eq.min()), "eq_max": float(eq.max()),485 "sd_q_first": float(np.sqrt(var_q.iloc[0])), "sd_q_last": float(np.sqrt(var_q.iloc[-1]))}486487488def cum(series: pd.Series, d0: pd.Timestamp, a: int, b: int) -> float:489 idx = series.index490 pos = idx.searchsorted(d0)491 if pos >= len(idx) or idx[pos] != d0:492 return np.nan493 lo, hi = pos + a, pos + b494 if lo < 0 or hi >= len(idx):495 return np.nan496 return float(series.iloc[lo:hi + 1].sum())497498499ev_rows, ev_records = [], []500for etype, ds in ann.items():501 for d0 in ds:502 rec = {"type": etype, "date": str(d0.date()),503 "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),504 "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),505 "eq_level": float(eq.loc[d0]) if d0 in eq.index else np.nan}506 ev_records.append(rec)507 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)} & "508 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)} \\\\")509STATS["events"]["records"] = ev_records510511# announcement-day vs other-day comparison512is_ann = d_eq.index.isin(all_ann)513is_wkend = d_eq.index.dayofweek >= 5514515516def compare(x: pd.Series, mask: np.ndarray, label: str) -> dict:517 a, b = x[mask].abs(), x[~mask].abs()518 u = stats.mannwhitneyu(a, b, alternative="greater")519 lev = stats.levene(x[mask], x[~mask], center="median")520 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()),521 "median_abs_a": float(a.median()), "median_abs_b": float(b.median()),522 "ratio_mean": float(a.mean() / b.mean()), "sd_a": float(x[mask].std()), "sd_b": float(x[~mask].std()),523 "var_ratio": float(x[mask].var() / x[~mask].var()), "mw_p": float(u.pvalue), "bf_stat": float(lev.statistic), "bf_p": float(lev.pvalue),524 "mean_a": float(x[mask].mean()), "mean_b": float(x[~mask].mean())}525526527cmp_rows = []528COMP = {}529for lab, mask in [("All scheduled releases", is_ann), ("FOMC", d_eq.index.isin(fomc)), ("CPI", d_eq.index.isin(cpi)),530 ("Employment report", d_eq.index.isin(nfp))]:531 c1 = compare(d_eq, mask, lab); c2 = compare(d_pnc, mask, lab)532 COMP[lab] = {"eq": c1, "pnc": c2}533 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)} & "534 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)} \\\\")535# weekend vs weekday (non-announcement days only)536wk_mask = is_wkend & ~is_ann537c1 = compare(d_eq[~is_ann], wk_mask[~is_ann], "Weekend"); c2 = compare(d_pnc[~is_ann], wk_mask[~is_ann], "Weekend")538COMP["Weekend"] = {"eq": c1, "pnc": c2}539# TNX benchmark: |Δ10y yield| on same announcement days540tnx_d = (tnx.diff() * 100).dropna()541tnx_d = tnx_d[(tnx_d.index > fed_main.index.min()) & (tnx_d.index <= fed_main.index.max())]542ct = compare(tnx_d, tnx_d.index.isin(all_ann), "TNX all releases")543COMP["TNX"] = ct544STATS["announcement"] = COMP545546tex = rf"""\begin{{table}}[htbp]547\centering\footnotesize548\begin{{threeparttable}}549\caption{{Scheduled macroeconomic releases and the risk-neutral expected policy change}}550\label{{tab:event_study}}551\setlength{{\tabcolsep}}{{4pt}}552\begin{{tabular}}{{l l r rrr rrr}}553\toprule554\multicolumn{{9}}{{l}}{{\textit{{Panel A: individual releases in the sample}}}} \\555\addlinespace[2pt]556 & & $\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)}} \\557\cmidrule(lr){{4-6}}\cmidrule(lr){{7-9}}558Release & Date & level (bps) & $[0]$ & $[-1,+1]$ & $[-3,+3]$ & $[0]$ & $[-1,+1]$ & $[-3,+3]$ \\559\midrule560{chr(10).join(ev_rows)}561\bottomrule562\end{{tabular}}\par\bigskip563\begin{{tabular}}{{l r rrrrr rrrrr}}564\toprule565\multicolumn{{12}}{{l}}{{\textit{{Panel B: absolute daily changes on release days versus all other days}}}} \\566\addlinespace[2pt]567 & & \multicolumn{{5}}{{c}}{{$|\Delta \mathbb{{E}}^{{\mathbb{{Q}}}}[\Delta r]|$ (bps)}} & \multicolumn{{5}}{{c}}{{$|\Delta p_{{\text{{No change}}}}|$ (pp)}} \\568\cmidrule(lr){{3-7}}\cmidrule(lr){{8-12}}569Days & $K$ & Ann. & Other & VR & $p_{{BF}}$ & $p_{{MW}}$ & Ann. & Other & VR & $p_{{BF}}$ & $p_{{MW}}$ \\570\midrule571"""572for lab in ["All scheduled releases", "FOMC", "CPI", "Employment report"]:573 c1, c2 = COMP[lab]["eq"], COMP[lab]["pnc"]574 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)} & "575 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")576c1, c2 = COMP["Weekend"]["eq"], COMP["Weekend"]["pnc"]577tex += (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)} & "578 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")579tex += rf"""\bottomrule580\end{{tabular}}581\begin{{tablenotes}}582\small583\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)}$).584\end{{tablenotes}}585\end{{threeparttable}}586\end{{table}}"""587save_table(tex, "table4_event_study")588589# ----------------------------------------------------------------------------590# 7. Granger causality / VAR on trading-day grid (Tables 5, 6)591# ----------------------------------------------------------------------------592from statsmodels.tsa.api import VAR # noqa: E402593import statsmodels.api as sm # noqa: E402594595grid = spy.index[(spy.index >= fed_main.index.min()) & (spy.index <= fed_main.index.max())]596eq_grid = eq.reindex(grid).ffill()597pnc_grid = fed_main["No change"].reindex(grid).ffill()598panel = pd.DataFrame({599 "Fed": eq_grid.diff(), # Δ E_Q[Δr], bps600 "TNX": (tnx.reindex(grid).diff() * 100), # Δ 10-year yield, bps601 "SPY": np.log(spy.reindex(grid)).diff() * 100,602 "VIX": np.log(vix.reindex(grid)).diff() * 100,603 "BTC": np.log(btc.reindex(grid)).diff() * 100,604}).dropna()605STATS["granger"] = {"T": int(len(panel)), "first": str(panel.index.min().date()), "last": str(panel.index.max().date())}606sel = VAR(panel).select_order(maxlags=5)607STATS["granger"]["aic_lag"] = int(sel.aic); STATS["granger"]["bic_lag"] = int(sel.bic); STATS["granger"]["hqic_lag"] = int(sel.hqic)608p_main = max(int(sel.aic), 1)609610611def gc_pair(y: str, x: str, p: int):612 """Bivariate VAR(p) Granger F-test of x → y with HAC-robust Wald as well."""613 df = panel[[y, x]]614 r = VAR(df).fit(p)615 t = r.test_causality(y, [x], kind="f")616 # HAC-robust Wald on single equation617 Y = df[y].iloc[p:]618 X = pd.concat([df.shift(k).add_suffix(f"_L{k}") for k in range(1, p + 1)], axis=1).iloc[p:]619 X = sm.add_constant(X)620 m = sm.OLS(Y, X).fit(cov_type="HAC", cov_kwds={"maxlags": 5})621 R = np.zeros((p, X.shape[1]))622 for k in range(1, p + 1):623 R[k - 1, list(X.columns).index(f"{x}_L{k}")] = 1624 w = m.wald_test(R, use_f=True, scalar=True)625 return float(t.test_statistic), float(t.pvalue), float(w.statistic), float(w.pvalue)626627628names = {"Fed": r"$\Delta\mathbb{E}^{\mathbb{Q}}[\Delta r]$", "TNX": "10-year yield", "SPY": "SPY", "VIX": "VIX", "BTC": "BTC"}629g_rows, G = [], {}630pairs = [("Fed", v) for v in ["TNX", "SPY", "VIX", "BTC"]] + [(v, "Fed") for v in ["TNX", "SPY", "VIX", "BTC"]]631for x, y in pairs:632 res = {}633 for p in [p_main, 1, 2, 5]:634 res[p] = gc_pair(y, x, p)635 G[f"{x}->{y}"] = {str(k): v for k, v in res.items()}636 Fm, pm, Wm, pwm = res[p_main]637 g_rows.append(f" {names[x]} $\\to$ {names[y]} & {f(Fm,2)}{stars(pm)} & {f(pm,3)} & {f(pwm,3)} & "638 f"{f(res[1][1],3)} & {f(res[2][1],3)} & {f(res[5][1],3)} \\\\")639STATS["granger"]["pairs"] = G640# old-style (min-p over lags 1..5, misaligned returns) for the changelog641from statsmodels.tsa.stattools import grangercausalitytests # noqa: E402642old = {}643for x, y in [("Fed", "VIX"), ("VIX", "Fed"), ("Fed", "SPY")]:644 r = grangercausalitytests(panel[[y, x]].values, maxlag=5, verbose=False)645 best = min(r, key=lambda k: r[k][0]["ssr_ftest"][1])646 old[f"{x}->{y}"] = {"minp_lag": int(best), "F": float(r[best][0]["ssr_ftest"][0]), "p": float(r[best][0]["ssr_ftest"][1])}647STATS["granger"]["min_p_over_lags_aligned"] = old648649tex = rf"""\begin{{table}}[htbp]650\centering\small651\begin{{threeparttable}}652\caption{{Granger causality between the Fed state-price panel and U.S. asset returns}}653\label{{tab:granger}}654\setlength{{\tabcolsep}}{{4pt}}655\begin{{tabular}}{{l rrr rrr}}656\toprule657 & \multicolumn{{3}}{{c}}{{VAR($p^*$), $p^*={p_main}$ (AIC)}} & \multicolumn{{3}}{{c}}{{$p$-value at fixed lag}} \\658\cmidrule(lr){{2-4}}\cmidrule(lr){{5-7}}659Direction & $F$ & $p$ & $p$ (HAC) & $p=1$ & $p=2$ & $p=5$ \\660\midrule661{chr(10).join(g_rows)}662\bottomrule663\end{{tabular}}664\begin{{tablenotes}}665\small666\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\%.667\end{{tablenotes}}668\end{{threeparttable}}669\end{{table}}"""670save_table(tex, "table5_granger")671672# VAR(1) coefficient table with HAC t-stats (equation by equation)673var_rows = []674Xv = sm.add_constant(panel.shift(1).iloc[1:])675coef, tval = {}, {}676for y in panel.columns:677 m = sm.OLS(panel[y].iloc[1:], Xv).fit(cov_type="HAC", cov_kwds={"maxlags": 5})678 coef[y] = m.params; tval[y] = m.tvalues679 STATS["granger"].setdefault("var1_r2", {})[y] = float(m.rsquared)680for reg in ["const"] + list(panel.columns):681 lab = "Constant" if reg == "const" else f"{names[reg]}$_{{t-1}}$"682 cells, tcells = [], []683 for y in panel.columns:684 c, t = coef[y][reg], tval[y][reg]685 p = 2 * (1 - stats.norm.cdf(abs(t)))686 cells.append(f"{f(c,3)}{stars(p)}")687 tcells.append(f"({f(t,2)})")688 var_rows.append(f" {lab} & " + " & ".join(cells) + " \\\\")689 var_rows.append(" & " + " & ".join(tcells) + " \\\\")690r2 = " & ".join(f(STATS["granger"]["var1_r2"][y], 3) for y in panel.columns)691tex = rf"""\begin{{table}}[htbp]692\centering\small693\begin{{threeparttable}}694\caption{{VAR(1) estimates on the trading-day grid}}695\label{{tab:var}}696\setlength{{\tabcolsep}}{{4pt}}697\begin{{tabular}}{{l rrrrr}}698\toprule699 & \multicolumn{{5}}{{c}}{{Dependent variable (day $t$)}} \\700\cmidrule(lr){{2-6}}701Regressor & {names['Fed']} & 10-year yield & SPY & VIX & BTC \\702\midrule703{chr(10).join(var_rows)}704\midrule705$R^2$ & {r2} \\706Observations & \multicolumn{{5}}{{c}}{{{len(panel)-1}}} \\707\bottomrule708\end{{tabular}}709\begin{{tablenotes}}710\small711\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\%.712\end{{tablenotes}}713\end{{threeparttable}}714\end{{table}}"""715save_table(tex, "table6_var")716717# contemporaneous correlations (for text)718STATS["granger"]["contemp_corr"] = panel.corr().round(3).to_dict()719720# ----------------------------------------------------------------------------721# 8. Volatility: GARCH / GJR (QML) and HAR (appendix tables 7, 8)722# ----------------------------------------------------------------------------723from arch import arch_model # noqa: E402724725r_nc = np.log(fed_main["No change"]).diff().dropna() * 100 # percent log returns, calendar days726STATS["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())}727# ARCH-LM728from statsmodels.stats.diagnostic import het_arch # noqa: E402729lm = het_arch(r_nc - r_nc.mean(), nlags=5)730STATS["garch"]["arch_lm_stat"] = float(lm[0]); STATS["garch"]["arch_lm_p"] = float(lm[1])731lb = __import__("statsmodels.stats.diagnostic", fromlist=["acorr_ljungbox"]).acorr_ljungbox(r_nc ** 2, lags=[10], return_df=True)732STATS["garch"]["lb_sq_10_stat"] = float(lb["lb_stat"].iloc[0]); STATS["garch"]["lb_sq_10_p"] = float(lb["lb_pvalue"].iloc[0])733lb1 = __import__("statsmodels.stats.diagnostic", fromlist=["acorr_ljungbox"]).acorr_ljungbox(r_nc, lags=[10], return_df=True)734STATS["garch"]["lb_ret_10_stat"] = float(lb1["lb_stat"].iloc[0]); STATS["garch"]["lb_ret_10_p"] = float(lb1["lb_pvalue"].iloc[0])735736garch_fits = {}737for label, o, dist in [("GARCH(1,1), normal", 0, "normal"), ("GJR-GARCH(1,1), normal", 1, "normal"),738 ("GARCH(1,1), Student-$t$", 0, "t"), ("GJR-GARCH(1,1), Student-$t$", 1, "t")]:739 am = arch_model(r_nc, vol="GARCH", p=1, o=o, q=1, dist=dist, mean="Constant")740 res = am.fit(disp="off", cov_type="robust")741 garch_fits[label] = res742 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()},743 "pvalues": {k: float(v) for k, v in res.pvalues.items()},744 "loglik": float(res.loglikelihood), "aic": float(res.aic), "bic": float(res.bic),745 "persistence": float(res.params.get("alpha[1]", 0) + res.params.get("beta[1]", 0) + 0.5 * res.params.get("gamma[1]", 0))}746 # Ljung-Box on squared standardised residuals747 z = res.std_resid.dropna()748 lbz = __import__("statsmodels.stats.diagnostic", fromlist=["acorr_ljungbox"]).acorr_ljungbox(z ** 2, lags=[10], return_df=True)749 STATS["garch"][label]["lb_z2_10_p"] = float(lbz["lb_pvalue"].iloc[0])750 STATS["garch"][label]["kurt_z"] = float(stats.kurtosis(z, fisher=False))751# LR test GARCH vs GJR (normal)752lr = 2 * (garch_fits["GJR-GARCH(1,1), normal"].loglikelihood - garch_fits["GARCH(1,1), normal"].loglikelihood)753STATS["garch"]["lr_gjr_vs_garch"] = float(lr); STATS["garch"]["lr_p"] = float(1 - stats.chi2.cdf(lr, 1))754755param_order = ["mu", "omega", "alpha[1]", "gamma[1]", "beta[1]", "nu"]756pnames = {"mu": r"$\mu$", "omega": r"$\omega$", "alpha[1]": r"$\alpha_1$", "gamma[1]": r"$\gamma_1$", "beta[1]": r"$\beta_1$", "nu": r"$\nu$"}757g_rows = []758labels = list(garch_fits.keys())759for pn in param_order:760 cells, secells = [], []761 for lab in labels:762 res = garch_fits[lab]763 if pn in res.params:764 cells.append(f"{f(res.params[pn],3)}{stars(res.pvalues[pn])}"); secells.append(f"({f(res.std_err[pn],3)})")765 else:766 cells.append("--"); secells.append("")767 g_rows.append(f" {pnames[pn]} & " + " & ".join(cells) + " \\\\")768 g_rows.append(" & " + " & ".join(secells) + " \\\\")769ll = " & ".join(f(garch_fits[l].loglikelihood, 1) for l in labels)770aic = " & ".join(f(garch_fits[l].aic, 1) for l in labels)771bic = " & ".join(f(garch_fits[l].bic, 1) for l in labels)772pers = " & ".join(f(STATS["garch"][l]["persistence"], 3) for l in labels)773lbz = " & ".join(f(STATS["garch"][l]["lb_z2_10_p"], 3) for l in labels)774tex = rf"""\begin{{table}}[htbp]775\centering\small776\begin{{threeparttable}}777\caption{{GARCH-family estimates for the ``No change'' contract (daily log returns, percent)}}778\label{{tab:garch}}779\setlength{{\tabcolsep}}{{4pt}}780\begin{{tabular}}{{l rrrr}}781\toprule782 & GARCH & GJR-GARCH & GARCH & GJR-GARCH \\783Parameter & normal & normal & Student-$t$ & Student-$t$ \\784\midrule785{chr(10).join(g_rows)}786\midrule787Log-likelihood & {ll} \\788AIC & {aic} \\789BIC & {bic} \\790Persistence $\alpha_1+\beta_1+\gamma_1/2$ & {pers} \\791Ljung--Box $p$, $z_t^2$ (10 lags) & {lbz} \\792Observations & \multicolumn{{4}}{{c}}{{{len(r_nc)}}} \\793\bottomrule794\end{{tabular}}795\begin{{tablenotes}}796\small797\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\%.798\end{{tablenotes}}799\end{{threeparttable}}800\end{{table}}"""801save_table(tex, "table7_garch")802803# HAR with standard component definitions804rv_d = r_nc.abs()805har = pd.DataFrame({"rv_d": rv_d, "rv_w": rv_d.rolling(5).mean(), "rv_m": rv_d.rolling(22).mean()})806y = rv_d.shift(-1)807dat = pd.concat([y.rename("y"), har], axis=1).dropna()808hm = sm.OLS(dat["y"], sm.add_constant(dat[["rv_d", "rv_w", "rv_m"]])).fit(cov_type="HAC", cov_kwds={"maxlags": 5})809# squared-return version810rv2 = r_nc ** 2811har2 = pd.DataFrame({"rv_d": rv2, "rv_w": rv2.rolling(5).mean(), "rv_m": rv2.rolling(22).mean()})812dat2 = pd.concat([rv2.shift(-1).rename("y"), har2], axis=1).dropna()813hm2 = sm.OLS(dat2["y"], sm.add_constant(dat2[["rv_d", "rv_w", "rv_m"]])).fit(cov_type="HAC", cov_kwds={"maxlags": 5})814STATS["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),815 "F_p": float(hm.f_pvalue)},816 "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),817 "F_p": float(hm2.f_pvalue)}}818hn = {"const": "Constant", "rv_d": r"$RV^{(d)}_t$", "rv_w": r"$RV^{(w)}_t$", "rv_m": r"$RV^{(m)}_t$"}819h_rows = []820for k in ["const", "rv_d", "rv_w", "rv_m"]:821 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)}) \\\\")822tex = rf"""\begin{{table}}[htbp]823\centering\small824\begin{{threeparttable}}825\caption{{HAR regressions for the ``No change'' contract}}826\label{{tab:har_rv}}827\begin{{tabular}}{{l rr rr}}828\toprule829 & \multicolumn{{2}}{{c}}{{$RV_t=|r_t|$}} & \multicolumn{{2}}{{c}}{{$RV_t=r_t^2$}} \\830\cmidrule(lr){{2-3}}\cmidrule(lr){{4-5}}831Regressor & Coef. & (SE) & Coef. & (SE) \\832\midrule833{chr(10).join(h_rows)}834\midrule835$R^2$ & \multicolumn{{2}}{{c}}{{{f(hm.rsquared,3)}}} & \multicolumn{{2}}{{c}}{{{f(hm2.rsquared,3)}}} \\836Adjusted $R^2$ & \multicolumn{{2}}{{c}}{{{f(hm.rsquared_adj,3)}}} & \multicolumn{{2}}{{c}}{{{f(hm2.rsquared_adj,3)}}} \\837Joint $F$-test $p$-value & \multicolumn{{2}}{{c}}{{{f(hm.f_pvalue,3)}}} & \multicolumn{{2}}{{c}}{{{f(hm2.f_pvalue,3)}}} \\838Observations & \multicolumn{{2}}{{c}}{{{int(hm.nobs)}}} & \multicolumn{{2}}{{c}}{{{int(hm2.nobs)}}} \\839\bottomrule840\end{{tabular}}841\begin{{tablenotes}}842\small843\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\%.844\end{{tablenotes}}845\end{{threeparttable}}846\end{{table}}"""847save_table(tex, "table8_har_rv")848849# ----------------------------------------------------------------------------850# 9. Roll spread851# ----------------------------------------------------------------------------852def roll_spread(p: pd.Series, window: int = 20) -> pd.Series:853 dp = p.diff()854 cov = dp.rolling(window).apply(lambda x: np.cov(x[:-1], x[1:])[0, 1] if len(x) > 2 else np.nan, raw=True)855 return 2 * np.sqrt(np.maximum(-cov, 0))856857858roll = {c: roll_spread(fed_main[c]) for c in STATES}859rs = roll["No change"].dropna()860STATS["roll"] = {"nochange_mean": float(rs.mean()), "nochange_median": float(rs.median()), "nochange_first30_mean": float(rs.iloc[:30].mean()),861 "nochange_last30_mean": float(rs.iloc[-30:].mean()), "share_zero": float((rs == 0).mean()),862 "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())}863864# ----------------------------------------------------------------------------865# 10. Election panel statistics866# ----------------------------------------------------------------------------867tr, ha, bi = election["Donald Trump"], election.get("Kamala Harris"), election.get("Joe Biden")868STATS["election"] = {"trump_last": float(tr.dropna().iloc[-1]), "trump_max": float(tr.max()), "trump_min": float(tr.min()),869 "trump_mean": float(tr.mean()), "harris_last": float(ha.dropna().iloc[-1]), "harris_max": float(ha.max()),870 "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]),871 "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]),872 "trump_jul12": float(tr.loc[pd.Timestamp("2024-07-12")]), "trump_jul14": float(tr.loc[pd.Timestamp("2024-07-14")]),873 "trump_sep9": float(tr.loc[pd.Timestamp("2024-09-09")]), "trump_sep11": float(tr.loc[pd.Timestamp("2024-09-11")]),874 "trump_oct_peak": float(tr.loc["2024-10"].max()), "trump_oct_peak_date": str(tr.loc["2024-10"].idxmax().date()),875 "trump_nov5": float(tr.loc[pd.Timestamp("2024-11-05")]), "harris_nov5": float(ha.loc[pd.Timestamp("2024-11-05")]),876 "others_sum_nov5": float(el_complete.loc[pd.Timestamp("2024-11-05")].drop(["Donald Trump", "Kamala Harris"]).sum()),877 "n_contracts": int(election.shape[1]), "T": int(len(election)), "T_complete": int(len(el_complete))}878879# ----------------------------------------------------------------------------880# 11. Figures881# ----------------------------------------------------------------------------882# fig1 state prices883fig, ax = plt.subplots(figsize=(6.5, 3.6))884for i, c in enumerate(STATES):885 ax.plot(fed_main.index, fed_main[c], color=PAL[i], ls=LS[i], label=c)886for d0 in fomc:887 ax.axvline(d0, color="0.6", lw=0.6, ls=":")888ax.set_ylabel("State price (USDC per \\$1 payoff)"); ax.set_ylim(-0.02, 1.02)889ax.legend(ncol=3, loc="center left", frameon=False)890ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))891save_fig(fig, "fig1_state_prices")892893# fig2 Σ with bounds894fig, ax = plt.subplots(figsize=(6.5, 3.4))895ax.plot(sigma.index, sigma.values, color=PAL[0], lw=1.2, label=r"$\Sigma_t$ (sum of closing prices)")896ax.plot(sigma.index, lower_fed, color=PAL[1], lw=1.0, ls="--", label=r"Lower bound $e^{-r_t\tau_t}$ (DFF)")897ax.axhline(1.0, color="black", lw=0.8, label="Upper bound 1")898ax.fill_between(sigma.index, lower_fed - 0.005, 1.005, color="0.85", alpha=0.5, label=r"Bounds $\pm 0.5$ pp")899ax.set_ylim(0.95, 1.12); ax.set_ylabel(r"$\Sigma_t$")900ax.legend(loc="upper right", frameon=False, ncol=2)901ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))902save_fig(fig, "fig2_completeness_test")903904# fig3 RN density snapshots905fig, ax = plt.subplots(figsize=(6.5, 3.4))906idx = q.index907snap = [idx[0], idx[len(idx) // 4], idx[len(idx) // 2], idx[3 * len(idx) // 4], idx[-1]]908w = 3.5909for i, d0 in enumerate(snap):910 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)911ax.set_xticks(STATE_BPS); ax.set_xticklabels(["$\\leq-50$", "$-25$", "0", "$+25$", "$\\geq+50$"])912ax.set_xlabel("Change in the target range upper bound, June 2026 meeting (bps)"); ax.set_ylabel("Risk-neutral probability $q_{j,t}$")913ax.legend(frameon=False)914save_fig(fig, "fig3_risk_neutral_density")915916# fig4 E_Q path with release markers917fig, axes = plt.subplots(2, 1, figsize=(6.5, 5.2), sharex=True, gridspec_kw={"height_ratios": [2, 1.3]})918axes[0].plot(eq.index, eq.values, color=PAL[0], lw=1.2)919for lab, ds, mk, col in [("FOMC", fomc, "v", PAL[1]), ("CPI", cpi, "o", PAL[2]), ("Employment report", nfp, "s", PAL[3])]:920 axes[0].scatter(ds, eq.reindex(ds).values, marker=mk, s=28, color=col, zorder=3, label=lab)921axes[0].set_ylabel(r"$\mathbb{E}^{\mathbb{Q}}_t[\Delta r]$ (bps)"); axes[0].legend(frameon=False, loc="lower left")922axes[1].bar(d_eq.index, d_eq.values, width=1, color=np.where(d_eq.index.isin(all_ann), PAL[1], "0.6"))923axes[1].set_ylabel(r"$\Delta\mathbb{E}^{\mathbb{Q}}_t[\Delta r]$ (bps)")924axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))925save_fig(fig, "fig4_expected_policy_change")926927# fig5 announcement-day distributions928fig, axes = plt.subplots(1, 2, figsize=(6.5, 3.0))929groups = [("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)])]930for k, (lab, s) in enumerate(groups):931 jitter = (np.random.RandomState(k).rand(len(s)) - 0.5) * 0.3932 axes[0].scatter(np.full(len(s), k) + jitter, s.abs().values, s=10, color=PAL[k], alpha=0.8)933 axes[0].hlines(s.abs().mean(), k - 0.3, k + 0.3, color="black", lw=1.2)934axes[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)")935axes[0].set_title("A. Expected policy change", loc="left")936groups2 = [("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)])]937for k, (lab, s) in enumerate(groups2):938 jitter = (np.random.RandomState(k).rand(len(s)) - 0.5) * 0.3939 axes[1].scatter(np.full(len(s), k) + jitter, s.abs().values, s=10, color=PAL[k], alpha=0.8)940 axes[1].hlines(s.abs().mean(), k - 0.3, k + 0.3, color="black", lw=1.2)941axes[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)")942axes[1].set_title("B. ``No change'' state price", loc="left")943fig.tight_layout()944save_fig(fig, "fig5_event_study_cars")945946# fig7 Granger p-value heatmap (AIC lag)947mat = pd.DataFrame(np.nan, index=panel.columns, columns=panel.columns)948for x in panel.columns:949 for y in panel.columns:950 if x != y:951 mat.loc[y, x] = gc_pair(y, x, p_main)[1]952fig, ax = plt.subplots(figsize=(4.6, 3.8))953im = ax.imshow(mat.values.astype(float), cmap="Blues_r", vmin=0, vmax=0.5)954ax.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")955ax.set_yticks(range(5)); ax.set_yticklabels([names[c].replace("$\\Delta\\mathbb{E}^{\\mathbb{Q}}[\\Delta r]$", "Fed") for c in panel.columns])956ax.set_xlabel("Causing variable (lagged)"); ax.set_ylabel("Caused variable"); ax.grid(False)957for i in range(5):958 for j in range(5):959 v = mat.values[i, j]960 if not np.isnan(v):961 ax.text(j, i, f"{v:.2f}", ha="center", va="center", fontsize=8, color="white" if v < 0.2 else "black")962fig.colorbar(im, ax=ax, label=f"Granger $p$-value, VAR({p_main})", shrink=0.8)963save_fig(fig, "fig7_granger_heatmap")964STATS["granger"]["full_matrix_p"] = mat.astype(float).round(4).to_dict()965966# fig9 volatility967fig, axes = plt.subplots(2, 1, figsize=(6.5, 4.6), sharex=True)968axes[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")969for lab, col, ls_ in [("GARCH(1,1), normal", PAL[0], "-"), ("GJR-GARCH(1,1), Student-$t$", PAL[1], "--")]:970 cv = garch_fits[lab].conditional_volatility971 axes[1].plot(cv.index, cv.values, color=col, ls=ls_, label=lab)972axes[1].set_ylabel("Conditional SD (%)"); axes[1].legend(frameon=False); axes[1].set_title("B. Conditional volatility", loc="left")973axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))974save_fig(fig, "fig9_volatility")975976# fig10 liquidity: Roll spread + order-book depth977fig, axes = plt.subplots(1, 2, figsize=(6.5, 3.2))978axes[0].plot(roll["No change"].index, roll["No change"].values * 100, color=PAL[0], label="No change")979axes[0].plot(roll["25 bps decrease"].index, roll["25 bps decrease"].values * 100, color=PAL[1], ls="--", label="25 bps decrease")980axes[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")981axes[0].xaxis.set_major_formatter(mdates.DateFormatter("%b"))982for i, c in enumerate(STATES):983 path = DATA_RAW_POLY / f"order_book_{[m for m in fed_mk if m['name']==c][0]['yes_token'][:16]}.json"984 d = json.load(open(path))985 bids = sorted([(float(b["price"]), float(b["size"])) for b in d["bids"]], reverse=True)986 asks = sorted([(float(a["price"]), float(a["size"])) for a in d["asks"]])987 mid = ob[c]["mid"]988 xb = [(mid - p) * 100 for p, _ in bids]; yb = np.cumsum([p * s for p, s in bids]) / 1e3989 xa = [(p - mid) * 100 for p, _ in asks]; ya = np.cumsum([p * s for p, s in asks]) / 1e3990 axes[1].step(xb, yb, where="post", color=PAL[i], ls="--", lw=1)991 axes[1].step(xa, ya, where="post", color=PAL[i], ls="-", lw=1, label=c)992axes[1].set_xscale("log"); axes[1].set_xlim(0.05, 100); axes[1].set_yscale("log")993axes[1].set_xlabel("Distance from midpoint (cents)"); axes[1].set_ylabel("Cumulative depth (USD 000)")994axes[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")995fig.tight_layout()996save_fig(fig, "fig10_liquidity")997998# fig11 election999fig, axes = plt.subplots(2, 1, figsize=(6.5, 5.0), sharex=True, gridspec_kw={"height_ratios": [2, 1]})1000for i, c in enumerate(["Donald Trump", "Kamala Harris", "Joe Biden", "Michelle Obama", "Robert F. Kennedy Jr."]):1001 axes[0].plot(election.index, election[c], color=PAL[i], ls=LS[i], label=c)1002for lab, d0 in ELECTION_KEY_DATES.items():1003 axes[0].axvline(pd.Timestamp(d0), color="0.6", lw=0.6, ls=":")1004 axes[0].text(pd.Timestamp(d0), 0.98, lab, rotation=90, fontsize=6, va="top", ha="right", color="0.3")1005axes[0].set_ylabel("Price (USDC)"); axes[0].legend(frameon=False, ncol=2, loc="upper left", fontsize=7); axes[0].set_ylim(0, 1.0)1006axes[1].plot(sigma_el.index, sigma_el.values, color=PAL[0]); axes[1].axhline(1, color="black", lw=0.8)1007axes[1].plot(sigma_el.index, lower_el, color=PAL[1], ls="--", lw=0.9)1008axes[1].set_ylabel(r"$\Sigma_t$, 17 contracts"); axes[1].set_ylim(0.9, 1.1)1009axes[1].xaxis.set_major_locator(mdates.MonthLocator()); axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%b"))1010axes[1].set_xlabel("2024")1011save_fig(fig, "fig11_election_case_study")10121013# fig12 deviation distribution1014fig, axes = plt.subplots(1, 2, figsize=(6.5, 2.9))1015dev = sigma - 11016axes[0].hist(dev.values * 100, bins=30, color=PAL[0], edgecolor="white"); axes[0].axvline(0, color="black", lw=0.8)1017axes[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")1018axes[0].set_title("A. Fed June 2026", loc="left")1019dev_el = sigma_el - 11020axes[1].hist(dev_el.values * 100, bins=30, color=PAL[0], edgecolor="white"); axes[1].axvline(0, color="black", lw=0.8)1021axes[1].set_xlabel(r"$\Sigma_t-1$ (pp), election group"); axes[1].set_title("B. Election 2024", loc="left")1022fig.tight_layout()1023save_fig(fig, "fig12_arbitrage_violations")102410251026# ----------------------------------------------------------------------------1027# 11b. Sensitivity checks reported in Section 71028# ----------------------------------------------------------------------------1029SENS = {}1030# (a) announcement tests excluding the two largest FOMC days1031big2 = d_pnc.abs().nlargest(2).index1032mask_ex = ~d_pnc.index.isin(big2)1033c_ex = compare(d_pnc[mask_ex], is_ann[mask_ex], "All releases ex 2 largest FOMC days")1034c_ex_eq = compare(d_eq[mask_ex], is_ann[mask_ex], "All releases ex 2 largest FOMC days")1035SENS["ex_two_largest_days"] = {"days": [str(x.date()) for x in big2], "pnc": c_ex, "eq": c_ex_eq}1036# (b) cost allowance of 1 cent on the complete set1037tau_f = (RESOLUTION_FED - sigma.index).days.values / 365.01038r_f = dff.reindex(sigma.index).ffill().bfill().values / 100.01039low_f = np.exp(-r_f * tau_f)1040SENS["share_inside_band_1c"] = float(((sigma.values >= low_f - 0.01) & (sigma.values <= 1.01)).mean())1041SENS["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())1042# (f) lower bound with the target-range upper bound (3.75%) instead of DFF1043low_ub = np.exp(-0.0375 * tau_f)1044SENS["max_abs_diff_lower_bound_ub_vs_dff"] = float(np.abs(low_ub - low_f).max())1045# (e) ±75 valuation of open-ended states1046eq75 = pd.Series(q.values @ np.array([-75.0, -25.0, 0.0, 25.0, 75.0]), index=q.index)1047d_eq75 = eq75.diff().dropna()1048SENS["eq75_corr_with_eq50"] = float(np.corrcoef(d_eq75.values, d_eq.values)[0, 1])1049SENS["eq75_sd_ratio"] = float(d_eq75.std() / d_eq.std())1050SENS["eq75_all_releases"] = compare(d_eq75, is_ann, "All releases, ±75")1051# (c) Granger with Δp(No change) as the Polymarket variable1052panel_pnc = panel.copy(); panel_pnc["Fed"] = (pnc_grid.diff() * 100).reindex(panel.index)1053panel_pnc = panel_pnc.dropna()1054def gc_pair_df(df, y, x, p):1055 r = VAR(df[[y, x]]).fit(p); t = r.test_causality(y, [x], kind="f"); return float(t.test_statistic), float(t.pvalue)1056SENS["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}1057# (d) half-sample Granger at p=11058h = len(panel) // 21059SENS["granger_halves_p1"] = {}1060for lab, df in [("first_half", panel.iloc[:h]), ("second_half", panel.iloc[h:])]:1061 SENS["granger_halves_p1"][lab] = {f"{x}->{y}": gc_pair_df(df, y, x, 1) for x, y in pairs}1062 SENS["granger_halves_p1"][lab]["T"] = int(len(df))1063# (g) election facts1064SENS["trump_max_date"] = str(tr.idxmax().date()); SENS["trump_max"] = float(tr.max())1065SENS["nochange_mar17"] = float(fed_main["No change"].loc["2026-03-17"]); SENS["nochange_mar18"] = float(fed_main["No change"].loc["2026-03-18"])1066SENS["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"])1067SENS["n_release_days"] = int(is_ann.sum()); SENS["n_other_days"] = int((~is_ann).sum())1068STATS["sensitivity"] = SENS10691070# ----------------------------------------------------------------------------1071# 12. Save stats1072# ----------------------------------------------------------------------------1073def _clean(o):1074 if isinstance(o, dict):1075 return {str(k): _clean(v) for k, v in o.items()}1076 if isinstance(o, (list, tuple)):1077 return [_clean(v) for v in o]1078 if isinstance(o, (np.floating, float)):1079 return None if np.isnan(o) else float(o)1080 if isinstance(o, (np.integer,)):1081 return int(o)1082 if isinstance(o, (pd.Timestamp,)):1083 return str(o)1084 return o108510861087(ROOT / "results" / "revision_stats.json").write_text(json.dumps(_clean(STATS), indent=2))1088fed_main.assign(state_price_sum=sigma, EQ_bps=eq).to_parquet(ROOT / "data" / "processed" / "fed_june2026_state_prices_aligned.parquet")1089election.to_parquet(ROOT / "data" / "processed" / "election_2024_panel_aligned.parquet")1090panel.to_parquet(ROOT / "data" / "processed" / "trading_day_panel.parquet")1091print(json.dumps(_clean({k: STATS[k] for k in ["fed", "completeness", "orderbook", "announcement", "granger", "garch", "har", "roll", "election", "events"]}), indent=1)[:20000])1092