options: OPRA daily bars, a quote-quality check and a 2024-02 backfill

The indicative feed's quotes are derived, not OPRA quotes (Alpaca docs);
the OPRA BBO needs the $99/month plan. Option bars are computed from OPRA
prints and are free beyond the last 15 minutes, back to 2024-02.

- snapshot also stores OPRA daily bars for every contract expiring within
  0-100 days (active and same-day inactive, so 0DTE volume is kept), and
  bases option volume, os_volume and pc_volume on them.
- snapshot samples three near-the-money contracts per name and compares
  the last 15:45-16:00 OPRA print with the indicative quote (share inside,
  IV error). status gives the pre-set verdict after 20 sessions: adequate
  if >= 90% inside and median |IV error| <= 1 vol point.
- backfill fetches daily bars since 2024-02-01 for every S&P 500 member
  since then, one resumable file per underlying, pausing 16:10-16:55 New
  York; volume-panel builds the daily option-volume panel.
- Bar requests send a timestamp end capped at now - 16 min: a date-only
  end of today returns 403. Adjusted (non-OCC) roots are skipped.

The daily path and the backfill agree contract by contract on AAPL
2026-10-09 (1,068 contracts, identical volume).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

CLAUDE.md

@@ -110,7 +110,9 @@ PYTHONPATH=. uv run --with numba python -m research.lab.dip_buying   # buy-the-d
 PYTHONPATH=. uv run python -m research.lab.dip_leverage_overlay     # +1x SPY only inside SPY RSI(2) dip windows, on the overnight_book books (preregistered; STOP 2026-10-12: Sharpe 0.98->0.99 / 1.38->1.33, placebo rank 77-88%; the books already earn the bounce at ~2x SPY beta)
 PYTHONPATH=. uv run --with numba python -m research.lab.factor_census [--part construction]  # every lab factor (price, EDGAR, insider, Lazy Prices, EAR, SVR, FTD) measured one way on the current PIT panel: IC, quintile spread, top-20 alpha vs EW (2026-10-12: 1 A low_svr21, 2 B accruals/rev_5, 69 N; no |t|>=3)
 PYTHONPATH=. uv run --with numba python -m research.lab.svr_combo     # low_svr21 blended into overnight_recovery_top10 / as a fourth book sleeve, 2018-26 (preregistered; STOP 2026-10-12: blend Sharpe 0.81->0.87 but MaxDD -35%->-41%; book 1.49->1.42)
-PYTHONPATH=. uv run python -m research.lab.fetch_options snapshot|status  # daily S&P 500 option chains after the close (Alpaca indicative feed, forward-only from 2026-10-12; data/options/us/; timer quant-trading-options)
+PYTHONPATH=. uv run python -m research.lab.fetch_options snapshot|status  # daily S&P 500 option chains after the close: indicative quotes (forward-only from 2026-10-12) + OPRA daily bars + print-vs-quote quality sample; data/options/us/; timer quant-trading-options
+PYTHONPATH=. uv run python -m research.lab.fetch_options backfill        # OPRA daily option bars since 2024-02 for every S&P 500 member since then (resumable, ~2-3 h; pauses 16:10-16:55 NY); then builds volume_panel.csv.gz
+PYTHONPATH=. uv run python -m research.lab.fetch_options volume-panel    # date x underlying option volume (all 0-100 DTE strikes and the 5-100 DTE / 0.7-1.3 band), os_volume, pc_volume
 uv run python -m research.lab.fetch_us_fundamentals --universe mid      # EDGAR PIT facts for S&P 400 ever-members (name-checked CIKs)
 uv run python -m research.lab.fetch_us_fundamentals --universe small    # the same for S&P 600 (IJR)
 PYTHONPATH=. uv run python -m research.lab.midcap_sleeve [--only h1..h4] [--universe small]  # S&P 400 sleeves vs IJH; S&P 600 replication vs IJR (all REJECT 2026-10-08)
@@ -220,13 +222,15 @@ Four `systemd --user` services and one timer, defined in `deploy/systemd/`, run
 | `quant-trading-api` | `deploy/start-web-api.sh` | Web API on `127.0.0.1:8787` |
 | `quant-trading-vite` | `deploy/start-web-vite.sh` | Vite on `127.0.0.1:5174`, behind nginx on `:5443` |
 | `quant-trading-overnight-paper` | `deploy/start-overnight-paper.sh` | `overnight_paper.py daemon`: overnight-only momentum on the Alpaca paper account (from 2026-10-12) |
-| `quant-trading-options` (timer, 16:25 New York Mon-Fri) | `deploy/start-options-snapshot.sh` | `research.lab.fetch_options snapshot`: S&P 500 option chains + per-name IV/skew/spread summary into `data/options/us/` (~11 MB/day) |
+| `quant-trading-options` (timer, 16:25 New York Mon-Fri) | `deploy/start-options-snapshot.sh` | `research.lab.fetch_options snapshot`: S&P 500 option chains (indicative quotes), OPRA daily bars, per-name IV/skew/spread/volume summary and a quote-quality sample into `data/options/us/` (~15 MB/day); `status` prints the indicative-vs-OPRA verdict after 20 sessions |
 
 ```bash
 cp deploy/systemd/*.service ~/.config/systemd/user/
 systemctl --user daemon-reload
 systemctl --user enable --now quant-trading-monitor quant-trading-api quant-trading-vite
 systemctl --user enable --now quant-trading-options.timer
+# one-off, resumable: OPRA option-bar history since 2024-02 (rerun the same command if it stops)
+systemd-run --user --unit=quant-trading-options-backfill --working-directory=$PWD --setenv=PYTHONPATH=. uv run python -m research.lab.fetch_options backfill
 journalctl --user -u quant-trading-monitor -f
 ```
 

research/lab/fetch_options.py

