WP40 RAID: manuscript polish — unified math notation and aligned definitions, TikZ architecture, nine journal-style figures (paper_style.py, make_paper_figures.py)
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
24 changed files +464 −61
modified
docs/PROGRESS.md
+14 −0
@@ -210,3 +210,17 @@ Known cosmetic overfull boxes in two listings; no placeholder, no hand-typed num | ||
| 210 | 210 | `key_numbers.json`. Secret/data scan: `scripts/check_no_secrets.sh` clean; only synthetic fixtures tracked. |
| 211 | 211 | - Release checklist: `docs/RELEASE-CHECKLIST.md` (gaps stated: Linux build not executed, pip wheel not built from a |
| 212 | 212 | clean environment, E1–E6 need a FinData key). Tag `v0.1.0-wp40`. |
| 213 | + | |
| 214 | +## Phase 7b — Manuscript polish (2026-11-02, 22:30–23:30 ET) | |
| 215 | +Owner request: nicer mathematics and journal-quality figures. Done: unified notation (`\DeclareMathOperator` for RV, BPV, | |
| 216 | +RQ, RSV, TriPV, RCov, Var; bold vectors/matrices; indicator and absolute-value macros), `mathtools`/`bm`; estimator | |
| 217 | +section rewritten with aligned, numbered definitions (returns, realized measures, OLS/HC/HAC/cluster covariances, HAR, | |
| 218 | +HARQ, QLIKE, GARCH/GJR/EGARCH, QML log-likelihood, QMLE sandwich, forecasts, bootstrap, VAR/GFEVD/DY, Ledoit–Wolf); | |
| 219 | +Metal section with the TwoSum error-free transformation and the input-rounding bound as equations; architecture | |
| 220 | +diagram redrawn in TikZ. Figures: new module `experiments/scripts/paper_style.py` (serif/STIX, thin spines, light | |
| 221 | +grids, panel labels "Panel A.", grayscale-safe palette, 300 dpi vector PDF) and `make_paper_figures.py` producing nine | |
| 222 | +figures from the artifacts — primitives (throughput + speedup with bootstrap CIs), E1 (throughput, speedup heat map, | |
| 223 | +fusion gain), precision vs segment length with the theoretical bound, E6 (stage breakdown + engines with CIs), E7 | |
| 224 | +(speedup vs GPU cores + cold/warm compile), E2 (QLIKE distributions + runtime), E3 (convergence, persistence, QLIKE), | |
| 225 | +E5 (coverage, bias/RMSE), E4 (conditioning, shrinkage, connectedness). Manuscript 36 pages, 0 undefined, no large | |
| 226 | +overfull boxes; figure hashes appended to `experiments/manifests/paper_assets.json`. | |
modified
experiments/manifests/paper_assets.json
+9 −2
@@ -14,7 +14,7 @@ | ||
| 14 | 14 | "tables/tab_primitives.tex": "dfccbd67b68adc9c451a2f7f196a70de42cb612ecc9826285e2bc7fadb1db0ff", |
| 15 | 15 | "figures/fig_primitives_throughput.pdf": "e9ce50655262712762064350c514e8dfd09fc2fce5c8db2027ed0e74b32b09eb", |
| 16 | 16 | "tables/tab_e1.tex": "c5a0de42b79895abd2dc0cffba23b92f9a2048093579fa3529335e24c1a90a97", |
| 17 | − "figures/fig_e1.pdf": "53e26cdf3689df17adaebf18f1f95441f2bac33685ffe6e496efe223014a521d", | |
| 17 | + "figures/fig_e1.pdf": "4680c8db79a8c66c306fc78991d17d55f6fd2ec8398064d7b947cb42e2486bf3", | |
| 18 | 18 | "tables/tab_e2_accuracy.tex": "a337d1e20c1b1298f83e55984296b06ebfc41834c723e0270236e78022bb6f65", |
| 19 | 19 | "tables/tab_e2_runtime.tex": "c43af1813c748886e5dc5c5f26257d3b14a4cc1dab5ea9cee5672e996f739abf", |
| 20 | 20 | "tables/tab_e3_fits.tex": "d5d6adfb8f9644524e5deac5a37ce775df6cf32b7df9cbcad3ffb15a1a63587e", |
@@ -22,7 +22,14 @@ | ||
| 22 | 22 | "tables/tab_e5.tex": "69b53b60302201eed7a2ffba95827f88fa44b7ba8b3e642f820ee696652495d7", |
| 23 | 23 | "tables/tab_e6.tex": "d06f84a429329af2a4f1e7831cff4cd4371a4189806fea8ba226e73890c47aee", |
| 24 | 24 | "tables/tab_e7.tex": "26125ba7860235e8ec9bc7cbf9094525ae3badfa7684b0fdeceff0c444db27bb", |
| 25 | − "figures/fig_e7.pdf": "9e32e7a7c8d4fb0e1688ba282d05650b48cc53349a14aee834f4b19d5b306819" | |
| 25 | + "figures/fig_e7.pdf": "d4e0e29cf24f331d1de4ed0a757a36c4725191c1c8f424e5af60e393821c0890", | |
| 26 | + "figures/fig_primitives.pdf": "89297ebe238fd09e1686350b3fb3385cf0a29c8674999adb19c73524463afb91", | |
| 27 | + "figures/fig_precision.pdf": "26daf13f3f7107990eda67cb46055c11ea34b479a6b1a86eb76e77edbc6aea0f", | |
| 28 | + "figures/fig_e6.pdf": "51ecf937b9e6a527d867c9cdee24c4a44f31e3f6e007b39e3e1907cbfdd85f6b", | |
| 29 | + "figures/fig_e2.pdf": "09364703b2496bd14c6d08a2fa2bc79aa039eebde9627ae2531cabc9db03686b", | |
| 30 | + "figures/fig_e3.pdf": "d5d1f28b390703a45da27ebab1e8c233c547c35acdbe7d8d319c6fc95f3d6f42", | |
| 31 | + "figures/fig_e5.pdf": "6bac89369d81bf1ab8f461b4a59c792bf49f34c7c278858fb60b5d48073820a5", | |
| 32 | + "figures/fig_e4.pdf": "79bce2a904cb37f9b222bc4eee0e87b2b9ff275556b47de2673e653e2f3b78dd" | |
| 26 | 33 | }, |
| 27 | 34 | "key_numbers_sha256": "4424672bebb9b3b34b9c9465c2afa92a10901aedfd967a1d7dc843855a0db242" |
| 28 | 35 | } |
| \ No newline at end of file | ||
added
experiments/scripts/make_paper_figures.py
+181 −0
@@ -0,0 +1,181 @@ | ||
| 1 | +#!/usr/bin/env python3 | |
| 2 | +"""Journal-style figures for the manuscript, generated from the tracked artifacts. Writes paper/figures/fig_*.pdf and | |
| 3 | +appends their SHA-256 to experiments/manifests/paper_assets.json. Run after make_paper_assets.py.""" | |
| 4 | +import os, sys, json, collections, hashlib, statistics | |
| 5 | +import numpy as np | |
| 6 | +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) | |
| 7 | +import paper_style as ps; ps.setup() | |
| 8 | +import matplotlib.pyplot as plt | |
| 9 | +from matplotlib.ticker import FuncFormatter, LogLocator | |
| 10 | +ROOT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", ".."); RES = os.path.join(ROOT, "experiments", "results"); RAW = os.path.join(RES, "raw") | |
| 11 | +FIG = os.path.join(ROOT, "paper", "figures"); os.makedirs(FIG, exist_ok=True) | |
| 12 | +key = json.load(open(os.path.join(RES, "key_numbers.json"))); rng = np.random.default_rng(20261102); outputs = {} | |
| 13 | +def sha(p): return hashlib.sha256(open(p, "rb").read()).hexdigest() | |
| 14 | +def load(name): | |
| 15 | + p = os.path.join(RAW, name); return [json.loads(l) for l in open(p) if l.strip()] if os.path.exists(p) else [] | |
| 16 | +def ci(x, n=2000): x = np.asarray(x, float); m = [np.median(rng.choice(x, x.size, replace=True)) for _ in range(n)]; return float(np.median(x)), float(np.percentile(m, 2.5)), float(np.percentile(m, 97.5)) | |
| 17 | +def save(fig, name): p = os.path.join(FIG, name); fig.savefig(p, bbox_inches="tight", metadata={"CreationDate": None}); plt.close(fig); outputs["figures/" + name] = sha(p); print("wrote", name) | |
| 18 | +W = 6.3 # text width in inches (letter, 1in margins) | |
| 19 | + | |
| 20 | +# ---------------------------------------------------------------- Fig 1: primitives — throughput (A) and speedup with CIs (B) | |
| 21 | +prim = [r for r in load("bench_primitives.jsonl") if r.get("kernel") != "shader_library_load"] | |
| 22 | +g = collections.defaultdict(list) | |
| 23 | +for r in prim: g[(r["kernel"], r["S"], r["n"], r["backend"], r["threads"] == 1)].append(r["wall_seconds"]) | |
| 24 | +med = {k: ci(v) for k, v in g.items()} | |
| 25 | +kernels = [("realized_fused", "Realized"), ("regression_stats", "Regression"), ("rolling_sumsq", "Rolling sum"), ("garch_nll_P64", "GARCH NLL")] | |
| 26 | +fig, axes = plt.subplots(2, 4, figsize=(W, 4.2), gridspec_kw={"hspace": 0.6, "wspace": 0.42}) | |
| 27 | +for j, (k, title) in enumerate(kernels): | |
| 28 | + ax = axes[0, j] | |
| 29 | + for be, one in [("cpu", True), ("cpu", False), ("metal", False)]: | |
| 30 | + tag = "cpu1" if one else ("cpuN" if be == "cpu" else "metal") | |
| 31 | + pts = sorted([(S * n, S * n / med[(kk, S, n, b, o)][0]) for (kk, S, n, b, o) in med if kk == k and b == be and o == one]) | |
| 32 | + if pts: ax.plot([p[0] for p in pts], [p[1] for p in pts], marker=ps.MARK[tag], ms=2.6, lw=0.8, color=ps.PALETTE[tag], label=ps.LABEL[tag], markerfacecolor="white" if tag == "cpu1" else None, markeredgewidth=0.6) | |
| 33 | + ax.set_xscale("log"); ax.set_yscale("log"); ps.panel(ax, "ABCD"[j], title); ps.light_grid(ax) | |
| 34 | + if j == 0: ax.set_ylabel("Elements per second") | |
| 35 | + ax = axes[1, j] | |
| 36 | + pts = [] | |
| 37 | + for (kk, S, n, b, o) in med: | |
| 38 | + if kk == k and b == "metal": | |
| 39 | + c = med.get((k, S, n, "cpu", False)); m = med[(kk, S, n, b, o)] | |
| 40 | + if c: pts.append((S * n, c[0] / m[0], c[0] / m[2], c[0] / m[1])) | |
| 41 | + pts.sort() | |
| 42 | + if pts: | |
| 43 | + x = [p[0] for p in pts]; y = [p[1] for p in pts]; lo = [p[2] for p in pts]; hi = [p[3] for p in pts] | |
| 44 | + ax.errorbar(x, y, yerr=[np.subtract(y, lo), np.subtract(hi, y)], fmt="^", ms=3, color=ps.PALETTE["metal"], ecolor="#b2182b", elinewidth=0.6, capsize=1.5, lw=0) | |
| 45 | + ax.set_xscale("log"); ax.set_yscale("log"); ps.ref_line(ax, 1.0, None); ps.light_grid(ax); ax.set_xlabel("Elements") | |
| 46 | + if j == 0: ax.set_ylabel("Speedup, Metal / CPU 18 threads") | |
| 47 | + ps.panel(ax, "EFGH"[j], "Speedup") | |
| 48 | +h, l = axes[0, 0].get_legend_handles_labels(); fig.legend(h, l, loc="lower center", ncol=3, bbox_to_anchor=(0.5, -0.04), fontsize=7) | |
| 49 | +save(fig, "fig_primitives.pdf") | |
| 50 | + | |
| 51 | +# ---------------------------------------------------------------- Fig 2: E1 — (A) throughput real, (B) speedup grid heat, (C) fusion gain | |
| 52 | +e1 = [r for r in load("e1_batched_rv.jsonl") if r.get("kind") == "bench"] | |
| 53 | +g1 = collections.defaultdict(list) | |
| 54 | +for r in e1: g1[(r["source"], r["S"], r["segment"], r["backend"])].append(r["wall_seconds"]) | |
| 55 | +m1 = {k: ci(v) for k, v in g1.items()} | |
| 56 | +segs = ["day", "week", "month", "year"]; Ss = sorted({r["S"] for r in e1 if r["source"] == "real"}) | |
| 57 | +fig, axes = plt.subplots(1, 3, figsize=(W, 2.4), gridspec_kw={"wspace": 0.55, "width_ratios": [1.4, 1, 1]}) | |
| 58 | +ax = axes[0] | |
| 59 | +for tag, be in [("cpu1", "cpu1"), ("cpuN", "cpuN"), ("metal", "metal"), ("metal_unfused", "metal_unfused")]: | |
| 60 | + pts = sorted([(next(r["elements"] for r in e1 if (r["source"], r["S"], r["segment"], r["backend"]) == ("real", S, s, be)), next(r["elements"] for r in e1 if (r["source"], r["S"], r["segment"], r["backend"]) == ("real", S, s, be)) / m1[("real", S, s, be)][0]) for S in Ss for s in segs if ("real", S, s, be) in m1]) | |
| 61 | + ax.plot([p[0] for p in pts], [p[1] for p in pts], ls="none", marker=ps.MARK[tag], ms=3, color=ps.PALETTE[tag], label=ps.LABEL[tag], markerfacecolor="white" if tag == "cpu1" else None, markeredgewidth=0.6) | |
| 62 | +ax.set_xscale("log"); ax.set_yscale("log"); ps.light_grid(ax); ax.set_xlabel("Returns in the batch"); ax.set_ylabel("Returns per second"); ps.panel(ax, "A", "Throughput, real panel"); ax.legend(fontsize=6, loc="lower right", handletextpad=0.3) | |
| 63 | +ax = axes[1]; M = np.array([[m1[("real", S, s, "cpuN")][0] / m1[("real", S, s, "metal")][0] for s in segs] for S in Ss]) | |
| 64 | +im = ax.imshow(M, cmap="RdGy_r", vmin=0.6, vmax=1.4, aspect="auto"); ax.set_xticks(range(4)); ax.set_xticklabels(segs); ax.set_yticks(range(len(Ss))); ax.set_yticklabels([str(S) for S in Ss]); ax.set_xlabel("Segment length"); ax.set_ylabel("Panel width $S$") | |
| 65 | +for i in range(len(Ss)): | |
| 66 | + for j in range(4): ax.text(j, i, f"{M[i, j]:.2f}", ha="center", va="center", fontsize=6.5, color="white" if abs(M[i, j] - 1) > 0.25 else "black") | |
| 67 | +ps.panel(ax, "B", "Speedup vs. CPU"); ax.spines["top"].set_visible(True); ax.spines["right"].set_visible(True) | |
| 68 | +ax = axes[2]; F = np.array([[m1[("real", S, s, "metal_unfused")][0] / m1[("real", S, s, "metal")][0] for s in segs] for S in Ss]) | |
| 69 | +im2 = ax.imshow(F, cmap="Greys", vmin=1.0, vmax=3.5, aspect="auto"); ax.set_xticks(range(4)); ax.set_xticklabels(segs); ax.set_yticks(range(len(Ss))); ax.set_yticklabels([str(S) for S in Ss]); ax.set_xlabel("Segment length") | |
| 70 | +for i in range(len(Ss)): | |
| 71 | + for j in range(4): ax.text(j, i, f"{F[i, j]:.2f}", ha="center", va="center", fontsize=6.5, color="white" if F[i, j] > 2.2 else "black") | |
| 72 | +ps.panel(ax, "C", "Fusion gain"); ax.spines["top"].set_visible(True); ax.spines["right"].set_visible(True) | |
| 73 | +save(fig, "fig_e1.pdf") | |
| 74 | + | |
| 75 | +# ---------------------------------------------------------------- Fig 3: precision — Metal vs FP64 relative error vs segment length (E1) + theory | |
| 76 | +fig, ax = plt.subplots(figsize=(W * 0.6, 2.3)) | |
| 77 | +pts = collections.defaultdict(list) | |
| 78 | +for r in e1: | |
| 79 | + if r["backend"] == "metal" and r["source"] == "real": pts[r["segment"]].append((r["elements"] / r["segments"], r["max_rel_err_rv_vs_fp64"])) | |
| 80 | +xs, ys = [], [] | |
| 81 | +for s in segs: | |
| 82 | + for M_, e in pts[s]: xs.append(M_); ys.append(e) | |
| 83 | +ax.scatter(xs, ys, s=10, color=ps.PALETTE["metal"], marker="^", zorder=3, label="Observed maximum over segments (E1)") | |
| 84 | +Mgrid = np.logspace(2.4, 5.1, 50); u = 2 ** -24 | |
| 85 | +ax.plot(Mgrid, 2 * u * np.sqrt(3 / Mgrid), color="#4d4d4d", lw=0.8, ls=(0, (4, 2)), label="RMS prediction $2u\\sqrt{3/M}$") | |
| 86 | +ax.plot(Mgrid, 3.2 * 2 * u * np.sqrt(3 / Mgrid), color="#4d4d4d", lw=0.8, ls=(0, (1, 1.5)), label="Expected maximum over 926 segments ($3.2\\sigma$)") | |
| 87 | +ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlabel("Returns per segment $M$"); ax.set_ylabel("Relative error"); ps.light_grid(ax, "both"); ax.legend(fontsize=6.5, loc="upper right"); ax.set_title("Metal vs. FP64 realized variance", loc="left", fontsize=8.5, fontweight="bold", pad=6) | |
| 88 | +save(fig, "fig_precision.pdf") | |
| 89 | + | |
| 90 | +# ---------------------------------------------------------------- Fig 4: E6 — stacked stages (A) and engines with CIs (B) | |
| 91 | +e6 = load("e6_pipeline.jsonl"); e6b = load("e6_baselines.jsonl") | |
| 92 | +fig, axes = plt.subplots(1, 2, figsize=(W, 2.3), gridspec_kw={"wspace": 0.4, "width_ratios": [1, 1.2]}) | |
| 93 | +ax = axes[0]; stages = ["read_seconds", "prep_seconds", "batch_seconds", "estimate_seconds"]; names = ["Parquet read", "Session filter + returns", "Batch build", "Estimators"]; cols = ["#1a1a1a", "#4d4d4d", "#999999", "#b2182b"] | |
| 94 | +for i, be in enumerate(["cpu", "metal"]): | |
| 95 | + rs = [r for r in e6 if r["backend_used"] == be]; bottom = 0 | |
| 96 | + for st, nm, c in zip(stages, names, cols): | |
| 97 | + v = statistics.median(r[st] for r in rs); ax.bar(i, v, bottom=bottom, color=c, width=0.6, label=nm if i == 0 else None, edgecolor="white", linewidth=0.4); bottom += v | |
| 98 | +ax.set_xticks([0, 1]); ax.set_xticklabels(["RAID, CPU", "RAID, Metal"]); ax.set_xlim(-0.6, 1.6); ax.set_ylabel("Seconds per file (2.4M bars)"); ps.light_grid(ax); ax.legend(fontsize=6.5, loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2); ps.panel(ax, "A", "Stage breakdown") | |
| 99 | +ax = axes[1]; eng = [("RAID, CPU", [r["total_seconds"] for r in e6 if r["backend_used"] == "cpu"], ps.PALETTE["raid"]), ("RAID, Metal", [r["total_seconds"] for r in e6 if r["backend_used"] == "metal"], "#ef8a62"), ("Polars", [r["total_seconds"] for r in e6b if r["engine"] == "polars"], ps.PALETTE["polars"]), ("DuckDB", [r["total_seconds"] for r in e6b if r["engine"] == "duckdb"], ps.PALETTE["duckdb"]), ("pandas", [r["total_seconds"] for r in e6b if r["engine"] == "pandas"], ps.PALETTE["pandas"])] | |
| 100 | +for i, (nm, v, c) in enumerate(eng): | |
| 101 | + m, lo, hi = ci(v); ax.barh(i, m, color=c, height=0.6, edgecolor="white", linewidth=0.4); ax.errorbar(m, i, xerr=[[m - lo], [hi - m]], fmt="none", ecolor="#1a1a1a", elinewidth=0.6, capsize=2); ax.text(m * 1.08, i, f"{m:.3f} s", va="center", fontsize=7) | |
| 102 | +ax.set_yticks(range(len(eng))); ax.set_yticklabels([e[0] for e in eng]); ax.invert_yaxis(); ax.set_xscale("log"); ax.set_xlabel("Seconds per file (median, 95% CI)"); ax.set_xlim(0.05, 6); ps.light_grid(ax, "x"); ps.panel(ax, "B", "End-to-end engines") | |
| 103 | +save(fig, "fig_e6.pdf") | |
| 104 | + | |
| 105 | +# ---------------------------------------------------------------- Fig 5: E7 — speedup vs GPU cores (A) and cold/warm compile (B) | |
| 106 | +e7 = key["E7"]; order = sorted(e7, key=lambda n: e7[n]["gpu_cores"]) | |
| 107 | +fig, axes = plt.subplots(1, 2, figsize=(W, 2.6), gridspec_kw={"wspace": 0.45, "width_ratios": [1.4, 1]}) | |
| 108 | +ax = axes[0] | |
| 109 | +for k, lab, c, mk in [("realized_fused@1000x3900_speedup_MT", "Fused realized, $1000\\times3900$", ps.PALETTE["metal"], "^"), ("garch_nll_P64@1000x390_speedup_MT", "GARCH likelihood, $P=64$", "#1a1a1a", "o"), ("regression_stats@1000x3900_speedup_MT", "Regression statistics", "#2166ac", "s"), ("rolling_sumsq@1000x3900_speedup_MT", "Rolling sum of squares", "#999999", "v")]: | |
| 110 | + ax.scatter([e7[n]["gpu_cores"] for n in order], [e7[n][k] for n in order], s=14, color=c, marker=mk, label=lab, zorder=3) | |
| 111 | +for n in order: ax.annotate(e7[n]["chip"].replace("Apple ", ""), (e7[n]["gpu_cores"], e7[n]["realized_fused@1000x3900_speedup_MT"]), fontsize=5.5, xytext=(3, 3), textcoords="offset points", color="#4d4d4d") | |
| 112 | +ps.ref_line(ax, 1.0, None); ax.set_yscale("log"); ax.set_xlabel("GPU cores"); ax.set_ylabel("Speedup, Metal / CPU all threads"); ax.set_ylim(0.12, 12); ps.light_grid(ax, "y"); ax.legend(fontsize=6.3, loc="upper center", bbox_to_anchor=(0.5, -0.28), ncol=2); ps.panel(ax, "A", "Eleven Apple Silicon machines") | |
| 113 | +ax = axes[1]; y = np.arange(len(order)) | |
| 114 | +ax.barh(y - 0.18, [1e3 * e7[n]["cold_compile_s"] for n in order], height=0.36, color="#1a1a1a", label="Cold (salted source)"); ax.barh(y + 0.18, [1e3 * e7[n]["warm_compile_s"] for n in order], height=0.36, color="#999999", label="Warm (shader cache)") | |
| 115 | +ax.set_yticks(y); ax.set_yticklabels([e7[n]["chip"].replace("Apple ", "") for n in order], fontsize=6.5); ax.set_xscale("log"); ax.set_xlabel("Runtime MSL compilation (ms)"); ax.invert_yaxis(); ps.light_grid(ax, "x"); ax.legend(fontsize=6.3, loc="lower right"); ps.panel(ax, "B", "Shader compilation") | |
| 116 | +save(fig, "fig_e7.pdf") | |
| 117 | + | |
| 118 | +# ---------------------------------------------------------------- Fig 6: E2 — QLIKE distributions by horizon (A) and runtime (B) | |
| 119 | +e2 = load("e2_har.jsonl"); acc = [r for r in e2 if r.get("kind") == "accuracy"]; rt = [r for r in e2 if r.get("kind") == "runtime"] | |
| 120 | +fig, axes = plt.subplots(1, 2, figsize=(W, 2.4), gridspec_kw={"wspace": 0.4, "width_ratios": [1.5, 1]}) | |
| 121 | +ax = axes[0]; pos = 0; ticks = []; labels = [] | |
| 122 | +for h in [1, 5, 22]: | |
| 123 | + rs = [r for r in acc if r["h"] == h] | |
| 124 | + for m, nm, c in [("random_walk", "Random walk", "#999999"), ("har", "HAR-RV", "#1a1a1a"), ("loghar", "log-HAR", ps.PALETTE["metal"]), ("harq", "HARQ", "#2166ac")]: | |
| 125 | + v = [r[m]["qlike"] for r in rs if "qlike" in r.get(m, {})] | |
| 126 | + bp = ax.boxplot(v, positions=[pos], widths=0.6, showfliers=False, patch_artist=True, medianprops={"color": "white", "lw": 1}, boxprops={"facecolor": c, "edgecolor": c, "lw": 0.5}, whiskerprops={"color": c, "lw": 0.6}, capprops={"color": c, "lw": 0.6}) | |
| 127 | + if h == 1: bp["boxes"][0].set_label(nm) | |
| 128 | + pos += 1 | |
| 129 | + ticks.append(pos - 2.5); labels.append(f"$h={h}$"); pos += 1 | |
| 130 | +ax.set_xticks(ticks); ax.set_xticklabels(labels); ax.set_ylabel("QLIKE across 100 symbols"); ps.light_grid(ax); ax.legend(fontsize=6.5, loc="upper left", ncol=2); ps.panel(ax, "A", "Rolling causal forecasts, 2024") | |
| 131 | +ax = axes[1]; g2 = collections.defaultdict(list) | |
| 132 | +for r in rt: g2[(r["S"], r["mode"])].append(r["wall_seconds"]) | |
| 133 | +Ss2 = sorted({r["S"] for r in rt}) | |
| 134 | +for m, nm, c, mk in [("per_series_cpu1", "Per-series QR, 1 thread", "#8c8c8c", "s"), ("per_series_cpuMT", "Per-series QR, 18 threads", "#1a1a1a", "o"), ("batched_cpu", "Batched statistics, CPU", "#2166ac", "D"), ("batched_metal", "Batched statistics, Metal", ps.PALETTE["metal"], "^")]: | |
| 135 | + ax.plot(Ss2, [1e3 * ci(g2[(S, m)])[0] for S in Ss2], marker=mk, ms=3, color=c, label=nm, markerfacecolor="white" if m == "per_series_cpu1" else None, markeredgewidth=0.6) | |
| 136 | +ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlabel("Series ($T=252$ each)"); ax.set_ylabel("Milliseconds"); ps.light_grid(ax, "both"); ax.legend(fontsize=6, loc="upper left"); ps.panel(ax, "B", "HAR estimation runtime") | |
| 137 | +save(fig, "fig_e2.pdf") | |
| 138 | + | |
| 139 | +# ---------------------------------------------------------------- Fig 7: E3 — convergence (A), persistence distribution (B), QLIKE (C) | |
| 140 | +e3 = load("e3_garch.jsonl"); fits = [r for r in e3 if r.get("kind") == "fit"] | |
| 141 | +models = [("garch.n", "GARCH"), ("gjr.n", "GJR"), ("egarch.n", "EGARCH"), ("garch.t", "GARCH-$t$")] | |
| 142 | +fig, axes = plt.subplots(1, 3, figsize=(W, 2.3), gridspec_kw={"wspace": 0.5}) | |
| 143 | +ax = axes[0] | |
| 144 | +for i, (m, nm) in enumerate(models): | |
| 145 | + rs = [r for r in fits if r["model"] == m]; conv = 100 * np.mean([bool(r.get("converged")) for r in rs]); ax.bar(i, conv, color="#1a1a1a" if conv > 99 else ("#4d4d4d" if conv > 95 else ps.PALETTE["metal"]), width=0.6); ax.text(i, conv + 1, f"{conv:.1f}%", ha="center", fontsize=6.5) | |
| 146 | +ax.set_xticks(range(4)); ax.set_xticklabels([m[1] for m in models], rotation=25, ha="right"); ax.set_ylim(80, 103); ax.set_ylabel("Converged fits (%)"); ps.light_grid(ax); ps.panel(ax, "A", "Convergence, 888 series") | |
| 147 | +ax = axes[1] | |
| 148 | +for i, (m, nm) in enumerate(models): | |
| 149 | + v = [r["persistence"] for r in fits if r["model"] == m and r.get("converged") and r.get("persistence") is not None and 0 < r["persistence"] < 1.05] | |
| 150 | + ax.boxplot(v, positions=[i], widths=0.55, showfliers=False, patch_artist=True, medianprops={"color": "white", "lw": 1}, boxprops={"facecolor": "#4d4d4d", "edgecolor": "#4d4d4d", "lw": 0.5}, whiskerprops={"color": "#4d4d4d", "lw": 0.6}, capprops={"color": "#4d4d4d", "lw": 0.6}) | |
| 151 | +ax.set_xticks(range(4)); ax.set_xticklabels([m[1] for m in models], rotation=25, ha="right"); ax.set_ylabel("Persistence"); ps.light_grid(ax); ps.panel(ax, "B", "Persistence") | |
| 152 | +ax = axes[2] | |
| 153 | +for i, (m, nm) in enumerate(models): | |
| 154 | + v = [r["qlike_vs_rv"] for r in fits if r["model"] == m and r.get("converged") and r.get("qlike_vs_rv") is not None] | |
| 155 | + ax.boxplot(v, positions=[i], widths=0.55, showfliers=False, patch_artist=True, medianprops={"color": "white", "lw": 1}, boxprops={"facecolor": "#1a1a1a", "edgecolor": "#1a1a1a", "lw": 0.5}, whiskerprops={"color": "#1a1a1a", "lw": 0.6}, capprops={"color": "#1a1a1a", "lw": 0.6}) | |
| 156 | +rw = [r["qlike_rw_rv"] for r in fits if r["model"] == "garch.n" and r.get("qlike_rw_rv") is not None]; ax.axhline(np.median(rw), color=ps.PALETTE["metal"], lw=0.8, ls=(0, (4, 2))); ax.text(-0.45, np.median(rw) + 0.03, "RV random walk (median)", fontsize=6, color=ps.PALETTE["metal"], ha="left", va="bottom") | |
| 157 | +ax.set_xticks(range(4)); ax.set_xticklabels([m[1] for m in models], rotation=25, ha="right"); ax.set_ylabel("QLIKE vs. next-day RV"); ax.set_ylim(0, 1.5); ps.light_grid(ax); ps.panel(ax, "C", "One-day forecasts") | |
| 158 | +save(fig, "fig_e3.pdf") | |
| 159 | + | |
| 160 | +# ---------------------------------------------------------------- Fig 8: E5 — coverage (A) and bias/RMSE (B) | |
| 161 | +e5 = key["E5"]; fig, axes = plt.subplots(1, 2, figsize=(W, 2.5), gridspec_kw={"wspace": 0.45}) | |
| 162 | +ax = axes[0]; params = ["mu", "omega", "alpha", "beta"]; lab = ["$\\mu$", "$\\omega$", "$\\alpha$", "$\\beta$"]; x = np.arange(4) | |
| 163 | +ax.bar(x - 0.18, [e5["garch_n500"][p]["coverage95"] for p in params], width=0.36, color="#999999", label="GARCH(1,1), $n=500$"); ax.bar(x + 0.18, [e5["garch_n2000"][p]["coverage95"] for p in params], width=0.36, color="#1a1a1a", label="GARCH(1,1), $n=2000$") | |
| 164 | +ax.bar(np.arange(4) + 5 - 0.18, e5["har_h1"], width=0.36, color="#4393c3", label="HAR $h=1$ (HAC)"); ax.bar(np.arange(4) + 5 + 0.18, e5["har_h5"], width=0.36, color="#2166ac", label="HAR $h=5$ (HAC)") | |
| 165 | +ax.axhline(0.95, color=ps.PALETTE["metal"], lw=0.8, ls=(0, (4, 2))); ax.set_xticks(list(x) + list(np.arange(4) + 5)); ax.set_xticklabels(lab + ["$\\beta_0$", "$\\beta_d$", "$\\beta_w$", "$\\beta_m$"]); ax.set_ylim(0.8, 1.0); ax.set_ylabel("Empirical coverage of 95% intervals"); ps.light_grid(ax); ax.legend(fontsize=6, loc="upper center", bbox_to_anchor=(0.5, -0.22), ncol=2); ps.panel(ax, "A", "Coverage, $R=1000$") | |
| 166 | +ax = axes[1] | |
| 167 | +for i, (n, c) in enumerate([("garch_n500", "#999999"), ("garch_n2000", "#1a1a1a")]): | |
| 168 | + ax.errorbar(x + (i - 0.5) * 0.3, [e5[n][p]["bias"] for p in params], yerr=[e5[n][p]["rmse"] for p in params], fmt="o", ms=3, color=c, ecolor=c, elinewidth=0.7, capsize=2, label=f"$n={n[7:]}$: bias $\\pm$ RMSE") | |
| 169 | +ax.axhline(0, color="#4d4d4d", lw=0.6); ax.set_xticks(x); ax.set_xticklabels(lab); ax.set_ylabel("Bias (markers) and RMSE (bars)"); ps.light_grid(ax); ax.legend(fontsize=6.5); ps.panel(ax, "B", "GARCH(1,1) QML") | |
| 170 | +save(fig, "fig_e5.pdf") | |
| 171 | + | |
| 172 | +# ---------------------------------------------------------------- Fig 9: E4 — conditioning (A), shrinkage (B), connectedness and time (C) | |
| 173 | +e4 = key["E4"]; Ns = sorted(int(n) for n in e4) | |
| 174 | +fig, axes = plt.subplots(1, 3, figsize=(W, 2.2), gridspec_kw={"wspace": 0.6}) | |
| 175 | +ax = axes[0]; ax.plot(Ns, [min(e4[str(n)]["cond_median"], 1e9) for n in Ns], marker="o", ms=3, color="#1a1a1a", label="Median"); ax.plot(Ns, [min(e4[str(n)]["cond_max"], 1e9) for n in Ns], marker="s", ms=3, color="#999999", label="Maximum"); ax.set_yscale("log"); ax.set_xlabel("$N$ (assets)"); ax.set_ylabel("Condition number of daily RCov"); ax.axvline(78, color=ps.PALETTE["metal"], lw=0.7, ls=(0, (4, 2))); ax.text(80, 3e2, "$M=78$ bins", rotation=90, fontsize=6.5, color=ps.PALETTE["metal"], ha="left", va="bottom"); ps.light_grid(ax, "y"); ax.legend(fontsize=6.5, loc="upper left"); ps.panel(ax, "A", "Conditioning") | |
| 176 | +ax = axes[1]; ax.plot(Ns, [e4[str(n)]["lw"] for n in Ns], marker="^", ms=3, color="#2166ac"); ax.set_xlabel("$N$ (assets)"); ax.set_ylabel("Ledoit–Wolf shrinkage $\\delta^\\star$"); ps.light_grid(ax, "y"); ps.panel(ax, "B", "Shrinkage") | |
| 177 | +ax = axes[2]; ax.plot(Ns, [e4[str(n)]["dy_total"] for n in Ns], marker="o", ms=3, color="#1a1a1a", label="Total connectedness (%)"); ax.set_xlabel("$N$ (assets)"); ax.set_ylabel("Mean total connectedness (%)"); ps.light_grid(ax, "y"); ps.panel(ax, "C", "Connectedness") | |
| 178 | +for n in Ns: ax.annotate(f"{e4[str(n)]['dy_ms']:.1f} ms", (n, e4[str(n)]["dy_total"]), xytext=(5, 2) if n < 100 else (-28, -9), textcoords="offset points", fontsize=6, color="#4d4d4d") | |
| 179 | +save(fig, "fig_e4.pdf") | |
| 180 | + | |
| 181 | +mf = os.path.join(ROOT, "experiments", "manifests", "paper_assets.json"); man = json.load(open(mf)); man["outputs"].update(outputs); json.dump(man, open(mf, "w"), indent=1); print("manifest updated") | |
added
experiments/scripts/paper_style.py
+35 −0
@@ -0,0 +1,35 @@ | ||
| 1 | +"""Journal-style figure settings (serif text, Times-like math, thin spines, grayscale-safe palette, panel labels).""" | |
| 2 | +import matplotlib | |
| 3 | +matplotlib.use("Agg") | |
| 4 | +import matplotlib.pyplot as plt | |
| 5 | +from matplotlib import rcParams | |
| 6 | + | |
| 7 | +PALETTE = {"cpu1": "#8c8c8c", "cpuN": "#1a1a1a", "metal": "#b2182b", "metal_unfused": "#ef8a62", "pandas": "#999999", "polars": "#4393c3", "duckdb": "#2166ac", "raid": "#b2182b", "neutral": "#4d4d4d"} | |
| 8 | +MARK = {"cpu1": "s", "cpuN": "o", "metal": "^", "metal_unfused": "v"} | |
| 9 | +LABEL = {"cpu1": "CPU, 1 thread", "cpuN": "CPU, 18 threads", "metal": "Metal (fused)", "metal_unfused": "Metal (4 kernels)"} | |
| 10 | + | |
| 11 | +def setup(): | |
| 12 | + rcParams.update({ | |
| 13 | + "font.family": "serif", "font.serif": ["STIX Two Text", "Times New Roman", "Times", "DejaVu Serif"], | |
| 14 | + "mathtext.fontset": "stix", "font.size": 8.5, "axes.titlesize": 9, "axes.labelsize": 8.5, "legend.fontsize": 7.5, | |
| 15 | + "xtick.labelsize": 7.5, "ytick.labelsize": 7.5, "axes.linewidth": 0.6, "xtick.major.width": 0.5, "ytick.major.width": 0.5, | |
| 16 | + "xtick.major.size": 2.5, "ytick.major.size": 2.5, "xtick.direction": "out", "ytick.direction": "out", | |
| 17 | + "axes.spines.top": False, "axes.spines.right": False, "legend.frameon": False, "legend.handlelength": 1.8, | |
| 18 | + "lines.linewidth": 1.0, "lines.markersize": 3.5, "figure.dpi": 200, "savefig.dpi": 300, "pdf.fonttype": 42, | |
| 19 | + "axes.grid": False, "axes.prop_cycle": matplotlib.cycler(color=["#1a1a1a", "#b2182b", "#2166ac", "#4d4d4d", "#ef8a62"]), | |
| 20 | + }) | |
| 21 | + | |
| 22 | +def panel(ax, letter, title=None): | |
| 23 | + """Panel label in the Journal-of-Finance style: 'Panel A. Title' above the axes, left-aligned.""" | |
| 24 | + t = f"Panel {letter}." + (f" {title}" if title else "") | |
| 25 | + ax.set_title(t, loc="left", fontsize=8.5, fontweight="bold", pad=6) | |
| 26 | + | |
| 27 | +def light_grid(ax, axis="y"): | |
| 28 | + ax.grid(True, axis=axis, color="#d9d9d9", linewidth=0.4, linestyle="-", zorder=0); ax.set_axisbelow(True) | |
| 29 | + | |
| 30 | +def ref_line(ax, y=1.0, label=None): | |
| 31 | + ax.axhline(y, color="#4d4d4d", linewidth=0.6, linestyle=(0, (4, 2)), zorder=1) | |
| 32 | + if label: ax.annotate(label, xy=(1.0, y), xycoords=("axes fraction", "data"), xytext=(-2, 2), textcoords="offset points", ha="right", va="bottom", fontsize=7, color="#4d4d4d") | |
| 33 | + | |
| 34 | +def fmt_elements(ax): | |
| 35 | + ax.set_xscale("log"); ax.set_xlabel("Elements (returns)") | |
modified
paper/figures/fig_e1.pdf
+0 −0
Binary file not shown.
added
paper/figures/fig_e2.pdf
+0 −0
Binary file not shown.
added
paper/figures/fig_e3.pdf
+0 −0
Binary file not shown.
added
paper/figures/fig_e4.pdf
+0 −0
Binary file not shown.
added
paper/figures/fig_e5.pdf
+0 −0
Binary file not shown.
added
paper/figures/fig_e6.pdf
+0 −0
Binary file not shown.
modified
paper/figures/fig_e7.pdf
+0 −0
Binary file not shown.
added
paper/figures/fig_precision.pdf
+0 −0
Binary file not shown.
added
paper/figures/fig_primitives.pdf
+0 −0
Binary file not shown.
modified
paper/main.pdf
+0 −0
Binary file not shown.
modified
paper/main.tex
+30 −5
@@ -12,6 +12,10 @@ | ||
| 12 | 12 | \usepackage{setspace} |
| 13 | 13 | \onehalfspacing |
| 14 | 14 | \usepackage{amsmath,amssymb,amsthm} |
| 15 | +\usepackage{mathtools} | |
| 16 | +\usepackage{bm} | |
| 17 | +\usepackage{tikz} | |
| 18 | +\usetikzlibrary{positioning,arrows.meta,fit,backgrounds,calc} | |
| 15 | 19 | \usepackage{booktabs} |
| 16 | 20 | \usepackage{tabularx} |
| 17 | 21 | \usepackage{multirow} |
@@ -43,13 +47,34 @@ | ||
| 43 | 47 | \lstset{basicstyle=\ttfamily\footnotesize,breaklines=true,columns=fullflexible,frame=single,framerule=0.2pt} |
| 44 | 48 | \usepackage{url} |
| 45 | 49 | \newtheorem{proposition}{Proposition} |
| 46 | −\newcommand{\RV}{\mathrm{RV}} | |
| 47 | −\newcommand{\BPV}{\mathrm{BPV}} | |
| 48 | −\newcommand{\RQ}{\mathrm{RQ}} | |
| 49 | −\newcommand{\RSV}{\mathrm{RSV}} | |
| 50 | +% ---- notation ---- | |
| 51 | +\DeclareMathOperator{\RV}{RV} | |
| 52 | +\DeclareMathOperator{\BPV}{BPV} | |
| 53 | +\DeclareMathOperator{\RQ}{RQ} | |
| 54 | +\DeclareMathOperator{\RSV}{RSV} | |
| 55 | +\DeclareMathOperator{\TriPV}{TriPV} | |
| 56 | +\DeclareMathOperator{\RCov}{RCov} | |
| 57 | +\DeclareMathOperator{\Var}{Var} | |
| 58 | +\DeclareMathOperator{\diag}{diag} | |
| 59 | +\DeclareMathOperator{\tr}{tr} | |
| 60 | +\DeclareMathOperator*{\argmin}{arg\,min} | |
| 50 | 61 | \newcommand{\raid}{\textsc{raid}} |
| 51 | 62 | \newcommand{\E}{\mathbb{E}} |
| 52 | −\newcommand{\I}{\mathbb{1}} | |
| 63 | +\newcommand{\I}[1]{\mathbb{1}\{#1\}} | |
| 64 | +\newcommand{\abs}[1]{\lvert #1 \rvert} | |
| 65 | +\newcommand{\norm}[1]{\lVert #1 \rVert} | |
| 66 | +\newcommand{\bbeta}{\bm{\beta}} | |
| 67 | +\newcommand{\btheta}{\bm{\theta}} | |
| 68 | +\newcommand{\bSigma}{\bm{\Sigma}} | |
| 69 | +\newcommand{\bPhi}{\bm{\Phi}} | |
| 70 | +\newcommand{\bA}{\bm{A}} | |
| 71 | +\newcommand{\bX}{\bm{X}} | |
| 72 | +\newcommand{\by}{\bm{y}} | |
| 73 | +\newcommand{\bu}{\bm{u}} | |
| 74 | +\newcommand{\be}{\bm{e}} | |
| 75 | +\newcommand{\ulp}{\operatorname{ulp}} | |
| 76 | +\newcommand{\fl}{\operatorname{fl}} | |
| 77 | +\allowdisplaybreaks | |
| 53 | 78 | |
| 54 | 79 | % ============================================================================ |
| 55 | 80 | % METADATA |
modified
paper/sections/applications.tex
+21 −3
@@ -3,7 +3,13 @@ | ||
| 3 | 3 | |
| 4 | 4 | The applications are not forecasting studies; they show that the batched estimators deliver the same statistical content as per-series estimation on real data, and where batching pays. |
| 5 | 5 | |
| 6 | −\paragraph{E2 --- HAR forecasting and batched estimation.} For the first 100 symbols of the universe, HAR-RV, log-HAR and HARQ are re-estimated every day on a causal 120-session window and forecast the average realized variance over the next $h \in \{1, 5, 22\}$ sessions; the benchmark is the random-walk forecast $\widehat{\RV}_{t+h} = \RV_t$ (Table~\ref{tab:e2acc}). With only one year of data the windows are short, yet the median QLIKE of HAR-RV is 0.196, 0.146 and 0.134 at $h = 1, 5, 22$ against 0.251, 0.239 and 0.301 for the random walk, and HAR beats the random walk for 84\,\%, 86\,\% and 78\,\% of symbols; log-HAR is best at the longer horizons (0.138 and 0.112), consistent with the literature. Table~\ref{tab:e2rt} reports the runtime of estimating HAR-RV ($h = 1$, 252 daily observations per series) on 10, 100 and 781 series: per-series QR fits take 0.08--6.0\,ms on one thread and 0.05--0.68\,ms on 18 threads; the batched path (sufficient statistics on the backend, FP64 solve) takes 0.10--3.9\,ms on the CPU and 0.32--4.7\,ms on Metal. With 252 observations per series the designs are tiny ($252\times4$ per series), the GPU is launch-bound, and the per-series multithreaded QR is the fastest option; batching would pay only with long daily histories or many more series. Coefficients from the batched path agree with per-series QR to $10^{-9}$ (CPU) and $10^{-4}$ relative (Metal, FP32 sufficient statistics), the documented tolerance. | |
| 6 | +\paragraph{E2 --- HAR forecasting and batched estimation.} Figure~\ref{fig:e2} summarises the experiment. For the first 100 symbols of the universe, HAR-RV, log-HAR and HARQ are re-estimated every day on a causal 120-session window and forecast the average realized variance over the next $h \in \{1, 5, 22\}$ sessions; the benchmark is the random-walk forecast $\widehat{\RV}_{t+h} = \RV_t$ (Table~\ref{tab:e2acc}). With only one year of data the windows are short, yet the median QLIKE of HAR-RV is 0.196, 0.146 and 0.134 at $h = 1, 5, 22$ against 0.251, 0.239 and 0.301 for the random walk, and HAR beats the random walk for 84\,\%, 86\,\% and 78\,\% of symbols; log-HAR is best at the longer horizons (0.138 and 0.112), consistent with the literature. Table~\ref{tab:e2rt} reports the runtime of estimating HAR-RV ($h = 1$, 252 daily observations per series) on 10, 100 and 781 series: per-series QR fits take 0.08--6.0\,ms on one thread and 0.05--0.68\,ms on 18 threads; the batched path (sufficient statistics on the backend, FP64 solve) takes 0.10--3.9\,ms on the CPU and 0.32--4.7\,ms on Metal. With 252 observations per series the designs are tiny ($252\times4$ per series), the GPU is launch-bound, and the per-series multithreaded QR is the fastest option; batching would pay only with long daily histories or many more series. Coefficients from the batched path agree with per-series QR to $10^{-9}$ (CPU) and $10^{-4}$ relative (Metal, FP32 sufficient statistics), the documented tolerance. | |
| 7 | + | |
| 8 | +\begin{figure}[t] | |
| 9 | +\centering\includegraphics[width=\textwidth]{figures/fig_e2.pdf} | |
| 10 | +\caption{E2 --- HAR on the 2024 panel. Panel A: distribution across 100 symbols of the QLIKE of rolling causal forecasts (120-session window) for the random walk, HAR-RV, log-HAR and HARQ at horizons $h = 1, 5, 22$ (boxes: interquartile range; whiskers: 1.5 IQR; outliers omitted). Panel B: runtime of HAR-RV estimation on 10, 100 and 781 series of 252 observations by method (medians of 5 repetitions).} | |
| 11 | +\label{fig:e2} | |
| 12 | +\end{figure} | |
| 7 | 13 | |
| 8 | 14 | \begin{table}[t] |
| 9 | 15 | \centering\small |
@@ -18,7 +24,13 @@ The applications are not forecasting studies; they show that the batched estimat | ||
| 18 | 24 | \resizebox{\textwidth}{!}{\input{tables/tab_e2_runtime}} |
| 19 | 25 | \end{table} |
| 20 | 26 | |
| 21 | −\paragraph{E3 --- GARCH-family panel.} GARCH(1,1), GJR, EGARCH and GARCH-$t$ are fitted by QML to the daily close-to-close log returns of 888 symbols (2015--2024 daily bars, median 2\,459 training observations; the last 40 sessions of 2024 are held out). Table~\ref{tab:e3} reports convergence rates, iterations, time per fit, persistence, and the QLIKE of the fixed-parameter one-day-ahead variance forecasts against the next day's realized variance. GARCH and GARCH-$t$ converge for all 888 series, GJR for 98.4\,\% and EGARCH for 92.8\,\% (after up to two restarts; the failures are reported as such, not as estimates), in 2.9--18\,ms per fit; persistence medians are 0.956--0.973. Daily GARCH forecasts have a median QLIKE of 0.39--0.40 against 0.27--0.28 for the realized-variance random walk: a one-day-ahead forecast built from intraday data beats one built from daily returns, as the literature since \citet{AndersenBollerslevDieboldLabys2003} would predict, and the gap is not a numerical artefact. The GPU's contribution in E3 is the likelihood grid: evaluating 256 $(\omega,\alpha,\beta)$ combinations per series for all 888 series takes 23\,ms on Metal against 117\,ms on 18 CPU threads (5.1$\times$); the grid argmin is identical on both backends for all 888 series, with a maximum relative difference of $2\times10^{-3}$ in the likelihood at the least-likely grid points (near the non-stationarity boundary). | |
| 27 | +\paragraph{E3 --- GARCH-family panel.} Figure~\ref{fig:e3} summarises the experiment. GARCH(1,1), GJR, EGARCH and GARCH-$t$ are fitted by QML to the daily close-to-close log returns of 888 symbols (2015--2024 daily bars, median 2\,459 training observations; the last 40 sessions of 2024 are held out). Table~\ref{tab:e3} reports convergence rates, iterations, time per fit, persistence, and the QLIKE of the fixed-parameter one-day-ahead variance forecasts against the next day's realized variance. GARCH and GARCH-$t$ converge for all 888 series, GJR for 98.4\,\% and EGARCH for 92.8\,\% (after up to two restarts; the failures are reported as such, not as estimates), in 2.9--18\,ms per fit; persistence medians are 0.956--0.973. Daily GARCH forecasts have a median QLIKE of 0.39--0.40 against 0.27--0.28 for the realized-variance random walk: a one-day-ahead forecast built from intraday data beats one built from daily returns, as the literature since \citet{AndersenBollerslevDieboldLabys2003} would predict, and the gap is not a numerical artefact. The GPU's contribution in E3 is the likelihood grid: evaluating 256 $(\omega,\alpha,\beta)$ combinations per series for all 888 series takes 23\,ms on Metal against 117\,ms on 18 CPU threads (5.1$\times$); the grid argmin is identical on both backends for all 888 series, with a maximum relative difference of $2\times10^{-3}$ in the likelihood at the least-likely grid points (near the non-stationarity boundary). | |
| 28 | + | |
| 29 | +\begin{figure}[t] | |
| 30 | +\centering\includegraphics[width=\textwidth]{figures/fig_e3.pdf} | |
| 31 | +\caption{E3 --- GARCH-family QML on 888 daily return series (2015--2024). Panel A: share of converged fits by model. Panel B: distribution of the estimated persistence across converged fits. Panel C: distribution of the QLIKE of fixed-parameter one-day-ahead variance forecasts against next-day realized variance over the last 40 sessions of 2024; the dashed line is the median QLIKE of the realized-variance random walk.} | |
| 32 | +\label{fig:e3} | |
| 33 | +\end{figure} | |
| 22 | 34 | |
| 23 | 35 | \begin{table}[t] |
| 24 | 36 | \centering\small |
@@ -27,7 +39,13 @@ The applications are not forecasting studies; they show that the batched estimat | ||
| 27 | 39 | \resizebox{\textwidth}{!}{\input{tables/tab_e3_fits}} |
| 28 | 40 | \end{table} |
| 29 | 41 | |
| 30 | −\paragraph{E4 --- realized covariance and connectedness.} For $N \in \{10, 30, 60, 100\}$ symbols with at least 95\,\% regular-bar coverage, daily realized covariance matrices are formed from synchronised 5-minute returns (78 per day) on the days where every symbol has the full grid, and rolling Diebold--Yilmaz total connectedness is computed from a VAR(1) on the daily log-RV panel (120-day windows, $H = 10$). Table~\ref{tab:e4} shows the two constraints the charter asked for. Memory is not one of them: the rolling connectedness over 133 windows takes 0.5, 2.3, 8.9 and 29\,ms for $N = 10, 30, 60, 100$ on the CPU and never fails. Conditioning is: the median condition number of the daily realized covariance grows from 81 ($N = 10$) to 422 ($N = 30$) and $5.9\times10^3$ ($N = 60$), and the matrix is singular for $N = 100 > 78$ bins; the Ledoit--Wolf shrinkage intensity rises from 0.16 to 0.33 accordingly. Strict synchronisation also discards 70--164 of the 252 days as $N$ grows. Connectedness itself is not GPU work at these sizes; the GPU-relevant part is the realized covariance of large panels, which reduces to the batched sufficient statistics of E1. | |
| 42 | +\paragraph{E4 --- realized covariance and connectedness.} For $N \in \{10, 30, 60, 100\}$ symbols with at least 95\,\% regular-bar coverage, daily realized covariance matrices are formed from synchronised 5-minute returns (78 per day) on the days where every symbol has the full grid, and rolling Diebold--Yilmaz total connectedness is computed from a VAR(1) on the daily log-RV panel (120-day windows, $H = 10$). Figure~\ref{fig:e4} and Table~\ref{tab:e4} show the two constraints the charter asked for. Memory is not one of them: the rolling connectedness over 133 windows takes 0.5, 2.3, 8.9 and 29\,ms for $N = 10, 30, 60, 100$ on the CPU and never fails. Conditioning is: the median condition number of the daily realized covariance grows from 81 ($N = 10$) to 422 ($N = 30$) and $5.9\times10^3$ ($N = 60$), and the matrix is singular for $N = 100 > 78$ bins; the Ledoit--Wolf shrinkage intensity rises from 0.16 to 0.33 accordingly. Strict synchronisation also discards 70--164 of the 252 days as $N$ grows. Connectedness itself is not GPU work at these sizes; the GPU-relevant part is the realized covariance of large panels, which reduces to the batched sufficient statistics of E1. | |
| 43 | + | |
| 44 | +\begin{figure}[t] | |
| 45 | +\centering\includegraphics[width=\textwidth]{figures/fig_e4.pdf} | |
| 46 | +\caption{E4 --- realized covariance and connectedness versus the number of assets $N$. Panel A: median and maximum condition number of the daily realized covariance from synchronised 5-minute returns (78 bins; values above $10^{9}$ clipped). Panel B: mean Ledoit--Wolf shrinkage intensity. Panel C: mean total Diebold--Yilmaz connectedness over 133 rolling 120-day windows, annotated with the CPU time for all windows.} | |
| 47 | +\label{fig:e4} | |
| 48 | +\end{figure} | |
| 31 | 49 | |
| 32 | 50 | \begin{table}[t] |
| 33 | 51 | \centering\small |
modified
paper/sections/architecture.tex
+2 −18
@@ -3,23 +3,7 @@ | ||
| 3 | 3 | |
| 4 | 4 | \raid{} is a C++20 library (MIT licence) with four layers (Figure~\ref{fig:arch}): a \emph{core} of strong types and infrastructure, a \emph{data} layer that turns provider bars into validated, calendar-aware return panels, a \emph{statistics/estimators/models} layer that implements the registry, and a \emph{backend} layer that exposes a small set of batched numerical primitives with a CPU reference implementation and a Metal implementation. Python bindings (nanobind) and a command-line tool sit on top. |
| 5 | 5 | |
| 6 | −\begin{figure}[t] | |
| 7 | −\centering | |
| 8 | −\begin{lstlisting}[basicstyle=\ttfamily\tiny] | |
| 9 | − Python (nanobind, GIL released) raid.api / raid.findata CLI: raid devices|validate|estimate|data | |
| 10 | − ------------------------------------------------------------------------------------------------------ | |
| 11 | − core types (Timestamp us UTC, Date, Interval, Backend, Precision) | RaidError | log (redaction) | |
| 12 | − rng (xoshiro256++ streams) | ThreadPool | SHA-256 | Estimate / Provenance / ExecutionPlan / Manifest | |
| 13 | − data MarketCalendar (XNYS rules, DST, early closes) | BarSeries + QualityReport | returns (session-aware) | |
| 14 | − aggregate (calendar-aware bins) | io: CSV fixtures, Parquet/Arrow IPC (projection) | FinData client | |
| 15 | − estimators reductions (Neumaier/pairwise/Welford) | rolling | covariance (HC, HAC, cluster) | OLS/WLS/ridge | |
| 16 | − realized (RV, BPV, RSV, RQ, J, TriPV, RCov) | HAR family | GARCH family (QML) | bootstrap | VAR/GFEVD/DY | |
| 17 | − backend Batch = values + CSR offsets -> sum_squares | sum_abs_adjacent | semivariances | sum_fourth | |
| 18 | − rolling_sum_squares | regression_stats (X'X, X'y, y'y) | garch_nll_batch [CpuBackend FP64 | MetalBackend FP32] | |
| 19 | −\end{lstlisting} | |
| 20 | −\caption{Layers of \raid{}. Estimators consume the backend interface through an \texttt{ExecutionPlan} (\texttt{auto|cpu|metal}); every result is an \texttt{Estimate} carrying coefficients, loss, inference, convergence, warnings, runtime, backend, precision and provenance.} | |
| 21 | −\label{fig:arch} | |
| 22 | −\end{figure} | |
| 6 | +\input{sections/architecture_fig} | |
| 23 | 7 | |
| 24 | 8 | \paragraph{Core types.} Timestamps are microseconds since the Unix epoch in UTC and are never mixed with raw integers; dates use a proleptic-Gregorian day count; intervals are an enumeration. Every estimator returns an \texttt{Estimate}: named coefficients, a loss (residual sum of squares or negative log-likelihood, documented per model), an optional \texttt{InferenceResult} (standard errors, $t$-statistics, $p$-values, covariance, covariance type), a \texttt{Convergence} record (status, iterations, function evaluations, gradient norm, message), warnings, runtime, the backend and precision actually used, and a \texttt{Provenance} record (dataset identifier and version, request identifiers, input hashes, code version). \texttt{Estimate::ok()} requires convergence; a failed optimisation is never reported as a success. An \texttt{ExecutionPlan} selects the backend, the precision policy (FP32 chunk length, compensated accumulation, validation tolerance), determinism, threads and fallback behaviour. |
| 25 | 9 | |
@@ -31,4 +15,4 @@ | ||
| 31 | 15 | |
| 32 | 16 | \paragraph{FinData client.} The data provider is FinData~OS, a private lakehouse serving FirstRate Data bars through a versioned REST API (Section~\ref{sec:setup}). The client authenticates with a bearer key read from the environment only, redacts secrets from logs, applies a token bucket (5 requests per second against a 600 per minute limit), retries with exponential backoff and \texttt{Retry-After}, follows cursor pagination of Arrow IPC pages, chunks requests to the documented caps (31 days of 1-minute bars per call, 12 years per export), polls asynchronous Parquet exports, and keeps a checksum-verified cache keyed by dataset version. Every exchange's request identifier and dataset version are kept for provenance. |
| 33 | 17 | |
| 34 | −\paragraph{Python and CLI.} The nanobind module releases the GIL around every call, accepts NumPy arrays as views and returns NumPy arrays and JSON-like dictionaries. \texttt{raid.findata} exposes the client; \texttt{raid.api} exposes the estimators with a \texttt{backend} argument. The CLI offers \texttt{raid devices}, \texttt{raid validate}, \texttt{raid estimate realized|har|garch|ols} and \texttt{raid data catalog|bars|export|sessions|coverage|query}. | |
| 18 | +\paragraph{Python and CLI.} The nanobind module releases the GIL around every call, accepts NumPy arrays as views and returns NumPy arrays and JSON-like dictionaries. \texttt{raid.findata} exposes the client; \texttt{raid.api} exposes the estimators with a \texttt{backend} argument. The CLI offers \texttt{raid devices}, \texttt{raid validate}, \texttt{raid estimate} (realized, har, garch, ols) and \texttt{raid data} (catalog, bars, export, sessions, coverage, query). | |
added
paper/sections/architecture_fig.tex
+21 −0
@@ -0,0 +1,21 @@ | ||
| 1 | +\begin{figure}[t] | |
| 2 | +\centering | |
| 3 | +\begin{tikzpicture}[ | |
| 4 | + font=\sffamily\scriptsize, node distance=3mm, | |
| 5 | + layer/.style={draw=black!70, line width=0.4pt, rounded corners=1.5pt, fill=black!3, text width=0.93\textwidth, align=left, inner sep=4pt}, | |
| 6 | + tag/.style={font=\sffamily\bfseries\scriptsize, text=black!80}, | |
| 7 | + arr/.style={-{Latex[length=1.6mm]}, line width=0.5pt, black!70}] | |
| 8 | +\node[layer, fill=black!8] (py) {\textbf{Python (nanobind, GIL released)} \quad \texttt{raid.api}: ols $\cdot$ realized\_measures $\cdot$ har\_fit / har\_rolling / har\_batch $\cdot$ garch\_fit $\cdot$ bootstrap $\cdot$ var\_connectedness \quad \texttt{raid.findata}: bars $\cdot$ export $\cdot$ catalog \hfill \textbf{CLI} \texttt{raid devices $|$ validate $|$ estimate $|$ data}}; | |
| 9 | +\node[layer, below=of py] (est) {\textbf{Estimators and models} \quad reductions (Neumaier, pairwise, Welford) $\cdot$ rolling windows $\cdot$ covariance (classical, HC0--HC3, HAC, cluster) $\cdot$ OLS/WLS/ridge $\cdot$ realized measures (RV, BPV, RSV$^\pm$, RQ, $J$, TriPV, RCov) $\cdot$ HAR family $\cdot$ GARCH family (QML) $\cdot$ block bootstrap $\cdot$ VAR, GFEVD, Diebold--Yilmaz, Ledoit--Wolf}; | |
| 10 | +\node[layer, below=of est] (data) {\textbf{Data} \quad MarketCalendar (XNYS rules, DST, early closes) $\cdot$ BarSeries $+$ QualityReport $\cdot$ session-aware returns $\cdot$ calendar-aware aggregation $\cdot$ Parquet / Arrow IPC reader (column projection) $\cdot$ FinData client (bearer key, token bucket, retries, cursors, exports, checksummed cache, provenance)}; | |
| 11 | +\node[layer, below=of data] (be) {\textbf{Backend contract} \quad Batch $=$ values $+$ CSR offsets $\;\longrightarrow\;$ sum\_squares $\cdot$ sum\_abs\_adjacent $\cdot$ semivariances $\cdot$ sum\_fourth $\cdot$ rolling\_sum\_squares $\cdot$ regression\_stats ($X'X$, $X'y$, $y'y$) $\cdot$ garch\_nll\_batch}; | |
| 12 | +\node[draw=black!70, line width=0.4pt, rounded corners=1.5pt, fill=white, anchor=north west, text width=0.445\textwidth, align=left, inner sep=4pt] (cpu) at ($(be.south west)+(0,-3mm)$) {\textbf{CpuBackend} --- FP64, Neumaier per segment, parallel over segments, deterministic order}; | |
| 13 | +\node[draw=black!70, line width=0.4pt, rounded corners=1.5pt, fill=red!4, anchor=north east, text width=0.445\textwidth, align=left, inner sep=4pt] (metal) at ($(be.south east)+(0,-3mm)$) {\textbf{MetalBackend} --- FP32 double-float (TwoSum), 4\,096-element chunks, tree reductions, FP64 combine on the host; prebuilt metallib or runtime MSL; buffer pool}; | |
| 14 | +\node[layer, anchor=north west] (core) at ($(cpu.south west)+(0,-3mm)$) {\textbf{Core} \quad strong types (Timestamp $\mu$s UTC, Date, Interval) $\cdot$ RaidError $\cdot$ structured log with redaction $\cdot$ counter-seeded RNG streams $\cdot$ ThreadPool $\cdot$ SHA-256 $\cdot$ Estimate / InferenceResult / Convergence / Provenance $\cdot$ ExecutionPlan (auto$\,|\,$cpu$\,|\,$metal, PrecisionPolicy) $\cdot$ ExperimentManifest}; | |
| 15 | +\draw[arr] (py.south) -- (est.north); \draw[arr] (est.south) -- (data.north); \draw[arr] (data.south) -- (be.north); | |
| 16 | +\draw[arr] (be.south -| cpu.north) -- (cpu.north); \draw[arr] (be.south -| metal.north) -- (metal.north); | |
| 17 | +\draw[arr, dashed] (metal.west |- cpu.east) -- node[above, font=\sffamily\tiny]{fallback} (cpu.east); | |
| 18 | +\end{tikzpicture} | |
| 19 | +\caption{Layers of \raid{}. Estimators consume the backend contract through an \texttt{ExecutionPlan}; every result is an \texttt{Estimate} that carries coefficients, loss, inference, convergence, warnings, runtime, backend, precision and provenance. The Metal backend falls back to the CPU backend when no device is available or a case is unsupported.} | |
| 20 | +\label{fig:arch} | |
| 21 | +\end{figure} | |
modified
paper/sections/correctness.tex
+7 −1
@@ -32,7 +32,13 @@ Metal vs CPU FP64 & \raid{} CPU backend & RV, BPV, RSV$^\pm$ ($n$ = 390--39\,000 | ||
| 32 | 32 | \end{threeparttable} |
| 33 | 33 | \end{table} |
| 34 | 34 | |
| 35 | −\paragraph{Monte Carlo.} Table~\ref{tab:e5} (experiment E5, $R = 1\,000$ replications) reports bias, RMSE and the empirical coverage of 95\,\% intervals. For GARCH(1,1) with $(\omega,\alpha,\beta) = (0.02, 0.08, 0.90)$ and $n = 2\,000$, biases are at most $0.005$ in absolute value, coverage of the QMLE intervals is 0.93--0.97, and all 1\,000 fits converge; at $n = 500$ the familiar finite-sample downward bias of $\beta$ ($-0.026$) and under-coverage of $\alpha$ (0.90) appear. For HAR with overlapping $h = 5$ targets and HAC inference, coverage of the pseudo-true projection coefficients is 0.90--0.94 ($T = 1\,200$), the usual mild under-coverage of HAC. Realized variance converges at the $\sqrt{2/M}$ rate: the Monte Carlo RMSE is within 10\,\% of $\sqrt{2/M}\,\mathrm{IV}$ for $M \in \{39, 390, 3\,900\}$. The stationary-bootstrap interval for the mean of an AR(1) covers in 89\,\% of replications ($R = 100$, $n = 800$, $b = 10$). | |
| 35 | +\begin{figure}[t] | |
| 36 | +\centering\includegraphics[width=\textwidth]{figures/fig_e5.pdf} | |
| 37 | +\caption{E5 --- Monte Carlo inference ($R = 1\,000$). Panel A: empirical coverage of nominal 95\,\% intervals for the GARCH(1,1) QMLE parameters at $n = 500$ and $n = 2\,000$ and for the four HAR coefficients with HAC inference at $h = 1$ and $h = 5$ ($T = 1\,500$); the dashed line is the nominal level. Panel B: bias (markers) and RMSE (bars) of the GARCH(1,1) estimates.} | |
| 38 | +\label{fig:e5} | |
| 39 | +\end{figure} | |
| 40 | + | |
| 41 | +\paragraph{Monte Carlo.} Figure~\ref{fig:e5} and Table~\ref{tab:e5} (experiment E5, $R = 1\,000$ replications) reports bias, RMSE and the empirical coverage of 95\,\% intervals. For GARCH(1,1) with $(\omega,\alpha,\beta) = (0.02, 0.08, 0.90)$ and $n = 2\,000$, biases are at most $0.005$ in absolute value, coverage of the QMLE intervals is 0.93--0.97, and all 1\,000 fits converge; at $n = 500$ the familiar finite-sample downward bias of $\beta$ ($-0.026$) and under-coverage of $\alpha$ (0.90) appear. For HAR with overlapping $h = 5$ targets and HAC inference, coverage of the pseudo-true projection coefficients is 0.90--0.94 ($T = 1\,200$), the usual mild under-coverage of HAC. Realized variance converges at the $\sqrt{2/M}$ rate: the Monte Carlo RMSE is within 10\,\% of $\sqrt{2/M}\,\mathrm{IV}$ for $M \in \{39, 390, 3\,900\}$. The stationary-bootstrap interval for the mean of an AR(1) covers in 89\,\% of replications ($R = 100$, $n = 800$, $b = 10$). | |
| 36 | 42 | |
| 37 | 43 | \begin{table}[t] |
| 38 | 44 | \centering\small |
modified
paper/sections/endtoend.tex
+8 −2
@@ -12,11 +12,17 @@ | ||
| 12 | 12 | |
| 13 | 13 | \begin{figure}[t] |
| 14 | 14 | \centering\includegraphics[width=.95\textwidth]{figures/fig_e1.pdf} |
| 15 | −\caption{E1 --- throughput (returns per second) of the realized-measures estimator by backend on the real and synthetic panels.} | |
| 15 | +\caption{E1 --- batched realized measures on the 2024 panel. Panel A: throughput of the estimator path by backend for all widths and segment lengths. Panel B: Metal speedup relative to the 18-thread CPU by panel width and segment length. Panel C: gain of the fused kernel over the four-kernel Metal variant.} | |
| 16 | 16 | \label{fig:e1} |
| 17 | 17 | \end{figure} |
| 18 | 18 | |
| 19 | −\paragraph{E6 --- the Arrow/Parquet-to-estimator pipeline.} Table~\ref{tab:e6} compares, on four of the cached export files (2.2--2.5 million bars, 50 symbols, one year each), the \raid{} pipeline---Parquet read with column projection and threaded decoding, session filtering with the exact calendar, within-session log returns, and the five realized measures per session---with the same computation written in pandas, Polars (lazy scan) and DuckDB (SQL window functions). The baselines use an approximate early-close rule (fixed dates) because none of them ships an exchange calendar. \raid{} on the CPU processes a file in 0.118\,s (20 million bars per second), Polars in 0.165\,s and DuckDB in 0.164\,s (both 14 million bars per second), pandas in 2.56\,s. Within \raid{}, decoding the Parquet columns takes 0.065\,s, session filtering and returns 0.021\,s and the estimators 0.031\,s on the CPU; Metal makes the estimator stage slower (0.037\,s) because the stage is below the break-even size of E1 and the conversion is paid again. The answer to RQ3 is therefore moderate: a typed, projected, calendar-exact C++ pipeline is about 1.4$\times$ faster than the best columnar engines and 22$\times$ faster than pandas on this task, and the GPU plays no role in it. Before column projection was added, the same pipeline took 0.35\,s---reading eight of twenty columns instead of all of them, and avoiding the per-row string column, was worth more than any kernel. True zero-copy from Arrow to the estimators was possible for every file (single-chunk, null-free numeric columns), but the session filter necessarily materialises a selection; we count four column copies per file. | |
| 19 | +\paragraph{E6 --- the Arrow/Parquet-to-estimator pipeline.} Figure~\ref{fig:e6} and Table~\ref{tab:e6} compare, on four of the cached export files (2.2--2.5 million bars, 50 symbols, one year each), the \raid{} pipeline---Parquet read with column projection and threaded decoding, session filtering with the exact calendar, within-session log returns, and the five realized measures per session---with the same computation written in pandas, Polars (lazy scan) and DuckDB (SQL window functions). The baselines use an approximate early-close rule (fixed dates) because none of them ships an exchange calendar. \raid{} on the CPU processes a file in 0.118\,s (20 million bars per second), Polars in 0.165\,s and DuckDB in 0.164\,s (both 14 million bars per second), pandas in 2.56\,s. Within \raid{}, decoding the Parquet columns takes 0.065\,s, session filtering and returns 0.021\,s and the estimators 0.031\,s on the CPU; Metal makes the estimator stage slower (0.037\,s) because the stage is below the break-even size of E1 and the conversion is paid again. The answer to RQ3 is therefore moderate: a typed, projected, calendar-exact C++ pipeline is about 1.4$\times$ faster than the best columnar engines and 22$\times$ faster than pandas on this task, and the GPU plays no role in it. Before column projection was added, the same pipeline took 0.35\,s---reading eight of twenty columns instead of all of them, and avoiding the per-row string column, was worth more than any kernel. True zero-copy from Arrow to the estimators was possible for every file (single-chunk, null-free numeric columns), but the session filter necessarily materialises a selection; we count four column copies per file. | |
| 20 | + | |
| 21 | +\begin{figure}[t] | |
| 22 | +\centering\includegraphics[width=\textwidth]{figures/fig_e6.pdf} | |
| 23 | +\caption{E6 --- Parquet-to-estimator pipeline on one export file (2.4 million 1-minute bars, 50 symbols, one year). Panel A: stage breakdown of the \raid{} pipeline on the CPU and on Metal (medians). Panel B: end-to-end seconds per file by engine, medians with 95\,\% bootstrap intervals over 4 files $\times$ 3 repetitions (log scale).} | |
| 24 | +\label{fig:e6} | |
| 25 | +\end{figure} | |
| 20 | 26 | |
| 21 | 27 | \begin{table}[t] |
| 22 | 28 | \centering\small |
modified
paper/sections/estimators.tex
+110 −21
@@ -1,38 +1,127 @@ | ||
| 1 | 1 | \section{Estimator definitions and conventions} |
| 2 | 2 | \label{sec:est} |
| 3 | 3 | |
| 4 | −The method registry (\texttt{docs/methods/registry.md} in the repository) records, for each method, the equation, the inputs, the reference implementation and its options, the test tolerance and the validation status. We reproduce the Tier~A definitions here because conventions---not formulas---are where implementations diverge. | |
| 4 | +The method registry (\texttt{docs/methods/registry.md} in the repository) records, for each method, the equation, the inputs, the reference implementation and its options, the test tolerance and the validation status. We reproduce the Tier~A definitions here because conventions---not formulas---are where implementations diverge. Throughout, $\I{\cdot}$ is the indicator function, $\abs{\cdot}$ the absolute value, and $u$ the unit roundoff of the floating-point format in use. | |
| 5 | 5 | |
| 6 | −\paragraph{Returns and sessions.} With bar closes $p_{t,0},\dots,p_{t,M_t}$ in session $t$ (regular hours 09:30--16:00 ET inclusive of the 16:00 closing print), intraday log returns are $r_{t,j} = \log(p_{t,j}/p_{t,j-1})$, $j = 1,\dots,M_t$; a normal session has $M_t = 390$ and an early close $M_t = 210$ when no minute is missing. Zero-volume minutes are absent at the source and are left absent. | |
| 6 | +\subsection{Returns and sessions} | |
| 7 | 7 | |
| 8 | −\paragraph{Realized measures.} For one session with $M$ returns, | |
| 8 | +Let session $t$ contain the regular-hours bar closes $p_{t,0}, p_{t,1}, \dots, p_{t,M_t}$ (09:30 to 16:00~ET inclusive of the closing print). Intraday log returns are | |
| 9 | +\begin{equation} | |
| 10 | +r_{t,j} \;=\; \log \frac{p_{t,j}}{p_{t,j-1}}, \qquad j = 1, \dots, M_t, | |
| 11 | +\label{eq:returns} | |
| 12 | +\end{equation} | |
| 13 | +so that a complete normal session has $M_t = 390$ and a complete early-close session $M_t = 210$. Zero-volume minutes are absent at the source and are left absent; no return is formed across a session boundary. The overnight return $\log(p_{t,0}/p_{t-1,M_{t-1}})$ and the daily close-to-close return are separate outputs. | |
| 14 | + | |
| 15 | +\subsection{Realized measures} | |
| 16 | + | |
| 17 | +For one session with returns $r_1, \dots, r_M$ (the session index is dropped) and $\mu_p = \E\abs{Z}^p$ for $Z \sim N(0,1)$, so that $\mu_1 = \sqrt{2/\pi}$ and $\mu_{4/3} = 2^{2/3}\,\Gamma(7/6)/\Gamma(1/2)$, | |
| 9 | 18 | \begin{align} |
| 10 | −\RV &= \sum_{j=1}^{M} r_j^2, \qquad | |
| 11 | −\BPV = \mu_1^{-2}\,\frac{M}{M-1}\sum_{j=2}^{M} |r_j|\,|r_{j-1}|, \quad \mu_1 = \sqrt{2/\pi}, \label{eq:rv}\\ | |
| 12 | −\RSV^{-} &= \sum_{j} r_j^2\,\I\{r_j<0\}, \qquad \RSV^{+} = \sum_{j} r_j^2\,\I\{r_j>0\}, \qquad | |
| 13 | −\RQ = \frac{M}{3}\sum_{j=1}^{M} r_j^4, \label{eq:rq}\\ | |
| 14 | −J &= \max(\RV-\BPV, 0), \qquad | |
| 15 | −\mathrm{TriPV} = \mu_{4/3}^{-3}\,\frac{M}{M-2}\sum_{j=3}^{M} |r_j r_{j-1} r_{j-2}|^{4/3}, \label{eq:tripv} | |
| 19 | +\RV &= \sum_{j=1}^{M} r_j^{2}, | |
| 20 | +& | |
| 21 | +\BPV &= \mu_1^{-2}\,\frac{M}{M-1}\sum_{j=2}^{M} \abs{r_j}\,\abs{r_{j-1}}, | |
| 22 | +\label{eq:rv-bpv}\\[2pt] | |
| 23 | +\RSV^{-} &= \sum_{j=1}^{M} r_j^{2}\,\I{r_j < 0}, | |
| 24 | +& | |
| 25 | +\RSV^{+} &= \sum_{j=1}^{M} r_j^{2}\,\I{r_j > 0}, | |
| 26 | +\label{eq:rsv}\\[2pt] | |
| 27 | +\RQ &= \frac{M}{3}\sum_{j=1}^{M} r_j^{4}, | |
| 28 | +& | |
| 29 | +J &= \max\bigl(\RV - \BPV,\, 0\bigr), | |
| 30 | +\label{eq:rq-j}\\[2pt] | |
| 31 | +\TriPV &= \mu_{4/3}^{-3}\,\frac{M}{M-2}\sum_{j=3}^{M} \abs{r_j\, r_{j-1}\, r_{j-2}}^{4/3}. | |
| 32 | +\label{eq:tripv} | |
| 16 | 33 | \end{align} |
| 17 | −with $\mu_p = \E|Z|^p$ for $Z\sim N(0,1)$. The finite-sample factors are options. Two conventions of \texttt{highfrequency} differ: \texttt{rBPCov} with \texttt{makeReturns = FALSE} applies no $M/(M-1)$ factor, and \texttt{rQuar} uses $(M+1)/3$, the number of prices, instead of $M/3$; the reference tests compare after conversion. Realized covariance on synchronised returns is $\sum_j r_j r_j'$ without demeaning. | |
| 34 | +The finite-sample factors $M/(M-1)$ and $M/(M-2)$ are options. Two conventions of \texttt{highfrequency} differ from \eqref{eq:rv-bpv}--\eqref{eq:rq-j}: \texttt{rBPCov} with \texttt{makeReturns = FALSE} applies no $M/(M-1)$ factor, and \texttt{rQuar} uses $(M+1)/3$, the number of prices, instead of $M/3$; the reference tests compare after conversion. For $N$ assets with synchronised return vectors $\bm{r}_j \in \mathbb{R}^N$, the realized covariance is $\RCov = \sum_{j=1}^{M} \bm{r}_j \bm{r}_j'$ without demeaning, and the realized correlation is $\diag(\RCov)^{-1/2}\,\RCov\,\diag(\RCov)^{-1/2}$. | |
| 18 | 35 | |
| 19 | −\paragraph{Regression and inference.} OLS, WLS and ridge are solved by column-pivoting Householder QR (SVD when rank-deficient; ridge by augmentation). Covariances follow \texttt{statsmodels}: classical $s^2(X'X)^{-1}$; HC0--HC3 with the leverage corrections; HAC with Bartlett weights and $L$ lags (default $\lfloor 4(n/100)^{2/9}\rfloor$, no small-sample correction); one-way cluster with $\frac{G}{G-1}\frac{n-1}{n-k}$. $p$-values use $t(n-k)$. Rolling OLS refits each window. | |
| 36 | +\subsection{Regression and inference} | |
| 20 | 37 | |
| 21 | −\paragraph{HAR family.} With daily $\RV_t$, lags $(1,5,22)$ and horizon $h$, | |
| 38 | +For $\by \in \mathbb{R}^n$ and $\bX \in \mathbb{R}^{n \times k}$, OLS, WLS and ridge solve | |
| 22 | 39 | \begin{equation} |
| 23 | −y_{t+h} = \beta_0 + \beta_d \RV_t + \beta_w \overline{\RV}^{(5)}_t + \beta_m \overline{\RV}^{(22)}_t \;[+\,\beta_q(\sqrt{\RQ_t}-\sqrt{\overline{\RQ}})\RV_t] + \varepsilon_{t+h}, \qquad \overline{\RV}^{(k)}_t = \tfrac1k\sum_{j=0}^{k-1}\RV_{t-j}, | |
| 40 | +\hat{\bbeta} \;=\; \argmin_{\bbeta}\; (\by - \bX\bbeta)'\,\bm{W}\,(\by - \bX\bbeta) \;+\; \lambda\,\bbeta'\bm{P}\bbeta, | |
| 41 | +\label{eq:ols} | |
| 42 | +\end{equation} | |
| 43 | +with $\bm{W}$ diagonal ($\bm{W} = \bm{I}$ for OLS) and $\bm{P}$ the identity on the penalised coefficients, by column-pivoting Householder QR of the (augmented) weighted design; the SVD minimum-norm solution is used when the design is rank deficient, and the reciprocal condition number is reported. With residuals $\be$ and leverages $h_{ii}$, the covariance estimators follow \texttt{statsmodels}: | |
| 44 | +\begin{align} | |
| 45 | +\widehat{\Var}_{\mathrm{cl}}(\hat\bbeta) &= s^2 (\bX'\bX)^{-1}, \qquad s^2 = \frac{\be'\be}{n-k}, | |
| 46 | +\label{eq:cov-classical}\\ | |
| 47 | +\widehat{\Var}_{\mathrm{HC}_m}(\hat\bbeta) &= (\bX'\bX)^{-1}\Bigl(\sum_{i=1}^{n} \omega_{m,i}\, e_i^{2}\, \bm{x}_i\bm{x}_i'\Bigr)(\bX'\bX)^{-1}, | |
| 48 | +\label{eq:cov-hc}\\ | |
| 49 | +&\qquad \omega_{0,i} = 1,\qquad \omega_{1,i} = \frac{n}{n-k},\qquad \omega_{2,i} = \frac{1}{1-h_{ii}},\qquad \omega_{3,i} = \frac{1}{(1-h_{ii})^{2}}, \nonumber\\ | |
| 50 | +\widehat{\Var}_{\mathrm{HAC}}(\hat\bbeta) &= (\bX'\bX)^{-1}\Bigl(\bm{\Gamma}_0 + \sum_{\ell=1}^{L}\Bigl(1 - \frac{\ell}{L+1}\Bigr)\bigl(\bm{\Gamma}_\ell + \bm{\Gamma}_\ell'\bigr)\Bigr)(\bX'\bX)^{-1}, | |
| 51 | +\quad \bm{\Gamma}_\ell = \sum_{i=\ell+1}^{n} e_i e_{i-\ell}\, \bm{x}_i \bm{x}_{i-\ell}', | |
| 52 | +\label{eq:cov-hac}\\ | |
| 53 | +\widehat{\Var}_{\mathrm{CR}}(\hat\bbeta) &= \frac{G}{G-1}\,\frac{n-1}{n-k}\,(\bX'\bX)^{-1}\Bigl(\sum_{g=1}^{G} \bm{s}_g \bm{s}_g'\Bigr)(\bX'\bX)^{-1}, | |
| 54 | +\quad \bm{s}_g = \sum_{i \in g} e_i\, \bm{x}_i, | |
| 55 | +\label{eq:cov-cluster} | |
| 56 | +\end{align} | |
| 57 | +with the Bartlett lag $L = \lfloor 4(n/100)^{2/9} \rfloor$ by default and no small-sample correction for HAC. Two-sided $p$-values use the $t_{n-k}$ distribution. Rolling OLS refits each window by QR; batched OLS accumulates $(\bX'\bX, \bX'\by, \by'\by)$ per series on the backend and solves in FP64. | |
| 58 | + | |
| 59 | +\subsection{HAR family} | |
| 60 | + | |
| 61 | +With daily realized variances $\RV_t$, averaging windows $\overline{\RV}^{(k)}_t = k^{-1}\sum_{j=0}^{k-1}\RV_{t-j}$ for $k \in \{1, 5, 22\}$, and horizon $h \ge 1$, the HAR-RV regression is | |
| 62 | +\begin{equation} | |
| 63 | +y_{t+h} \;=\; \beta_0 \;+\; \beta_d\,\RV_t \;+\; \beta_w\,\overline{\RV}^{(5)}_t \;+\; \beta_m\,\overline{\RV}^{(22)}_t \;+\; \varepsilon_{t+h}, | |
| 64 | +\qquad | |
| 65 | +y_{t+h} = \frac{1}{h}\sum_{j=1}^{h} \RV_{t+j}, | |
| 24 | 66 | \label{eq:har} |
| 25 | 67 | \end{equation} |
| 26 | −where $y_{t+h}$ is the average of $\RV_{t+1},\dots,\RV_{t+h}$ (the single-day target is an option) and the bracketed HARQ term follows \citet{BollerslevPattonQuaedvlieg2016} in the demeaned form used by \texttt{highfrequency}; log-HAR applies the same regression to $\log\RV$ and back-transforms with $\exp$ (the smearing factor is reported, not applied). Overlapping targets for $h>1$ induce MA$(h-1)$ errors, so the HAC lag is at least $h-1$. Rolling estimation is strictly causal: the model fitted on observations up to $t$ forecasts $y_{t+h}$. Forecast losses are MSE and QLIKE $= y/\hat y - \log(y/\hat y) - 1$ \citep{Patton2011}. | |
| 68 | +where the single-day target $y_{t+h} = \RV_{t+h}$ is an option. HARQ \citep{BollerslevPattonQuaedvlieg2016} adds the term | |
| 69 | +\begin{equation} | |
| 70 | +\beta_q\,\bigl(\sqrt{\RQ_t} - \sqrt{\overline{\RQ}}\,\bigr)\,\RV_t, | |
| 71 | +\qquad \overline{\RQ} = \frac{1}{T}\sum_{t=1}^{T}\RQ_t \ \text{over the estimation sample}, | |
| 72 | +\label{eq:harq} | |
| 73 | +\end{equation} | |
| 74 | +in the demeaned form used by \texttt{highfrequency}; log-HAR applies \eqref{eq:har} to $\log\RV$ and back-transforms with $\exp$ (the smearing factor $T^{-1}\sum_t e^{\hat\varepsilon_t}$ is reported, not applied). Overlapping targets for $h > 1$ induce MA$(h-1)$ errors, so the HAC lag is at least $h-1$. Rolling estimation is strictly causal: the model fitted on observations up to $t$ forecasts $y_{t+h}$. Forecast losses are the mean squared error and | |
| 75 | +\begin{equation} | |
| 76 | +\mathrm{QLIKE}(y, \hat y) \;=\; \frac{y}{\hat y} \;-\; \log\frac{y}{\hat y} \;-\; 1, | |
| 77 | +\label{eq:qlike} | |
| 78 | +\end{equation} | |
| 79 | +the robust loss of \citet{Patton2011}. | |
| 80 | + | |
| 81 | +\subsection{GARCH family} | |
| 27 | 82 | |
| 28 | −\paragraph{GARCH family.} With $\varepsilon_t = r_t - \mu$, | |
| 83 | +With $\varepsilon_t = r_t - \mu$ and $z_t = \varepsilon_t/\sigma_t$, the conditional variance recursions are | |
| 29 | 84 | \begin{align} |
| 30 | −\text{GARCH}(1,1):\;& \sigma_t^2 = \omega + \alpha\varepsilon_{t-1}^2 + \beta\sigma_{t-1}^2, \qquad | |
| 31 | −\text{GJR}:\; \sigma_t^2 = \omega + (\alpha + \gamma\I\{\varepsilon_{t-1}<0\})\varepsilon_{t-1}^2 + \beta\sigma_{t-1}^2, \label{eq:garch}\\ | |
| 32 | −\text{EGARCH}:\;& \log\sigma_t^2 = \omega + \alpha(|z_{t-1}| - \E|z|) + \gamma z_{t-1} + \beta\log\sigma_{t-1}^2, \quad z = \varepsilon/\sigma. \label{eq:egarch} | |
| 85 | +\text{GARCH}(1,1):\qquad & \sigma_t^{2} = \omega + \alpha\,\varepsilon_{t-1}^{2} + \beta\,\sigma_{t-1}^{2}, | |
| 86 | +\label{eq:garch}\\ | |
| 87 | +\text{GJR}(1,1):\qquad & \sigma_t^{2} = \omega + \bigl(\alpha + \gamma\,\I{\varepsilon_{t-1} < 0}\bigr)\varepsilon_{t-1}^{2} + \beta\,\sigma_{t-1}^{2}, | |
| 88 | +\label{eq:gjr}\\ | |
| 89 | +\text{EGARCH}(1,1):\qquad & \log\sigma_t^{2} = \omega + \alpha\bigl(\abs{z_{t-1}} - \E\abs{z}\bigr) + \gamma\, z_{t-1} + \beta\,\log\sigma_{t-1}^{2}, | |
| 90 | +\label{eq:egarch} | |
| 33 | 91 | \end{align} |
| 34 | −Errors are Gaussian or standardised Student-$t$ with $\nu\in(2.05, 500)$ (the bounds of \texttt{arch}); $\sigma_0^2$ is the sample variance of $\varepsilon$ (as \texttt{rugarch}; \texttt{arch} uses a backcast, which explains the $5\times10^{-4}$ log-likelihood tolerance against it). Estimation is quasi-maximum likelihood by BFGS on transformed parameters (log for $\omega$, logistic for persistence and the $\alpha/\beta$ share, $\tanh$ for the EGARCH $\beta$), with central-difference gradients validated across step sizes; up to two alternative starting points are tried on non-convergence, and the number tried is reported. Standard errors are the Bollerslev--Wooldridge sandwich $H^{-1}\,\mathrm{OPG}\,H^{-1}$ by default (Hessian and OPG are options). Forecasts are analytical for GARCH/GJR and by seeded simulation beyond one step for EGARCH. The \texttt{rugarch} EGARCH maps its \texttt{alpha1} to our $\gamma$ and its \texttt{gamma1} to our $\alpha$. | |
| 92 | +with constraints $\omega > 0$, $\alpha, \beta \ge 0$, $\alpha + \beta < 1$ (GARCH), $\alpha + \gamma \ge 0$, $\alpha + \beta + \gamma/2 < 1$ (GJR) and $\abs{\beta} < 1$ (EGARCH). Errors are Gaussian or standardised Student-$t$ with $\nu \in (2.05, 500)$, the bounds of \texttt{arch}, and the quasi log-likelihood is | |
| 93 | +\begin{equation} | |
| 94 | +\ell_T(\btheta) \;=\; \sum_{t=1}^{T}\Bigl[\, c(\nu) \;-\; \tfrac12 \log\sigma_t^{2} \;-\; \tfrac{\nu+1}{2}\,\log\Bigl(1 + \frac{\varepsilon_t^{2}}{(\nu-2)\,\sigma_t^{2}}\Bigr)\Bigr], | |
| 95 | +\qquad | |
| 96 | +c(\nu) = \log\frac{\Gamma\bigl(\tfrac{\nu+1}{2}\bigr)}{\Gamma\bigl(\tfrac{\nu}{2}\bigr)\sqrt{\pi(\nu-2)}}, | |
| 97 | +\label{eq:garch-ll} | |
| 98 | +\end{equation} | |
| 99 | +with the Gaussian case $\ell_T(\btheta) = \sum_t\bigl[-\tfrac12\log 2\pi - \tfrac12\log\sigma_t^2 - \varepsilon_t^2/(2\sigma_t^2)\bigr]$. The recursion starts at $\sigma_0^2 = \widehat{\Var}(\varepsilon)$, the sample variance, as in \texttt{rugarch}; \texttt{arch} uses a backcast, which explains the $5\times10^{-4}$ log-likelihood tolerance against it. Estimation maximises \eqref{eq:garch-ll} by BFGS on transformed parameters (Appendix~B), with central-difference gradients validated across step sizes; up to two alternative starting points are tried on non-convergence, and the number tried is reported. Standard errors are the Bollerslev--Wooldridge sandwich | |
| 100 | +\begin{equation} | |
| 101 | +\widehat{\Var}(\hat\btheta) \;=\; \hat{\bm{H}}^{-1}\Bigl(\sum_{t=1}^{T} \hat{\bm{s}}_t \hat{\bm{s}}_t'\Bigr)\hat{\bm{H}}^{-1}, | |
| 102 | +\qquad \hat{\bm{H}} = -\frac{\partial^2 \ell_T}{\partial\btheta\,\partial\btheta'}\Big|_{\hat\btheta},\quad \hat{\bm{s}}_t = \frac{\partial \ell_t}{\partial\btheta}\Big|_{\hat\btheta}, | |
| 103 | +\label{eq:qmle-cov} | |
| 104 | +\end{equation} | |
| 105 | +by default (Hessian and OPG are options). Forecasts are analytical for GARCH and GJR, | |
| 106 | +\begin{equation} | |
| 107 | +\hat\sigma^2_{T+1} = \omega + \bigl(\alpha + \gamma\,\I{\varepsilon_T < 0}\bigr)\varepsilon_T^{2} + \beta\,\sigma_T^{2}, | |
| 108 | +\qquad | |
| 109 | +\hat\sigma^2_{T+h} = \omega + \bigl(\alpha + \tfrac{\gamma}{2} + \beta\bigr)\hat\sigma^2_{T+h-1},\quad h \ge 2, | |
| 110 | +\label{eq:garch-fc} | |
| 111 | +\end{equation} | |
| 112 | +and by seeded simulation beyond one step for EGARCH. The \texttt{rugarch} EGARCH parameterisation maps its \texttt{alpha1} to our $\gamma$ and its \texttt{gamma1} to our $\alpha$. | |
| 113 | + | |
| 114 | +\subsection{Bootstrap} | |
| 35 | 115 | |
| 36 | −\paragraph{Bootstrap.} The stationary bootstrap of \citet{PolitisRomano1994} draws blocks of geometric length (mean $b$, default $n^{1/3}$) starting at uniform indices with circular wrap; the circular block bootstrap uses fixed length; an iid option exists for comparison only. Replica $r$ uses the stream $(\text{seed}, r)$ of a counter-seeded xoshiro256++ generator, so replicas are independent of the thread that computes them. | |
| 116 | +The stationary bootstrap of \citet{PolitisRomano1994} draws blocks of geometric length with mean $b$ (default $b = n^{1/3}$) starting at uniform indices with circular wrap-around; the circular block bootstrap uses fixed length $b$; an iid option exists for comparison only. Replica $r$ uses the stream $(\text{seed}, r)$ of a counter-seeded xoshiro256++ generator, so replicas are independent of the thread that computes them and of the order of computation. | |
| 37 | 117 | |
| 38 | −\paragraph{Tier B.} VAR$(p)$ with intercept by multivariate least squares, generalized FEVD \citep{KoopPesaranPotter1996,PesaranShin1998} normalised row-wise, Diebold--Yilmaz total, from, to and net connectedness \citep{DieboldYilmaz2012,DieboldYilmaz2014} with rolling windows, and Ledoit--Wolf shrinkage toward the scaled identity \citep{LedoitWolf2004}. Realized GARCH, HEAVY, DCC, GARCH-MIDAS, state-space and change-point methods, pre-averaging, realized kernels and order-flow methods are \emph{not} implemented; the registry says so. | |
| 118 | +\subsection{Tier B} | |
| 119 | + | |
| 120 | +For a panel $\by_t \in \mathbb{R}^N$, VAR$(p)$ with intercept $\by_t = \bm{c} + \sum_{i=1}^{p}\bA_i\by_{t-i} + \bu_t$ is estimated by multivariate least squares, with $\bSigma_u = (T-Np-1)^{-1}\sum_t \hat\bu_t\hat\bu_t'$ and moving-average matrices $\bPhi_0 = \bm{I}$, $\bPhi_h = \sum_{i=1}^{\min(h,p)}\bA_i\bPhi_{h-i}$. The generalized forecast-error variance decomposition at horizon $H$ \citep{KoopPesaranPotter1996,PesaranShin1998} and its row normalisation are | |
| 121 | +\begin{equation} | |
| 122 | +\theta_{ij}(H) = \frac{\sigma_{jj}^{-1}\sum_{h=0}^{H-1}\bigl(\be_i'\bPhi_h\bSigma_u\be_j\bigr)^{2}}{\sum_{h=0}^{H-1}\be_i'\bPhi_h\bSigma_u\bPhi_h'\be_i}, | |
| 123 | +\qquad | |
| 124 | +\tilde\theta_{ij}(H) = \frac{\theta_{ij}(H)}{\sum_{j=1}^{N}\theta_{ij}(H)}, | |
| 125 | +\label{eq:gfevd} | |
| 126 | +\end{equation} | |
| 127 | +and the Diebold--Yilmaz measures \citep{DieboldYilmaz2012,DieboldYilmaz2014} are the total connectedness $C(H) = 100\,N^{-1}\sum_{i \ne j}\tilde\theta_{ij}(H)$, the from-others $C_{i\leftarrow\bullet} = 100\sum_{j \ne i}\tilde\theta_{ij}$, the to-others $C_{\bullet\leftarrow j} = 100\sum_{i \ne j}\tilde\theta_{ij}$ and the net $C_{\bullet\leftarrow j} - C_{j\leftarrow\bullet}$, computed on rolling windows. Ledoit--Wolf shrinkage \citep{LedoitWolf2004} is $\bSigma^{\star} = (1-\delta^{\star})\bm{S} + \delta^{\star}\,\mu\,\bm{I}$ with $\mu = \tr(\bm{S})/N$ and the analytical $\delta^{\star}$. Realized GARCH, HEAVY, DCC, GARCH-MIDAS, state-space and change-point methods, pre-averaging, realized kernels and order-flow methods are \emph{not} implemented; the registry says so. | |
modified
paper/sections/metal.tex
+22 −5
@@ -1,15 +1,32 @@ | ||
| 1 | 1 | \section{Metal execution and precision strategy} |
| 2 | 2 | \label{sec:metal} |
| 3 | 3 | |
| 4 | −\paragraph{The constraint.} Apple GPUs and the Metal Shading Language have no 64-bit floating-point type; a kernel declaring a \texttt{double} does not compile (``\texttt{'double' is not supported in Metal}''). Econometric accumulators---sums of $4\times10^5$ squared returns of magnitude $10^{-6}$, Gram matrices $X'X$ with regressors of very different scales, log-likelihoods summed over thousands of observations---lose several digits in naive FP32. The strategy has four parts. | |
| 4 | +\paragraph{The constraint.} Apple GPUs and the Metal Shading Language have no 64-bit floating-point type; a kernel declaring a \texttt{double} does not compile (``\texttt{'double' is not supported in Metal}''). Econometric accumulators---sums of $4\times10^{5}$ squared returns of magnitude $10^{-6}$, Gram matrices $\bX'\bX$ with regressors of very different scales, log-likelihoods summed over thousands of observations---lose several digits in naive FP32, whose unit roundoff is $u = 2^{-24} \approx 6\times10^{-8}$. The strategy has four parts. | |
| 5 | 5 | |
| 6 | −\paragraph{1. Double-float accumulation.} Every accumulator on the GPU is a pair $(h, \ell)$ of \texttt{float}s representing $h+\ell$ with $|\ell|\le\frac12\,\mathrm{ulp}(h)$. Additions use the error-free transformation TwoSum \citep{Dekker1971,Neumaier1974}: $s = h + v$, $e = (h - (s - (s - h))) + (v - (s - h))$, followed by a renormalisation. The pair carries about 48 significant bits; the error analysis of \citet{Higham1993} gives a bound of order $n\,u^2$ rather than $n\,u$ for a sum of $n$ terms with unit roundoff $u$. The Metal compiler must not reassociate: the kernels are compiled with \texttt{-fno-fast-math} offline and with \texttt{fastMathEnabled = NO} at runtime. Per-element inputs are still rounded to FP32 on the way in; for a sum of squares of $M$ independent returns this input rounding contributes a relative error of order $2u\sqrt{3/M}$, i.e. about $10^{-8}$ for $M = 390$, which is what we observe (Section~\ref{sec:correct}). | |
| 6 | +\paragraph{1. Double-float accumulation.} Every accumulator on the GPU is a pair $(h, \ell)$ of \texttt{float}s representing $h + \ell$ with $\abs{\ell} \le \tfrac12\ulp(h)$. Additions use the error-free transformation TwoSum \citep{Dekker1971,Neumaier1974}: for floats $a, b$, | |
| 7 | +\begin{equation} | |
| 8 | +s = \fl(a + b), \qquad e = \bigl(a - (s - b')\bigr) + \bigl(b - b'\bigr),\quad b' = \fl(s - a), \qquad a + b = s + e \ \text{exactly}, | |
| 9 | +\label{eq:twosum} | |
| 10 | +\end{equation} | |
| 11 | +followed by a renormalisation of the pair. The pair carries about 48 significant bits, and the analysis of \citet{Higham1993} bounds the relative error of a compensated sum of $n$ terms by $2u + O(nu^{2})$ instead of $(n-1)u$ for recursive summation. The Metal compiler must not reassociate: the kernels are compiled with \texttt{-fno-fast-math} offline and \texttt{fastMathEnabled = NO} at runtime. The inputs are still rounded to FP32 on the way in; for a sum of squares of $M$ independent returns with rounding errors $\delta_j$, $\abs{\delta_j} \le u$, | |
| 12 | +\begin{equation} | |
| 13 | +\frac{\fl(\RV) - \RV}{\RV} \;\approx\; \frac{2\sum_j \delta_j\, r_j^{2}}{\sum_j r_j^{2}}, \qquad | |
| 14 | +\sqrt{\E\Bigl[\Bigl(\tfrac{2\sum_j \delta_j r_j^{2}}{\sum_j r_j^{2}}\Bigr)^{2}\Bigr]} \;\lesssim\; 2u\,\sqrt{\frac{\sum_j r_j^{4}}{\bigl(\sum_j r_j^{2}\bigr)^{2}}} \;\approx\; 2u\sqrt{\frac{3}{M}} | |
| 15 | +\label{eq:rounding} | |
| 16 | +\end{equation} | |
| 17 | +for Gaussian returns, i.e. about $10^{-8}$ for $M = 390$ and decreasing in $M$. Figure~\ref{fig:precision} shows that the maximum error observed across the segments of each E1 case follows this $M^{-1/2}$ law and lies within a factor of about two of the expected maximum of 926 Gaussian deviations ($\approx 3.2$ root-mean-square units), the excess at short segments reflecting the heavy tails of real returns. | |
| 7 | 18 | |
| 8 | −\paragraph{2. Chunked partial sums combined in FP64.} A segment longer than the chunk length (4\,096 elements by default, a \texttt{PrecisionPolicy} parameter) is split into chunks; one threadgroup of 256 threads reduces each chunk with a fixed-order tree reduction of double-float values and writes a \texttt{float2}; the host adds the chunks in FP64 in a fixed order and adds the chunk-boundary cross terms ($|r_j||r_{j-1}|$ across a boundary) in FP64. The combination order is deterministic, so results are bitwise reproducible across runs for a given chunk length and threadgroup width (verified in Section~\ref{sec:correct}). | |
| 19 | +\begin{figure}[t] | |
| 20 | +\centering\includegraphics[width=.62\textwidth]{figures/fig_precision.pdf} | |
| 21 | +\caption{Precision of the Metal realized variance. Maximum relative difference between the Metal (FP32 double-float) and the CPU FP64 realized variance across the segments of each E1 case, versus the number of returns per segment, with the root-mean-square prediction of \eqref{eq:rounding} and its expected maximum over 926 segments.} | |
| 22 | +\label{fig:precision} | |
| 23 | +\end{figure} | |
| 9 | 24 | |
| 10 | −\paragraph{3. Mixed precision by design.} Quantities that are cheap and precision-critical stay in FP64 on the CPU: the finite-sample constants, the initial variance of a GARCH recursion, the Student-$t$ density constants (Metal has no \texttt{lgamma}), the regression solve from the GPU-accumulated sufficient statistics, and all inference. The batched GARCH kernel scales each series to unit variance on the host (FP64) and un-scales the likelihood by $n\log s$, which conditions the FP32 recursion for returns of magnitude $10^{-2}$. | |
| 25 | +\paragraph{2. Chunked partial sums combined in FP64.} A segment longer than the chunk length ($4\,096$ elements by default, a \texttt{PrecisionPolicy} parameter) is split into chunks; one threadgroup of 256 threads reduces each chunk with a fixed-order tree reduction of double-float values and writes a \texttt{float2}; the host adds the chunks in FP64 in a fixed order and adds the chunk-boundary cross terms ($\abs{r_j}\abs{r_{j-1}}$ across a boundary) in FP64. The combination order is deterministic, so results are bitwise reproducible across runs for a given chunk length and threadgroup width (Section~\ref{sec:correct}). | |
| 11 | 26 | |
| 12 | −\paragraph{4. Kernel fusion.} The realized kernel computes the five accumulators of \eqref{eq:rv}--\eqref{eq:rq} in a single pass over the data; the four-kernel variant exists only for the ablation. Regression statistics accumulate the upper triangle of $X'X$, $X'y$ and $y'y$ ($k\le 8$) in registers. The GARCH kernel runs one recursion per thread for one (series, parameter vector) pair---the parallelism comes from the number of parameter vectors, not from the recursion. | |
| 27 | +\paragraph{3. Mixed precision by design.} Quantities that are cheap and precision-critical stay in FP64 on the CPU: the finite-sample constants, the initial variance of a GARCH recursion, the Student-$t$ constant $c(\nu)$ and $\E\abs{z}$ in \eqref{eq:garch-ll}--\eqref{eq:egarch} (Metal has no \texttt{lgamma}), the regression solve from the GPU-accumulated sufficient statistics, and all inference. The batched GARCH kernel scales each series to unit variance on the host, $r_t' = r_t/s$, $\mu' = \mu/s$, $\omega' = \omega/s^{2}$ (EGARCH: $\omega' = \omega + (1-\beta)\log s^{-2}$), and the host un-scales the log-likelihood by $-T\log s$, which conditions the FP32 recursion for returns of magnitude $10^{-2}$. | |
| 28 | + | |
| 29 | +\paragraph{4. Kernel fusion.} The realized kernel computes the five accumulators of \eqref{eq:rv-bpv}--\eqref{eq:rq-j} in a single pass over the data; the four-kernel variant exists only for the ablation. Regression statistics accumulate the upper triangle of $\bX'\bX$, $\bX'\by$ and $\by'\by$ ($k \le 8$) in registers. The GARCH kernel runs one recursion per thread for one (series, parameter vector) pair---the parallelism comes from the number of parameter vectors, not from the recursion. | |
| 13 | 30 | |
| 14 | 31 | \paragraph{Host side.} The backend owns one device, one command queue, one library and one pipeline cache per process; buffers come from a pool keyed by power-of-two size class in shared storage (unified memory). The shader library is loaded from a prebuilt \texttt{.metallib} when present and compiled from the embedded MSL source otherwise; both paths are measured. GPU time is read from the command buffer's \texttt{GPUStartTime}/\texttt{GPUEndTime}; wall-clock includes FP64$\to$FP32 conversion into the shared buffer, encoding, commit and completion. Metal does not expose the GPU core count; the manifest records it from \texttt{system\_profiler}. |
| 15 | 32 | |
modified
paper/sections/microbenchmarks.tex
+3 −3
@@ -11,8 +11,8 @@ | ||
| 11 | 11 | \end{table} |
| 12 | 12 | |
| 13 | 13 | \begin{figure}[t] |
| 14 | −\centering\includegraphics[width=\textwidth]{figures/fig_primitives_throughput.pdf} | |
| 15 | −\caption{Throughput of the primitives (elements per second, medians) versus problem size on the laptop.} | |
| 14 | +\centering\includegraphics[width=\textwidth]{figures/fig_primitives.pdf} | |
| 15 | +\caption{Primitive kernels on the laptop (M5 Max). Panels A--D: throughput (elements per second, medians of 5 repetitions) of the single-thread CPU, the 18-thread CPU and the Metal backend versus problem size, for every $(S, n)$ shape. Panels E--H: Metal speedup relative to the 18-thread CPU with 95\,\% bootstrap intervals of the median; the dashed line is break-even.} | |
| 16 | 16 | \label{fig:prim} |
| 17 | 17 | \end{figure} |
| 18 | 18 | |
@@ -31,6 +31,6 @@ Four findings answer RQ1 at the kernel level. First, the fused realized kernel i | ||
| 31 | 31 | |
| 32 | 32 | \begin{figure}[t] |
| 33 | 33 | \centering\includegraphics[width=.8\textwidth]{figures/fig_e7.pdf} |
| 34 | −\caption{E7 --- Metal/CPU-MT speedup versus GPU core count for three kernels on eleven Apple Silicon machines.} | |
| 34 | +\caption{E7 --- generation sweep. Panel A: Metal speedup relative to the CPU with all threads, at $1\,000\times3\,900$ (realized, regression, rolling) and $1\,000\times390$, $P = 64$ (GARCH), versus GPU core count on eleven machines (3 repetitions). Panel B: runtime compilation of the embedded Metal source, cold (salted source, system cache defeated) and warm.} | |
| 35 | 35 | \label{fig:e7} |
| 36 | 36 | \end{figure} |
modified
paper/sections/reproducibility.tex
+1 −1
@@ -1,7 +1,7 @@ | ||
| 1 | 1 | \section{Reproducibility and licensing} |
| 2 | 2 | \label{sec:repro} |
| 3 | 3 | |
| 4 | −The repository (\texttt{wp40\_raid}, git remote on the author's SPB Git server) contains the library, kernels, bindings, CLI, tests, benchmark drivers, every raw benchmark record (\texttt{experiments/results/raw/*.jsonl}) with its manifest (host, OS, compiler, library versions, load, thermal state, seeds, input hashes), the scripts that generate every table and figure of this paper from those records (\texttt{experiments/scripts/make\_paper\_assets.py}; the SHA-256 of each input and output is written to \texttt{experiments/manifests/paper\_assets.json}), and \texttt{RESEARCH\_STATUS.md}, which lists what is implemented, verified, experimental and deliberately not implemented. No number in the text was typed by hand: the prose quotes \texttt{experiments/results/key\_numbers.json}. | |
| 4 | +The repository (\path{wp40_raid}, git remote on the author's SPB Git server) contains the library, kernels, bindings, CLI, tests, benchmark drivers, every raw benchmark record (\path{experiments/results/raw/*.jsonl}) with its manifest (host, OS, compiler, library versions, load, thermal state, seeds, input hashes), the scripts that generate every table and figure of this paper from those records (\path{experiments/scripts/make_paper_assets.py}; the SHA-256 of each input and output is written to \path{experiments/manifests/paper_assets.json}), and \path{RESEARCH_STATUS.md}, which lists what is implemented, verified, experimental and deliberately not implemented. No number in the text was typed by hand: the prose quotes \path{experiments/results/key_numbers.json}. | |
| 5 | 5 | |
| 6 | 6 | A clean checkout builds with CMake presets (\texttt{macos-release} with Metal and Arrow; \texttt{cpu-only} for Linux) and passes 55 test cases; the Python package builds with scikit-build-core. The synthetic fixtures and the reference values produced by \texttt{statsmodels}, \texttt{arch}, \texttt{exchange\_calendars}, \texttt{highfrequency} and \texttt{rugarch} are committed, so the validation suite runs without network access or vendor data. Reproducing E1--E6 requires an API key for the FinData~OS lakehouse; the universe selection is seeded and the dataset version is pinned, so the same bars are served to a reader with access. E7 requires only the bundle and SSH access to the machines. |
| 7 | 7 | |
| 8 | 8 | |