@@ -1,63 +1,85 @@
 """
-Daily option-chain snapshots for the S&P 500, collected forward from
+Option data for the S&P 500, for options-implied stock factors:
-2026-10-12. Options-implied stock factors have no free history, so the
-record starts with this collector:
   - IV skew (Xing, Zhang & Zhao 2010)
   - the call-put IV spread (Cremers & Weinbaum 2010)
   - option/stock volume (Johnson & So 2012)
   - implied-vol level and changes
 
-    uv run python -m research.lab.fetch_options snapshot            # after the close (systemd timer)
+    uv run python -m research.lab.fetch_options snapshot          # after the close (systemd timer)
     uv run python -m research.lab.fetch_options snapshot --limit 5 --force --out /tmp/x   # smoke test
-    uv run python -m research.lab.fetch_options status
+    uv run python -m research.lab.fetch_options backfill          # OPRA daily bars since 2024-02 (resumable)
+    uv run python -m research.lab.fetch_options volume-panel      # daily option volume per underlying
+    uv run python -m research.lab.fetch_options status            # runs, backfill progress, quote-quality verdict
 
-Sources (keys ALPACA_API_KEY / ALPACA_API_SECRET from .env, never logged):
+Sources: Alpaca. Keys are ALPACA_API_KEY / ALPACA_API_SECRET from .env and
-  - stock close and volume: Alpaca SIP daily bars, unadjusted, so they
+are never logged. The account has the free Basic plan.
-    match the strikes. The free plan serves bars older than 15 minutes.
+  - Option quotes come from /v1beta1/options/snapshots, feed=indicative,
-  - option quotes: /v1beta1/options/snapshots/{underlying}, feed=indicative.
+    with Alpaca's IV and greeks. Alpaca's docs say indicative quotes "are not
-    This is Alpaca's free derived OPRA feed. It gives quote, last trade,
+    actual OPRA quotes". The OPRA BBO needs Algo Trader Plus ($99/month)
-    today's bar, and Alpaca's own IV and greeks. The real OPRA feed returns
+    plus a signed OPRA agreement. Quotes have no free history, so they are
-    403 until an OPRA agreement is signed.
+    collected forward from 2026-10-12.
-  - open interest: /v2/options/contracts (OCC, as of the prior day).
+  - Option bars come from /v1beta1/options/bars. They are computed from
-  The filter keeps expiries 5-100 calendar days out and strikes at
+    OPRA trades, so they are real prints. History starts in 2024-02, and
-  0.7-1.3 x the close. The universe is data/universe_sp500.json plus SPY,
+    Basic cannot read the last 15 minutes. Volume factors use these bars,
-  QQQ and IWM.
+    backfilled.
+  - Open interest comes from /v2/options/contracts (OCC, prior day); a
+    null means none was reported. Stock prices are SIP daily bars,
+    unadjusted, so they match the strikes.
 
-When: 16:25 America/New_York, Mon-Fri (deploy/systemd/quant-trading-options.timer).
+snapshot (16:25 America/New_York, Mon-Fri, deploy/systemd/quant-trading-options.timer)
-The Alpaca calendar decides the session. The script records the latest
+  The Alpaca calendar decides the session: the latest one that closed at
-session that closed at least 10 minutes ago. It skips the run if the next
+  least 10 minutes ago. The run is skipped if the next session has opened
-session has already opened (the quotes would belong to another day) and if
+  or the session's chain file exists, unless --force.
-the session's file exists, unless --force.
+  universe     data/universe_sp500.json + SPY, QQQ, IWM
+  quotes       expiries 5-100 days, strikes 0.7-1.3 x close
+  bars + OI    every contract expiring within 0-100 days, all strikes
+  quality      for each underlying, the ATM call, the ATM put and the
+               0.95-moneyness put of the expiry nearest 30 days (14-45 days)
+               are compared with OPRA prints. The print is the last 1-minute
+               bar from 15:45-16:00, and the reference is the indicative bid,
+               ask and mid. Measured: share of prints inside the quote,
+               position in the spread, IV(print) - IV(mid).
+               Pre-set rule: once >= 20 sessions are in, indicative quotes are
+               adequate if >= 90% of prints are inside the quote and the
+               median |IV(print) - IV(mid)| <= 1 vol point. Otherwise
+               consider the OPRA plan. `status` prints the verdict.
+backfill
+  - Universe: every name that was in the S&P 500 at some point since
+    2024-02-01 (data/sp500_history.json), plus the three ETFs, so names that
+    left the index are kept.
+  - Contracts: expired and active, expiring from 2024-02-01 to end+100 days.
+  - Bars: daily, from 100 days before each expiry to its expiry, for every
+    strike.
+  - One file per underlying, so the run resumes. It pauses 16:10-16:55 New
+    York on weekdays so it does not compete with snapshot for the
+    200/minute limit.
 
 Files (data/options/us/):
-  chains/YYYY/chain_<session>.csv.gz   one row per contract: quote, last trade,
+  chains/YYYY/chain_<session>.csv.gz   quotes, OPRA bar, OI, Alpaca IV/greeks, close,
-                                        day volume, OI, Alpaca IV/greeks, stock
+                                        r, implied forward, own IVs (iv_fwd, iv_spot)
-                                        close/volume, r, implied forward F, and
+  summary/summary_<session>.csv        per underlying: iv_atm30, skew_xzz, cp_spread,
-                                        two IVs from this script's own solver
+                                        option volume (OPRA), os_volume, pc ratios
-  summary/summary_<session>.csv        one row per underlying: iv_atm30,
+  bars/YYYY/bars_<session>.csv.gz      OPRA daily bars, 0-100 day expiries
-                                        skew_xzz, cp_spread, os_volume,
+  bars_hist/<UNDERLYING>.csv.gz        the backfill, same columns
-                                        pc_volume, pc_oi, quote quality
+  stock/stock_<session>.csv, stock_raw.csv.gz   unadjusted stock close/volume
-  manifest.jsonl                       one line per run
+  quality/quality_<session>.csv        the print-vs-quote sample
-The IVs are European. iv_fwd is Black-76 on the forward implied by put-call
+  volume_panel.csv.gz                  date x underlying option volume, built from bars
-parity at the three strikes nearest the close, so dividends drop out. iv_spot
+  manifest.jsonl, backfill_manifest.jsonl
-is Black-Scholes on the close with no dividends. cp_spread uses iv_spot,
+IVs are European. iv_fwd is Black-76 on the forward implied by put-call
-which is what Cremers-Weinbaum measure, so a dividend payer shows a negative
+parity at the three strikes nearest the close, so dividends drop out.
-spread of about its yield; adjust for that when testing. American early
+iv_spot is Black-Scholes on the close with no dividends; cp_spread uses it
-exercise is ignored, so the summary uses near-the-money strikes. The raw
+(as Cremers-Weinbaum do), so dividend payers read negative. American early
-chains are kept so a better solver can be run later.
+exercise is ignored, so the summaries use near-the-money strikes.
-
-NEEDS EVIDENCE: whether the indicative quotes match the OPRA NBBO closely
-enough for skew and spread. Compare a sample against another source before
-the first factor test. A cross-sectional test needs about 2-3 years of
-snapshots.
 """
 
 from __future__ import annotations
 
 import argparse
 import datetime as dt
+import glob
 import json
 import os
+import re
 import time
 from zoneinfo import ZoneInfo
 
@@ -73,28 +95,47 @@ DATA = "https://data.alpaca.markets"
 TRADE = "https://paper-api.alpaca.markets"
 OUT = os.path.join("data", "options", "us")
 UNIVERSE = os.path.join("data", "universe_sp500.json")
+HISTORY = os.path.join("data", "sp500_history.json")
 RF_PATH = os.path.join("data", "research", "ext", "fred_dtb3.csv")
 EXTRA = ["SPY", "QQQ", "IWM"]
 DTE_MIN, DTE_MAX = 5, 100
 MNY_LO, MNY_HI = 0.7, 1.3
-MIN_INTERVAL = 0.32                      # <= ~190 requests a minute (limit 200)
+HIST_START = "2024-02-01"
+QUALITY_MIN_SESSIONS, QUALITY_INSIDE, QUALITY_IV = 20, 0.90, 0.01
+OCC = re.compile(r"^[A-Z]{1,5}\d{6,7}[CP]\d{8}$")   # what the bars endpoint accepts; adjusted roots like 1AAPL fail
+BAR_COLS = ["date", "underlying", "symbol", "type", "expiry", "strike", "o", "h", "l", "c", "v", "n", "vw"]
 
 
 class Api:
-    def __init__(self) -> None:
+    """Paced GET with retries. ``yield_to_snapshot`` sleeps through 16:10-16:55
+    New York on weekdays so a backfill leaves the rate limit to the daily run."""
+
+    def __init__(self, min_interval: float = 0.32, yield_to_snapshot: bool = False) -> None:
         self.s = requests.Session()
         self.s.headers.update(_keys())
+        self.min_interval = min_interval
+        self.yield_to_snapshot = yield_to_snapshot
         self.next = 0.0
         self.n = 0
 
+    def _quiet_window(self) -> None:
+        now = dt.datetime.now(NY)
+        start = now.replace(hour=16, minute=10, second=0, microsecond=0)
+        end = now.replace(hour=16, minute=55, second=0, microsecond=0)
+        if now.weekday() < 5 and start <= now < end:
+            print(f"  pausing until 16:55 New York for the daily snapshot", flush=True)
+            time.sleep((end - now).total_seconds())
+
     def get(self, url: str, params: dict) -> dict:
-        for attempt in range(6):
+        for attempt in range(7):
+            if self.yield_to_snapshot:
+                self._quiet_window()
             wait = self.next - time.time()
             if wait > 0:
                 time.sleep(wait)
-            self.next = time.time() + MIN_INTERVAL
+            self.next = time.time() + self.min_interval
             try:
-                r = self.s.get(url, params=params, timeout=60)
+                r = self.s.get(url, params=params, timeout=90)
             except requests.RequestException:
                 time.sleep(2 ** attempt)
                 continue
@@ -119,7 +160,7 @@ class Api:
 
 
 # ---------------------------------------------------------------------------
-# session and inputs
+# inputs
 # ---------------------------------------------------------------------------
 
 def pick_session(api: Api, now: dt.datetime, force: bool) -> dt.date | None:
@@ -149,36 +190,65 @@ def universe() -> list[str]:
     return sorted({t.replace("-", ".") for t in ticks} | set(EXTRA))
 
 
-def stock_bars(api: Api, syms: list[str], session: dt.date, now: dt.datetime) -> dict[str, dict]:
+def pit_universe(start: str, end: str) -> list[str]:
-    end = min(now - dt.timedelta(minutes=16), dt.datetime.combine(session + dt.timedelta(days=1), dt.time(0), NY))
+    """Every name in the S&P 500 at some point in [start, end], plus EXTRA."""
-    out = {}
+    with open(HISTORY) as fh:
+        iv = json.load(fh)["intervals"]
+    keep = {t for t, spans in iv.items()
+            if any((a or "0000") <= end and (b is None or b >= start) for a, b in spans)}
+    return sorted({t.replace("-", ".") for t in keep} | set(universe()))
+
+
+def stock_bars(api: Api, syms: list[str], start: str, end: dt.datetime) -> pd.DataFrame:
+    rows = []
     for i in range(0, len(syms), 100):
-        params = {"symbols": ",".join(syms[i:i + 100]), "timeframe": "1Day", "start": session.isoformat(),
+        params = {"symbols": ",".join(syms[i:i + 100]), "timeframe": "1Day", "start": start,
                   "end": end.isoformat(), "feed": "sip", "adjustment": "raw", "limit": 10000}
         for bars in api.pages(f"{DATA}/v2/stocks/bars", params, "bars"):
             for sym, bl in bars.items():
-                for b in bl:
+                rows += [{"date": b["t"][:10], "symbol": sym, "close": b["c"], "volume": b["v"]} for b in bl]
-                    if b["t"][:10] == session.isoformat():
+    return pd.DataFrame(rows, columns=["date", "symbol", "close", "volume"])
-                        out[sym] = {"close": b["c"], "volume": b["v"]}
-    return out
 
 
-def open_interest(api: Api, syms: list[str], session: dt.date) -> dict[str, dict]:
+def list_contracts(api: Api, syms: list[str], gte: str, lte: str, status: str = "active") -> list[dict]:
-    out = {}
+    out = []
-    lo, hi = session + dt.timedelta(days=DTE_MIN), session + dt.timedelta(days=DTE_MAX)
     for i in range(0, len(syms), 10):
-        params = {"underlying_symbols": ",".join(syms[i:i + 10]), "expiration_date_gte": lo.isoformat(),
+        params = {"underlying_symbols": ",".join(syms[i:i + 10]), "status": status,
-                  "expiration_date_lte": hi.isoformat(), "limit": 10000}
+                  "expiration_date_gte": gte, "expiration_date_lte": lte, "limit": 10000}
         try:
             for page in api.pages(f"{TRADE}/v2/options/contracts", params, "option_contracts"):
                 for c in page:
-                    out[c["symbol"]] = {"oi": float(c["open_interest"]) if c.get("open_interest") else np.nan,
+                    out.append({"symbol": c["symbol"], "underlying": c["underlying_symbol"],
-                                        "oi_date": c.get("open_interest_date")}
+                                "type": "C" if c["type"] == "call" else "P", "expiry": c["expiration_date"],
+                                "strike": float(c["strike_price"]),
+                                "oi": float(c["open_interest"]) if c.get("open_interest") else np.nan,
+                                "oi_date": c.get("open_interest_date")})
         except RuntimeError as exc:
-            print(f"open interest {syms[i]}..: {exc}", flush=True)
+            print(f"contracts {syms[i]}..: {exc}", flush=True)
     return out
 
 
+def option_bars(api: Api, symbols: list[str], start: str, end: str, timeframe: str = "1Day") -> list[dict]:
+    """OPRA-derived bars. Contracts adjusted for corporate actions (non-standard
+    roots) are skipped: the endpoint rejects their symbols."""
+    symbols = [x for x in symbols if OCC.match(x)]
+    # Basic cannot read the last 15 minutes, and a date-only end of today counts
+    # as all of today (HTTP 403 "OPRA agreement is not signed"): send a timestamp
+    # no later than 16 minutes ago.
+    e = (dt.datetime.combine(dt.date.fromisoformat(end), dt.time(0), NY) if len(end) == 10
+         else dt.datetime.fromisoformat(end))
+    end = min(e, dt.datetime.now(NY) - dt.timedelta(minutes=16)).isoformat()
+    rows = []
+    for i in range(0, len(symbols), 1000):
+        params = {"symbols": ",".join(symbols[i:i + 1000]), "timeframe": timeframe, "start": start, "end": end,
+                  "limit": 10000}
+        for bars in api.pages(f"{DATA}/v1beta1/options/bars", params, "bars"):
+            for sym, bl in bars.items():
+                rows += [{"symbol": sym, "t": b["t"], "o": b["o"], "h": b["h"], "l": b["l"], "c": b["c"],
+                          "v": b["v"], "n": b.get("n"), "vw": b.get("vw")} for b in bl]
+    return rows
+
+
 def chain(api: Api, sym: str, close: float, session: dt.date) -> list[dict]:
     params = {"feed": "indicative", "limit": 1000,
               "expiration_date_gte": (session + dt.timedelta(days=DTE_MIN)).isoformat(),
@@ -193,7 +263,8 @@ def chain(api: Api, sym: str, close: float, session: dt.date) -> list[dict]:
                          "expiry": f"20{occ[-15:-13]}-{occ[-13:-11]}-{occ[-11:-9]}",
                          "strike": int(occ[-8:]) / 1000, "bid": q.get("bp"), "ask": q.get("ap"),
                          "bid_size": q.get("bs"), "ask_size": q.get("as"), "quote_time": q.get("t"),
-                         "last": t.get("p"), "last_time": t.get("t"), "day_volume": b.get("v", 0) if today else 0,
+                         "ind_last": t.get("p"), "ind_last_time": t.get("t"),
+                         "ind_day_volume": b.get("v", 0) if today else 0,
                          "alp_iv": s.get("impliedVolatility"), "delta": g.get("delta"), "gamma": g.get("gamma"),
                          "theta": g.get("theta"), "vega": g.get("vega")})
     return rows
@@ -220,6 +291,8 @@ def black76(F, K, T, sig, call, disc):
 
 
 def implied_vol(price, F, K, T, call, disc, iters: int = 80) -> np.ndarray:
+    price, F, K, T, call, disc = (np.asarray(x, dtype=float if i != 4 else bool)
+                                  for i, x in enumerate((price, F, K, T, call, disc)))
     lo, hi = np.full(len(price), 1e-4), np.full(len(price), 5.0)
     intrinsic = disc * np.where(call, np.maximum(F - K, 0), np.maximum(K - F, 0))
     ok = np.isfinite(price) & (price > intrinsic + 1e-6) & (price < black76(F, K, T, hi, call, disc))
@@ -237,7 +310,7 @@ def add_iv(df: pd.DataFrame, r: float) -> pd.DataFrame:
     df["T"] = df["dte"] / 365.0
     disc = np.exp(-r * df["T"])
     df["fwd"] = np.nan
-    for (_, _), g in df.groupby(["underlying", "expiry"]):
+    for _, g in df.groupby(["underlying", "expiry"]):
         piv = g.pivot_table(index="strike", columns="type", values="mid")
         if {"C", "P"} <= set(piv.columns):
             piv = piv.dropna()
@@ -251,19 +324,18 @@ def add_iv(df: pd.DataFrame, r: float) -> pd.DataFrame:
     call = (df["type"] == "C").to_numpy()
     args = (df["strike"].to_numpy(), df["T"].to_numpy(), call, disc.to_numpy())
     df["iv_fwd"] = implied_vol(df["mid"].to_numpy(), df["fwd"].to_numpy(), *args)
-    spot_f = (df["stock_close"] * np.exp(r * df["T"])).to_numpy()
+    df["iv_spot"] = implied_vol(df["mid"].to_numpy(), (df["stock_close"] * np.exp(r * df["T"])).to_numpy(), *args)
-    df["iv_spot"] = implied_vol(df["mid"].to_numpy(), spot_f, *args)
     return df
 
 
-def summarize(g: pd.DataFrame) -> dict:
+def summarize(g: pd.DataFrame, vol: dict) -> dict:
     S = g["stock_close"].iloc[0]
     v = g[g["iv_fwd"].notna()]
     out = {"underlying": g["underlying"].iloc[0], "stock_close": S, "stock_volume": g["stock_volume"].iloc[0],
            "n_contracts": len(g), "n_valid": len(v),
            "med_rel_spread": float(((v["ask"] - v["bid"]) / v["mid"]).median()) if len(v) else np.nan}
     atm = {}
-    for e, x in v.groupby("expiry"):
+    for _, x in v.groupby("expiry"):
         F = x["fwd"].iloc[0]
         k = x["strike"].iloc[np.argmin(np.abs(x["strike"].to_numpy() - F))]
         atm[int(x["dte"].iloc[0])] = float(x.loc[x["strike"] == k, "iv_fwd"].mean())
@@ -297,7 +369,7 @@ def summarize(g: pd.DataFrame) -> dict:
                     else pd.Series(1.0, index=pv.index)
                 w = w if w.sum() > 0 else pd.Series(1.0, index=pv.index)
                 out["cp_spread"] = float(((pv[("iv_spot", "C")] - pv[("iv_spot", "P")]) * w).sum() / w.sum())
-    cv, pv_ = g.loc[g["type"] == "C", "day_volume"].sum(), g.loc[g["type"] == "P", "day_volume"].sum()
+    cv, pv_ = vol.get("C", np.nan), vol.get("P", np.nan)
     co, po = g.loc[g["type"] == "C", "oi"].sum(), g.loc[g["type"] == "P", "oi"].sum()
     out.update({"call_volume": cv, "put_volume": pv_, "call_oi": co, "put_oi": po,
                 "pc_volume": pv_ / cv if cv else np.nan, "pc_oi": po / co if co else np.nan,
@@ -305,10 +377,65 @@ def summarize(g: pd.DataFrame) -> dict:
     return out
 
 
+def quality(api: Api, df: pd.DataFrame, session: dt.date) -> tuple[pd.DataFrame, dict]:
+    """Last OPRA print in 15:45-16:00 vs the indicative quote, near the money."""
+    x = df[df["dte"].between(14, 45) & df["mid"].notna() & df["iv_fwd"].notna()].copy()
+    if x.empty:
+        return x, {}
+    x["m_f"] = (x["strike"] / x["fwd"] - 1).abs()
+    x["m_95"] = (x["strike"] / x["stock_close"] - 0.95).abs()
+    pick = []
+    for _, g in x.groupby("underlying"):
+        e = g.loc[(g["dte"] - 30).abs().idxmin(), "expiry"]
+        g = g[g["expiry"] == e]
+        for typ, col in (("C", "m_f"), ("P", "m_f"), ("P", "m_95")):
+            h = g[g["type"] == typ]
+            if len(h):
+                pick.append(h[col].idxmin())
+    s = x.loc[sorted(set(pick))].copy()
+    t0 = dt.datetime.combine(session, dt.time(15, 30), NY).isoformat()
+    t1 = dt.datetime.combine(session, dt.time(16, 0), NY).isoformat()
+    bars = pd.DataFrame(option_bars(api, s["symbol"].tolist(), t0, t1, "1Min"))
+    if bars.empty:
+        return s.iloc[:0], {}
+    bars["t"] = pd.to_datetime(bars["t"], utc=True).dt.tz_convert(NY)
+    last = bars.sort_values("t").groupby("symbol").tail(1).set_index("symbol")
+    s = s.join(last[["c", "t"]].rename(columns={"c": "print", "t": "print_time"}), on="symbol")
+    cut = pd.Timestamp(dt.datetime.combine(session, dt.time(15, 45), NY))
+    s = s[s["print_time"].notna() & (s["print_time"] >= cut)].copy()
+    if s.empty:
+        return s, {}
+    s["inside"] = (s["print"] >= s["bid"] - 0.005) & (s["print"] <= s["ask"] + 0.005)
+    s["pos"] = (s["print"] - s["mid"]) / (s["ask"] - s["bid"])
+    s["iv_print"] = implied_vol(s["print"], s["fwd"], s["strike"], s["T"], s["type"] == "C", np.exp(-s["r"] * s["T"]))
+    s["iv_err"] = s["iv_print"] - s["iv_fwd"]
+    s["rel_spread"] = (s["ask"] - s["bid"]) / s["mid"]
+    summ = {"q_n": int(len(s)), "q_inside": float(s["inside"].mean()), "q_abs_pos": float(s["pos"].abs().median()),
+            "q_abs_iv_err": float(s["iv_err"].abs().median()), "q_rel_spread": float(s["rel_spread"].median())}
+    cols = ["session", "underlying", "symbol", "type", "expiry", "dte", "strike", "stock_close", "fwd", "bid", "ask",
+            "quote_time", "print", "print_time", "inside", "pos", "iv_fwd", "iv_print", "iv_err", "rel_spread"]
+    return s[cols], summ
+
+
 # ---------------------------------------------------------------------------
 # commands
 # ---------------------------------------------------------------------------
 
+def _atomic_csv(df: pd.DataFrame, path: str) -> None:
+    os.makedirs(os.path.dirname(path), exist_ok=True)
+    df.to_csv(path + ".tmp", index=False, compression="gzip" if path.endswith(".gz") else None)
+    os.replace(path + ".tmp", path)
+
+
+def _bars_frame(rows: list[dict], meta: pd.DataFrame) -> pd.DataFrame:
+    if not rows:
+        return pd.DataFrame(columns=BAR_COLS)
+    b = pd.DataFrame(rows)
+    b["date"] = b.pop("t").str[:10]
+    b = b.join(meta.set_index("symbol")[["underlying", "type", "expiry", "strike"]], on="symbol")
+    return b[BAR_COLS]
+
+
 def snapshot(args) -> None:
     out = args.out or OUT
     now = dt.datetime.now(NY)
@@ -322,8 +449,17 @@ def snapshot(args) -> None:
         return
     syms = universe()[: args.limit] if args.limit else universe()
     t0 = time.time()
-    px = stock_bars(api, syms, session, now)
+    end = min(now - dt.timedelta(minutes=16), dt.datetime.combine(session + dt.timedelta(days=1), dt.time(0), NY))
-    oi = open_interest(api, [s for s in syms if s in px], session)
+    st = stock_bars(api, syms, session.isoformat(), end)
+    st = st[st["date"] == session.isoformat()]
+    _atomic_csv(st.assign(session=session.isoformat()), os.path.join(out, "stock", f"stock_{session}.csv"))
+    px = st.set_index("symbol")[["close", "volume"]].to_dict("index")
+    have = [s for s in syms if s in px]
+    # contracts expiring on the session itself may already be inactive at 16:25,
+    # and 0DTE carries much of the day's volume: list both, as the backfill does
+    meta = pd.DataFrame(list_contracts(api, have, session.isoformat(), (session + dt.timedelta(days=DTE_MAX)).isoformat())
+                        + list_contracts(api, have, session.isoformat(), session.isoformat(), "inactive"))
+    meta = meta.drop_duplicates("symbol")
     rows, failed = [], []
     for k, sym in enumerate(syms):
         if sym not in px:
@@ -335,53 +471,179 @@ def snapshot(args) -> None:
             failed.append([sym, str(exc)[:120]])
             continue
         for r in rs:
-            r.update(stock_close=px[sym]["close"], stock_volume=px[sym]["volume"], **oi.get(r["symbol"], {}))
+            r.update(stock_close=px[sym]["close"], stock_volume=px[sym]["volume"])
         rows += rs
-        if (k + 1) % 50 == 0:
+        if (k + 1) % 100 == 0:
             print(f"  {k + 1}/{len(syms)} underlyings, {len(rows)} contracts, {api.n} requests "
                   f"[{time.time() - t0:.0f}s]", flush=True)
     if not rows:
         raise SystemExit("no contracts collected")
+    bars_err = None
+    try:
+        bars = _bars_frame(option_bars(api, meta["symbol"].tolist(), session.isoformat(),
+                                       (session + dt.timedelta(days=1)).isoformat()), meta)
+        bars = bars[bars["date"] == session.isoformat()]
+    except RuntimeError as exc:
+        bars, bars_err = pd.DataFrame(columns=BAR_COLS), str(exc)[:200]
+    _atomic_csv(bars, os.path.join(out, "bars", str(session.year), f"bars_{session}.csv.gz"))
     df = pd.DataFrame(rows)
-    for c in ("oi", "oi_date"):
+    df = df.join(meta.set_index("symbol")[["oi", "oi_date"]], on="symbol")
-        if c not in df:
+    df = df.join(bars.set_index("symbol")[["c", "v", "vw", "n"]].rename(
-            df[c] = np.nan
+        columns={"c": "opra_close", "v": "opra_volume", "vw": "opra_vwap", "n": "opra_trades"}), on="symbol")
     df.insert(0, "session", session.isoformat())
     df["dte"] = (pd.to_datetime(df["expiry"]) - pd.Timestamp(session)).dt.days
     r, r_date = rate(session)
     df["r"] = r
     df = add_iv(df, r)
-    os.makedirs(os.path.dirname(path), exist_ok=True)
+    try:
-    df.drop(columns=["mid", "T"]).to_csv(path + ".tmp", index=False, compression="gzip")
+        qdf, qsum = quality(api, df, session)
-    os.replace(path + ".tmp", path)
+        _atomic_csv(qdf, os.path.join(out, "quality", f"quality_{session}.csv"))
-    summ = pd.DataFrame([summarize(g) for _, g in df.groupby("underlying")])
+    except RuntimeError as exc:
+        qsum = {"q_error": str(exc)[:200]}
+    _atomic_csv(df.drop(columns=["mid", "T"]), path)
+    vol = bars.groupby(["underlying", "type"])["v"].sum()
+    summ = pd.DataFrame([summarize(g, vol.get(u, pd.Series(dtype=float)).to_dict())
+                         for u, g in df.groupby("underlying")])
     summ.insert(0, "session", session.isoformat())
-    spath = os.path.join(out, "summary", f"summary_{session}.csv")
+    _atomic_csv(summ, os.path.join(out, "summary", f"summary_{session}.csv"))
-    os.makedirs(os.path.dirname(spath), exist_ok=True)
-    summ.to_csv(spath, index=False)
     rec = {"session": session.isoformat(), "started": dt.datetime.fromtimestamp(t0, NY).isoformat(timespec="seconds"),
            "finished": dt.datetime.now(NY).isoformat(timespec="seconds"), "underlyings": len(syms),
            "with_chain": int(df["underlying"].nunique()), "contracts": len(df),
-           "valid_iv": int(df["iv_fwd"].notna().sum()), "requests": api.n, "r": r, "r_date": r_date,
+           "valid_iv": int(df["iv_fwd"].notna().sum()), "bars": len(bars), "bars_error": bars_err,
-           "feed": "indicative", "failed": failed[:50], "n_failed": len(failed), "forced": bool(args.force)}
+           "requests": api.n, "r": r, "r_date": r_date, "feed": "indicative", **qsum,
+           "failed": failed[:50], "n_failed": len(failed), "forced": bool(args.force)}
     with open(os.path.join(out, "manifest.jsonl"), "a") as fh:
         fh.write(json.dumps(rec) + "\n")
-    print(f"session {session}: {rec['with_chain']}/{len(syms)} underlyings, {len(df)} contracts "
+    print(f"session {session}: {rec['with_chain']}/{len(syms)} underlyings, {len(df)} quoted contracts "
-          f"({rec['valid_iv']} with IV), {api.n} requests, {time.time() - t0:.0f}s, failed {len(failed)} -> {path}")
+          f"({rec['valid_iv']} with IV), {len(bars)} OPRA bars, quality {qsum}, {api.n} requests, "
+          f"{time.time() - t0:.0f}s, failed {len(failed)}")
+
+
+def backfill(args) -> None:
+    out = args.out or OUT
+    api = Api(min_interval=0.4, yield_to_snapshot=True)       # ~150 requests a minute
+    end = args.end or (dt.datetime.now(NY).date() - dt.timedelta(days=1)).isoformat()
+    names = args.only or pit_universe(args.start, end)
+    hist = os.path.join(out, "bars_hist")
+    os.makedirs(hist, exist_ok=True)
+    raw = os.path.join(out, "stock_raw.csv.gz")
+    if not os.path.exists(raw) or args.only:
+        st = stock_bars(api, names, args.start, dt.datetime.combine(dt.date.fromisoformat(end), dt.time(23, 59), NY))
+        if not args.only:
+            _atomic_csv(st, raw)
+        print(f"stock bars: {st['symbol'].nunique()} symbols, {len(st)} rows", flush=True)
+    todo = [u for u in names if not os.path.exists(os.path.join(hist, f"{u}.csv.gz"))]
+    print(f"backfill {args.start}..{end}: {len(names)} underlyings, {len(todo)} to do", flush=True)
+    lte_active = (dt.date.fromisoformat(end) + dt.timedelta(days=DTE_MAX)).isoformat()
+    for k, u in enumerate(todo):
+        try:
+            backfill_one(api, u, args.start, end, lte_active, hist, out, k, len(todo))
+        except RuntimeError as exc:                  # no file is written, so the next run retries it
+            print(f"  [{k + 1}/{len(todo)}] {u}: failed, will retry on the next run: {exc}", flush=True)
+    if not args.only and not args.no_panel:
+        volume_panel(args)
+
+
+def backfill_one(api: Api, u: str, start: str, end: str, lte_active: str, hist: str, out: str, k: int, n: int) -> None:
+    t0, n0 = time.time(), api.n
+    meta = pd.DataFrame(list_contracts(api, [u], start, end, "inactive")
+                        + list_contracts(api, [u], start, lte_active, "active"))
+    frames = []
+    if len(meta):
+        meta = meta.drop_duplicates("symbol")
+        meta["month"] = meta["expiry"].str[:7]
+        for _, g in meta.groupby("month"):
+            first, last = g["expiry"].min(), g["expiry"].max()
+            a = max(start, (dt.date.fromisoformat(first) - dt.timedelta(days=DTE_MAX)).isoformat())
+            b = min(end, last)
+            if a > b:
+                continue
+            f = _bars_frame(option_bars(api, g["symbol"].tolist(), a, (dt.date.fromisoformat(b)
+                                                                      + dt.timedelta(days=1)).isoformat()), g)
+            lo = (pd.to_datetime(f["expiry"]) - pd.Timedelta(days=DTE_MAX)).dt.strftime("%Y-%m-%d")
+            frames.append(f[(f["date"] >= lo) & (f["date"] <= f["expiry"]) & (f["date"] <= end)])
+    bars = pd.concat(frames) if frames else pd.DataFrame(columns=BAR_COLS)
+    _atomic_csv(bars, os.path.join(hist, f"{u}.csv.gz"))
+    rec = {"underlying": u, "contracts": int(len(meta)), "bars": int(len(bars)), "requests": api.n - n0,
+           "seconds": round(time.time() - t0, 1), "start": start, "end": end,
+           "finished": dt.datetime.now(NY).isoformat(timespec="seconds")}
+    with open(os.path.join(out, "backfill_manifest.jsonl"), "a") as fh:
+        fh.write(json.dumps(rec) + "\n")
+    print(f"  [{k + 1}/{n}] {u}: {rec['contracts']} contracts, {rec['bars']} bars, "
+          f"{rec['requests']} requests, {rec['seconds']}s (total {api.n})", flush=True)
+
+
+def volume_panel(args) -> None:
+    """Daily option volume per underlying from the OPRA bars (backfill + daily
+    files): all strikes with 0-100 days to expiry, and the band the quote
+    snapshot covers (5-100 days, 0.7-1.3 x the unadjusted close)."""
+    out = args.out or OUT
+    stock = [pd.read_csv(os.path.join(out, "stock_raw.csv.gz"))] if os.path.exists(os.path.join(out, "stock_raw.csv.gz")) else []
+    stock += [pd.read_csv(f).drop(columns=["session"], errors="ignore") for f in glob.glob(os.path.join(out, "stock", "*.csv"))]
+    stock = pd.concat(stock).drop_duplicates(["date", "symbol"], keep="last").rename(columns={"symbol": "underlying"})
+    files = sorted(glob.glob(os.path.join(out, "bars_hist", "*.csv.gz")))
+    daily = sorted(glob.glob(os.path.join(out, "bars", "*", "*.csv.gz")))
+    bman = os.path.join(out, "backfill_manifest.jsonl")
+    hist_end = max((json.loads(x)["end"] for x in open(bman) if x.strip()), default="") if os.path.exists(bman) else ""
+    parts = []
+    for f in files + daily:
+        b = pd.read_csv(f)
+        if b.empty:
+            continue
+        if f in daily:
+            b = b[b["date"] > hist_end]           # the backfill wins where both exist
+        b = b.merge(stock[["date", "underlying", "close"]], on=["date", "underlying"], how="left")
+        dte = (pd.to_datetime(b["expiry"]) - pd.to_datetime(b["date"])).dt.days
+        band = dte.between(DTE_MIN, DTE_MAX) & (b["strike"] / b["close"]).between(MNY_LO, MNY_HI)
+        b["v_band"] = b["v"].where(band, 0)
+        g = b.groupby(["date", "underlying", "type"])[["v", "v_band"]].sum().unstack("type").fillna(0)
+        g.columns = [f"{a}_{t}" for a, t in g.columns]
+        parts.append(g.reset_index())
+    if not parts:
+        print("no bars yet")
+        return
+    P = pd.concat(parts).groupby(["date", "underlying"]).sum().reset_index()
+    for c in ("v_C", "v_P", "v_band_C", "v_band_P"):
+        if c not in P:
+            P[c] = 0.0
+    P = P.merge(stock[["date", "underlying", "close", "volume"]], on=["date", "underlying"], how="left")
+    P["os_volume"] = (P["v_C"] + P["v_P"]) * 100 / P["volume"]
+    P["os_volume_band"] = (P["v_band_C"] + P["v_band_P"]) * 100 / P["volume"]
+    P["pc_volume"] = P["v_P"] / P["v_C"].replace(0, np.nan)
+    _atomic_csv(P.sort_values(["date", "underlying"]), os.path.join(out, "volume_panel.csv.gz"))
+    print(f"volume panel: {P['underlying'].nunique()} underlyings, {P['date'].nunique()} days "
+          f"({P['date'].min()}..{P['date'].max()}), {len(P)} rows")
 
 
 def status(args) -> None:
     out = args.out or OUT
     man = os.path.join(out, "manifest.jsonl")
-    if not os.path.exists(man):
+    if os.path.exists(man):
-        print("no snapshots yet")
+        recs = [json.loads(x) for x in open(man) if x.strip()]
-        return
+        print(f"snapshots: {len(recs)} runs, {len({r['session'] for r in recs})} sessions, "
-    recs = [json.loads(x) for x in open(man) if x.strip()]
+              f"{recs[0]['session']}..{recs[-1]['session']}")
-    print(f"{len(recs)} runs, {len({r['session'] for r in recs})} sessions, first {recs[0]['session']}, "
+        for r in recs[-5:]:
-          f"last {recs[-1]['session']}")
+            print(f"  {r['session']} {r['with_chain']}/{r['underlyings']} underlyings, {r['contracts']} contracts, "
-    for r in recs[-10:]:
+                  f"{r.get('bars', '-')} OPRA bars, failed {r['n_failed']}, {r['started']} -> {r['finished']}")
-        print(f"  {r['session']} {r['with_chain']}/{r['underlyings']} underlyings {r['contracts']} contracts "
+        q = [r for r in recs if r.get("q_n")]
-              f"{r['valid_iv']} IV, failed {r['n_failed']}, {r['started']} -> {r['finished']}")
+        q = list({r["session"]: r for r in q}.values())
+        if q:
+            ins = float(np.average([r["q_inside"] for r in q], weights=[r["q_n"] for r in q]))
+            err = float(np.median([r["q_abs_iv_err"] for r in q]))
+            print(f"quote quality over {len(q)} sessions: {ins:.1%} of prints inside the indicative quote, "
+                  f"median |IV(print) - IV(mid)| {err * 100:.2f} vol points")
+            if len(q) >= QUALITY_MIN_SESSIONS:
+                ok = ins >= QUALITY_INSIDE and err <= QUALITY_IV
+                print("  verdict: indicative quotes adequate" if ok else "  verdict: consider the OPRA plan")
+            else:
+                print(f"  verdict after {QUALITY_MIN_SESSIONS} sessions")
+    else:
+        print("snapshots: none yet")
+    bman = os.path.join(out, "backfill_manifest.jsonl")
+    if os.path.exists(bman):
+        b = [json.loads(x) for x in open(bman) if x.strip()]
+        print(f"backfill: {len(b)} underlyings done, {sum(x['bars'] for x in b):,} bars, last {b[-1]['underlying']} "
+              f"at {b[-1]['finished']}")
 
 
 def main() -> None:
@@ -393,6 +655,16 @@ def main() -> None:
     s.add_argument("--session", help="label the snapshot with this session date (testing)")
     s.add_argument("--out", help="output directory (default data/options/us)")
     s.set_defaults(fn=snapshot)
+    b = sub.add_parser("backfill")
+    b.add_argument("--start", default=HIST_START)
+    b.add_argument("--end", help="last date (default yesterday)")
+    b.add_argument("--only", nargs="*", help="these underlyings only (testing; no panel)")
+    b.add_argument("--no-panel", action="store_true")
+    b.add_argument("--out")
+    b.set_defaults(fn=backfill)
+    v = sub.add_parser("volume-panel")
+    v.add_argument("--out")
+    v.set_defaults(fn=volume_panel)
     t = sub.add_parser("status")
     t.add_argument("--out")
     t.set_defaults(fn=status)