merge: #696 — event study around code-family exposures: pfs.event_study + notebooks/event_study.py (refs #696)

notebooks/event_study.py

@@ -0,0 +1,555 @@
+import marimo
+
+__generated_with = "0.23.13"
+app = marimo.App(width="medium")
+
+
+@app.cell(hide_code=True)
+def _():
+    import marimo as mo
+
+    return (mo,)
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md("""
+    # Event study: what happens around a code-family exposure
+
+    The physician fee schedule is law, and `pfs.code_exposure` is its exposure calendar:
+    the dates on which a code became payable, was revalued, changed status, ended, or was
+    listed for telehealth, each anchored to the Federal Register paragraph that says so.
+    This notebook puts pre/post windows around one exposure at a time on three outcome
+    series — services billed (`pfs.utilization`), clinical-staff wages (`bls.oews`) and
+    health-care employment (`bls.ces`) — and, for utilization, compares the treated code
+    with a never-treated control set drawn from codes of the same RVU status and decile
+    band that had no exposure of their own in the window. A code with no pre-period of
+    its own — a new code, or a successor such as 99439 — is measured through its family's
+    summed series, so G2058 → 99439 reads as one continuous line rather than a birth.
+
+    **What the numbers are.** The utilization estimate is a difference in log changes:
+    the treated code's log ratio of post-window to pre-window mean services, minus the mean
+    of the same quantity over its controls; the interval is a 95% percentile bootstrap over
+    the control set. It is a descriptive contrast around a dated legal event, reported with
+    its uncertainty. It is not a causal estimate, and the chat never states one from it.
+    """)
+    return
+
+
+@app.cell(hide_code=True)
+def _():
+    # ── 0. Setup ──
+    import altair as alt
+    import polars as pl
+
+    from conf import connect, path
+
+    connect.theme()
+    NOTES = {}
+    _primary = path("db.aco")
+    REPLICA_PATH = _primary.parent / "replica" / f"{_primary.stem}.ro.duckdb"
+
+    def _open_replica():
+        try:
+            return connect.duckdb("aco", read_only=True)
+        except Exception as e:  # noqa: BLE001 — degrade, never crash the page
+            NOTES["replica"] = f"replica unavailable: {e}"
+            return None
+
+    def _open_bib():
+        try:
+            return connect.bib()
+        except Exception as e:  # noqa: BLE001
+            NOTES["bib"] = f"bibliography unavailable: {e}"
+            return None
+
+    con = _open_replica()
+    store = _open_bib()
+
+    def fr_md(item_key, p_id):
+        if not item_key:
+            return "_unanchored_"
+        if store is None:
+            return f"{item_key} ¶{p_id}"
+        try:
+            from bib.frlink import md_link
+
+            return md_link(
+                f"p-{p_id}", store=store, item_key=item_key, text=f"{item_key} ¶{p_id}"
+            )
+        except Exception:  # noqa: BLE001
+            return f"{item_key} ¶{p_id}"
+
+    def not_built(cmd):
+        return f"_Not built yet — run `{cmd}` and republish the replica._"
+
+    return NOTES, REPLICA_PATH, alt, con, fr_md, not_built, pl, store
+
+
+@app.cell(hide_code=True)
+def _(con, mo, not_built):
+    # ── 1. Pick an exposure ──
+    from pfs.codetables import read_exposures as _read_exposures
+
+    #: The fixture exposures from #693/#696 — the CCM 2015 and APCM 2025 panels the
+    #: nightly run must render, plus the other dated events the calendar was built on.
+    FIXTURES = [
+        ("99490", 2015, "becomes-payable"),
+        ("99487", 2017, "becomes-payable"),
+        ("99490", 2022, "revalued"),
+        ("99439", 2021, "becomes-payable"),
+        ("G2058", 2021, "ends"),
+        ("99424", 2022, "becomes-payable"),
+        ("G2211", 2024, "becomes-payable"),
+        ("G0556", 2025, "becomes-payable"),
+        ("99441", 2025, "ends"),
+    ]
+    _rows = {}
+    if con is not None:
+        try:
+            for _c, _y, _k in FIXTURES:
+                for _r in _read_exposures(con, code=_c):
+                    if _r.year == _y and _r.kind == _k:
+                        _rows[(_c, _y, _k)] = _r
+        except Exception:  # noqa: BLE001 — not built yet
+            _rows = {}
+    labels = {
+        f"{c} {k} {y}"
+        + (f" [{_rows[(c, y, k)].family}]" if (c, y, k) in _rows else ""): (c, y, k)
+        for c, y, k in FIXTURES
+    }
+    exposure_picker = mo.ui.dropdown(
+        options=list(labels), value=list(labels)[0], label="Exposure"
+    )
+    _intro = mo.md(
+        "## 1. Pick an exposure\n\n"
+        "Eight dated events from the exposure calendar. Each is a row of "
+        "`pfs.code_exposure` with an effective date, the code's RVU status and decile "
+        "band that year, and the FR paragraph it is anchored to."
+        + ("" if _rows else "\n\n" + not_built("stack pfs exposure --write"))
+    )
+    mo.vstack([_intro, exposure_picker])
+    return FIXTURES, exposure_picker, labels
+
+
+@app.cell(hide_code=True)
+def _(exposure_picker, labels):
+    code, year, kind = labels[exposure_picker.value]
+    return code, kind, year
+
+
+@app.cell(hide_code=True)
+def _(alt, code, con, fr_md, kind, math_exp, mo, not_built, pl, year):
+    # ── 2. Utilization ──
+    from pfs.event_study import utilization_panel as _panel
+
+    _p = None
+    if con is not None:
+        try:
+            _p = _panel(con, code, year, kind)
+        except Exception as _e:  # noqa: BLE001 — a missing table is "not built yet"
+            _p = None
+    _head = (
+        f"## 2. Utilization around {code} {kind} ({year})\n\n"
+        "Services billed per year for the treated code (solid) against its never-treated "
+        "controls, each control indexed to its own pre-window mean so the shapes are "
+        "comparable; the shaded band is the 10th–90th percentile of the indexed controls. "
+        "The table gives the log-change contrast and its bootstrap interval."
+    )
+    if _p is None:
+        _view = mo.md(_head + "\n\n" + not_built("stack pfs utilization --write"))
+    else:
+        _e = _p.effect
+        _anchor = (
+            fr_md(_p.exposure.item_key, _p.exposure.p_id)
+            if _p.exposure
+            else "_no exposure row_"
+        )
+        _eff = _p.exposure.effective_date.isoformat() if _p.exposure else "—"
+
+        def _fmt(v, pct=False):
+            if v is None:
+                return "—"
+            return f"{v:+.1%}" if pct else f"{v:,.0f}"
+
+        _tbl = pl.DataFrame(
+            {
+                "quantity": [
+                    "effective date",
+                    "anchor",
+                    "treated series",
+                    "pre window",
+                    "post window",
+                    "treated pre mean (services)",
+                    "treated post mean (services)",
+                    "treated change",
+                    "control change (mean)",
+                    "contrast (treated − control)",
+                    "95% interval",
+                    "controls",
+                ],
+                "value": [
+                    _eff,
+                    _anchor,
+                    f"{code} alone"
+                    if _e.scope == "code"
+                    else f"{_p.exposure.family} family total (the code had no pre-period of its own)",
+                    ", ".join(map(str, _e.pre_years)) or "— (none on file)",
+                    ", ".join(map(str, _e.post_years)) or "— (none on file)",
+                    _fmt(_e.treated_pre),
+                    _fmt(_e.treated_post),
+                    _fmt(math_exp(_e.treated_change), True)
+                    if _e.treated_change is not None
+                    else "— (no pre-period: the code did not exist or was not billed)",
+                    _fmt(math_exp(_e.control_change), True)
+                    if _e.control_change is not None
+                    else "—",
+                    _fmt(_e.effect_pct, True),
+                    f"{math_exp(_e.ci_low):+.1%} to {math_exp(_e.ci_high):+.1%}"
+                    if _e.ci_low is not None
+                    else "— (needs ≥2 controls and a treated change)",
+                    str(_e.n_controls),
+                ],
+            }
+        )
+        _panels = [
+            mo.md(_head),
+            mo.ui.table(_tbl, label=f"pfs.event_study — {code} {year}"),
+        ]
+        if _p.treated or _p.controls:
+            _recs = []
+
+            def _indexed(series, pre_years):
+                base = [series[y] for y in pre_years if y in series]
+                b = sum(base) / len(base) if base else None
+                return {y: (v / b if b else None) for y, v in series.items()}
+
+            for y, v in sorted(_p.treated.items()):
+                _recs.append(
+                    {
+                        "year": y,
+                        "series": f"{code} (treated)",
+                        "services": v,
+                        "kind": "treated",
+                    }
+                )
+            _ctrl_idx = {
+                _c: _indexed(_s, _e.pre_years) for _c, _s in _p.controls.items()
+            }
+            _by_year = {}
+            for _c, _s in _ctrl_idx.items():
+                for _y, _v in _s.items():
+                    if _v is not None:
+                        _by_year.setdefault(_y, []).append(_v)
+            _band = [
+                {
+                    "year": y,
+                    "p10": sorted(v)[int(0.1 * (len(v) - 1))],
+                    "p90": sorted(v)[int(0.9 * (len(v) - 1))],
+                    "median": sorted(v)[len(v) // 2],
+                }
+                for y, v in sorted(_by_year.items())
+            ]
+            _t = (
+                pl.DataFrame(_recs)
+                if _recs
+                else pl.DataFrame(
+                    {"year": [], "series": [], "services": [], "kind": []}
+                )
+            )
+            _treated_chart = (
+                alt.Chart(_t.to_pandas())
+                .mark_line(point=True)
+                .encode(
+                    x=alt.X("year:O", title="Calendar year"),
+                    y=alt.Y("services:Q", title="Services (treated)"),
+                    tooltip=["year", "services"],
+                )
+                .properties(width=640, height=220, title=f"{code}: services billed")
+            )
+            _rule = (
+                alt.Chart(pl.DataFrame({"year": [year]}).to_pandas())
+                .mark_rule(strokeDash=[4, 4], color="gray")
+                .encode(x="year:O")
+            )
+            _panels.append(mo.ui.altair_chart(_treated_chart + _rule))
+            if _band:
+                _b = pl.DataFrame(_band)
+                _band_chart = (
+                    alt.Chart(_b.to_pandas())
+                    .mark_area(opacity=0.25)
+                    .encode(
+                        x=alt.X("year:O", title="Calendar year"),
+                        y=alt.Y("p10:Q", title="Indexed to pre-window mean"),
+                        y2="p90:Q",
+                    )
+                    .properties(
+                        width=640,
+                        height=220,
+                        title=f"Controls ({_e.n_controls}): indexed services, 10th–90th pct and median",
+                    )
+                )
+                _med = (
+                    alt.Chart(_b.to_pandas())
+                    .mark_line()
+                    .encode(x="year:O", y="median:Q")
+                )
+                _panels.append(mo.ui.altair_chart(_band_chart + _med + _rule))
+        else:
+            _panels.append(
+                mo.md(
+                    "_No utilization rows on file for this code or its controls yet — for a 2025 exposure that means the 2025 public use file has not been published._"
+                )
+            )
+        _view = mo.vstack(_panels)
+    _view
+    return
+
+
+@app.cell(hide_code=True)
+def _():
+    import math as _math
+
+    def math_exp(v):
+        return _math.exp(v) - 1 if v is not None else None
+
+    return (math_exp,)
+
+
+@app.cell(hide_code=True)
+def _(alt, code, con, kind, mo, not_built, pl, year):
+    # ── 3. Wages and employment ──
+    from bls.ces import read_series as _read_ces
+    from bls.oews import read_occupations as _read_oews
+    from pfs.event_study import pre_post as _pre_post
+
+    _occs = {
+        "31-9092": "Medical assistants",
+        "29-2061": "LPN/LVN",
+        "29-1141": "Registered nurses",
+        "29-1171": "Nurse practitioners",
+    }
+    _series_ids = {
+        "CES6562110001": "Offices of physicians, employment (k)",
+        "CES6562100001": "Ambulatory care, employment (k)",
+        "CES6562110003": "Offices of physicians, avg hourly earnings ($)",
+    }
+    _w, _c = [], []
+    if con is not None:
+        try:
+            _w = _read_oews(con, list(_occs))
+        except Exception:  # noqa: BLE001
+            _w = []
+        try:
+            _c = _read_ces(con, list(_series_ids), annual=True)
+        except Exception:  # noqa: BLE001
+            _c = []
+    _head = (
+        f"## 3. Clinical-staff wages and employment around {year}\n\n"
+        "These series are national, not code-specific, so there is no control set: the "
+        "table reports the pre- and post-window means and the log change for each, "
+        "which is context for the utilization contrast above, not an estimate of the "
+        "exposure's effect."
+    )
+    if not _w and not _c:
+        _view = mo.md(
+            _head
+            + "\n\n"
+            + not_built("stack bls oews --write` and `stack bls ces --write")
+        )
+    else:
+        _rows = []
+        for occ, name in _occs.items():
+            s = {
+                r.year: r.a_mean
+                for r in _w
+                if r.occ_code == occ and r.a_mean is not None
+            }
+            if s:
+                p = _pre_post(s, year)
+                _rows.append(
+                    {
+                        "series": f"{occ} {name}: mean annual wage ($)",
+                        "pre": ", ".join(map(str, p.pre_years)),
+                        "pre mean": p.pre_mean,
+                        "post": ", ".join(map(str, p.post_years)),
+                        "post mean": p.post_mean,
+                        "log change": p.change,
+                    }
+                )
+        for sid, name in _series_ids.items():
+            s = {
+                r.year: r.value
+                for r in _c
+                if r.series_id == sid and r.value is not None
+            }
+            if s:
+                p = _pre_post(s, year)
+                _rows.append(
+                    {
+                        "series": f"{sid} {name}",
+                        "pre": ", ".join(map(str, p.pre_years)),
+                        "pre mean": p.pre_mean,
+                        "post": ", ".join(map(str, p.post_years)),
+                        "post mean": p.post_mean,
+                        "log change": p.change,
+                    }
+                )
+        _tbl = pl.DataFrame(_rows) if _rows else pl.DataFrame({"series": []})
+        _wdf = pl.DataFrame(
+            {
+                "occupation": [f"{r.occ_code} {_occs[r.occ_code]}" for r in _w],
+                "year": [r.year for r in _w],
+                "mean annual wage": [r.a_mean for r in _w],
+            }
+        )
+        _rule = (
+            alt.Chart(pl.DataFrame({"year": [year]}).to_pandas())
+            .mark_rule(strokeDash=[4, 4], color="gray")
+            .encode(x="year:O")
+        )
+        _panels = [mo.md(_head), mo.ui.table(_tbl, label="pre/post window means")]
+        if len(_wdf):
+            _chart = (
+                alt.Chart(_wdf.to_pandas())
+                .mark_line(point=True)
+                .encode(
+                    x=alt.X("year:O", title="OEWS release year"),
+                    y=alt.Y("mean annual wage:Q", title="Mean annual wage ($)"),
+                    color=alt.Color("occupation:N", title="Occupation"),
+                    tooltip=["occupation", "year", "mean annual wage"],
+                )
+                .properties(width=640, height=260, title="Clinical-staff wages (OEWS)")
+            )
+            _panels.append(mo.ui.altair_chart(_chart + _rule))
+        _view = mo.vstack(_panels)
+    _view
+    return
+
+
+@app.cell(hide_code=True)
+def _(FIXTURES, con, fr_md, mo, not_built, pl):
+    # ── 4. All fixture exposures ──
+    from pfs.event_study import utilization_panel as _panel_all
+
+    _rows = []
+    if con is not None:
+        for _c, _y, _k in FIXTURES:
+            try:
+                _pn = _panel_all(con, _c, _y, _k)
+            except Exception:  # noqa: BLE001
+                continue
+            _pe = _pn.effect
+            _rows.append(
+                {
+                    "exposure": f"{_c} {_k} {_y}",
+                    "family": _pn.exposure.family if _pn.exposure else "",
+                    "treated": _pe.scope,
+                    "anchor": fr_md(_pn.exposure.item_key, _pn.exposure.p_id)
+                    if _pn.exposure
+                    else "—",
+                    "pre": ", ".join(map(str, _pe.pre_years)),
+                    "post": ", ".join(map(str, _pe.post_years)),
+                    "contrast": f"{_pe.effect_pct:+.1%}"
+                    if _pe.effect_pct is not None
+                    else "—",
+                    "95% interval": f"{__import__('math').exp(_pe.ci_low) - 1:+.1%} to {__import__('math').exp(_pe.ci_high) - 1:+.1%}"
+                    if _pe.ci_low is not None
+                    else "—",
+                    "controls": _pe.n_controls,
+                    "note": ""
+                    if _pe.effect is not None
+                    else (
+                        "no pre-period (new code)"
+                        if _pe.treated_change is None
+                        and _pe.pre_years
+                        and not _pe.treated_pre
+                        else (
+                            "no post-period on file yet"
+                            if not _pe.post_years or _pe.treated_post is None
+                            else "insufficient controls"
+                        )
+                    ),
+                }
+            )
+    _head = (
+        "## 4. All fixture exposures\n\n"
+        "The same contrast for every fixture event, with the FR anchor that dates it. A new "
+        "code has no pre-period by construction — its first years are the outcome, not a "
+        "change — and a 2025 exposure has no post-period until CMS publishes the 2025 file; "
+        "both are stated rather than filled in."
+    )
+    _view = (
+        mo.vstack(
+            [
+                mo.md(_head),
+                mo.ui.table(pl.DataFrame(_rows), label="event-study contrasts"),
+            ]
+        )
+        if _rows
+        else mo.md(
+            _head
+            + "\n\n"
+            + not_built(
+                "stack pfs exposure --write` and `stack pfs utilization --write"
+            )
+        )
+    )
+    _view
+    return
+
+
+@app.cell(hide_code=True)
+def _(NOTES, REPLICA_PATH, con, mo, pl):
+    # ── 5. Provenance ──
+    import datetime as _dt
+
+    if REPLICA_PATH.exists():
+        _mtime = _dt.datetime.fromtimestamp(REPLICA_PATH.stat().st_mtime).isoformat(
+            timespec="seconds"
+        )
+        _line = f"`{REPLICA_PATH}` — last published {_mtime}"
+    else:
+        _line = f"`{REPLICA_PATH}` — not found"
+
+    def _q(sql):
+        if con is None:
+            return pl.DataFrame()
+        try:
+            return con.execute(sql).pl()
+        except Exception as e:  # noqa: BLE001
+            NOTES[sql[:40]] = str(e)
+            return pl.DataFrame()
+
+    _counts = _q(
+        "SELECT 'pfs.code_exposure' t, count(*) n FROM pfs.code_exposure "
+        "UNION ALL SELECT 'pfs.utilization', count(*) FROM pfs.utilization "
+        "UNION ALL SELECT 'bls.oews', count(*) FROM bls.oews "
+        "UNION ALL SELECT 'bls.ces', count(*) FROM bls.ces"
+    )
+    _log = _q(
+        "SELECT ingested_at, module, table_name, rows, pincite_key FROM cms.ingest_log "
+        "WHERE module IN ('pfs.utilization', 'bls.oews', 'bls.ces') ORDER BY ingested_at DESC LIMIT 30"
+    )
+    _notes = "\n".join(f"- {k}: {v}" for k, v in NOTES.items())
+    mo.vstack(
+        [
+            mo.md(
+                "## 5. Provenance\n\n"
+                f"Replica: {_line}. Every outcome row carries the bib key of the dataset "
+                "release it came from (`item_key`), and each release is a Source item "
+                "tagged with its table; `cms.ingest_log` records the ingest with the source "
+                "file's sha256. Market and industry series (sector ETFs, telehealth and "
+                "care-management tickers) are not on file: no keyless quote source is "
+                "reachable from the rack, so that panel waits on a data source (#696 follow-up)."
+                + (f"\n\n{_notes}" if _notes else "")
+            ),
+            mo.ui.table(_counts, label="row counts")
+            if len(_counts)
+            else mo.md("_row counts unavailable_"),
+            mo.ui.table(_log, label="cms.ingest_log — outcome series")
+            if len(_log)
+            else mo.md("_no ingest log rows_"),
+        ]
+    )
+    return

src/pfs/event_study.py

@@ -0,0 +1,265 @@
+"""Event study around code-family exposures (#696): pre/post windows on an
+outcome series with never-treated control sets.
+
+The fee schedule is law; ``pfs.code_exposure`` is the exposure calendar
+(P50). For an exposure in *year*, the window is ``pre`` file years before
+it and ``post`` years from it on (the effective date is 1 January, so
+the exposure year is the first treated year). The estimate is a
+difference in log changes::
+
+    effect = log(post_mean/pre_mean)[treated] − mean_c log(post_mean/pre_mean)[control c]
+
+with a 95% percentile-bootstrap interval from resampling the control set
+(``_BOOT`` draws, fixed seed — the same inputs always give the same
+interval). It is a descriptive contrast, not a causal estimate: the
+notebook says so and the chat only links here.
+
+    windows(year, available, pre, post)  → (pre_years, post_years)
+    log_change(series, pre, post)        → log ratio of window means
+    did_effect(treated, controls, …)     → Effect
+    pre_post(series, year, …)            → PrePost (single-series, no controls)
+    utilization_panel(con, code, year, kind) → Panel over pfs.utilization
+    family_series(con, family)           → member codes summed per year
+"""
+
+from __future__ import annotations
+
+import math
+import random
+from dataclasses import dataclass
+from typing import Any, Iterable, Mapping, Sequence
+
+_BOOT = 500
+_SEED = 20260922
+
+
+@dataclass(frozen=True)
+class Effect:
+    code: str
+    year: int
+    kind: str
+    outcome: str
+    pre_years: tuple[int, ...]
+    post_years: tuple[int, ...]
+    treated_pre: float | None
+    treated_post: float | None
+    treated_change: float | None
+    control_change: float | None
+    effect: float | None
+    ci_low: float | None
+    ci_high: float | None
+    n_controls: int
+    #: ``"code"`` when the treated series is the code's own utilization;
+    #: ``"family"`` when the code had no pre-period of its own (a new
+    #: code, or a successor) and the family's summed series stood in —
+    #: G2058 → 99439 is then one continuous series, not a birth.
+    scope: str = "code"
+
+    @property
+    def effect_pct(self) -> float | None:
+        return math.exp(self.effect) - 1 if self.effect is not None else None
+
+
+@dataclass(frozen=True)
+class PrePost:
+    year: int
+    pre_years: tuple[int, ...]
+    post_years: tuple[int, ...]
+    pre_mean: float | None
+    post_mean: float | None
+    change: float | None
+
+
+@dataclass(frozen=True)
+class Panel:
+    exposure: Any  # ExposureRow | None
+    effect: Effect
+    treated: dict[int, float]
+    controls: dict[str, dict[int, float]]
+
+
+# ── pure pieces ───────────────────────────────────────────────────────
+
+
+def windows(
+    year: int, *, available: Iterable[int], pre: int = 2, post: int = 2
+) -> tuple[tuple[int, ...], tuple[int, ...]]:
+    have = set(available)
+    return (
+        tuple(y for y in range(year - pre, year) if y in have),
+        tuple(y for y in range(year, year + post) if y in have),
+    )
+
+
+def _mean(vals: Sequence[float]) -> float | None:
+    return sum(vals) / len(vals) if vals else None
+
+
+def _window_means(
+    series: Mapping[int, float], pre: Sequence[int], post: Sequence[int]
+) -> tuple[float | None, float | None]:
+    pv = [series[y] for y in pre if series.get(y) is not None]
+    qv = [series[y] for y in post if series.get(y) is not None]
+    return _mean(pv), _mean(qv)
+
+
+def log_change(
+    series: Mapping[int, float], pre: Sequence[int], post: Sequence[int]
+) -> float | None:
+    a, b = _window_means(series, pre, post)
+    if a is None or b is None or a <= 0 or b <= 0:
+        return None
+    return math.log(b / a)
+
+
+def did_effect(
+    treated: Mapping[int, float],
+    controls: Mapping[str, Mapping[int, float]],
+    year: int,
+    *,
+    code: str,
+    kind: str,
+    outcome: str,
+    pre: int = 2,
+    post: int = 2,
+) -> Effect:
+    years = set(treated) | {y for s in controls.values() for y in s}
+    pre_years, post_years = windows(year, available=years, pre=pre, post=post)
+    t_pre, t_post = _window_means(treated, pre_years, post_years)
+    t_change = log_change(treated, pre_years, post_years)
+    changes = [
+        c
+        for c in (log_change(s, pre_years, post_years) for s in controls.values())
+        if c is not None
+    ]
+    c_change = _mean(changes)
+    effect = (
+        t_change - c_change if t_change is not None and c_change is not None else None
+    )
+    lo = hi = None
+    if effect is not None and len(changes) >= 2:
+        rng = random.Random(_SEED)
+        n = len(changes)
+        draws = sorted(
+            t_change - sum(rng.choice(changes) for _ in range(n)) / n
+            for _ in range(_BOOT)
+        )
+        lo, hi = draws[int(0.025 * _BOOT)], draws[int(0.975 * _BOOT) - 1]
+    return Effect(
+        code=code,
+        year=year,
+        kind=kind,
+        outcome=outcome,
+        pre_years=pre_years,
+        post_years=post_years,
+        treated_pre=t_pre,
+        treated_post=t_post,
+        treated_change=t_change,
+        control_change=c_change,
+        effect=effect,
+        ci_low=lo,
+        ci_high=hi,
+        n_controls=len(changes),
+    )
+
+
+def pre_post(
+    series: Mapping[int, float], year: int, *, pre: int = 2, post: int = 2
+) -> PrePost:
+    """A single series' window means around *year* — for outcomes that
+    have no code-level control set (occupation wages, employment)."""
+    pre_years, post_years = windows(year, available=series, pre=pre, post=post)
+    a, b = _window_means(series, pre_years, post_years)
+    return PrePost(
+        year, pre_years, post_years, a, b, log_change(series, pre_years, post_years)
+    )
+
+
+# ── panels over the replica ───────────────────────────────────────────
+
+
+def utilization_panel(
+    con: Any,
+    code: str,
+    year: int,
+    kind: str,
+    *,
+    outcome: str = "n_services",
+    pre: int = 2,
+    post: int = 2,
+    window: int = 2,
+    max_controls: int = 60,
+) -> Panel:
+    """The exposure row, its ``pfs.utilization`` series, the never-treated
+    control series (``pfs.exposure.control_codes``) and the estimate."""
+    from pfs.codetables import read_exposures
+    from pfs.exposure import control_codes
+    from pfs.utilization import read_series
+
+    code = code.upper()
+    exposure = next(
+        (
+            r
+            for r in read_exposures(con, code=code)
+            if r.year == year and r.kind == kind
+        ),
+        None,
+    )
+    try:
+        ctrl = control_codes(con, code, year, window=window, limit=max_controls)
+    except Exception:  # noqa: BLE001 — no exposure calendar yet
+        ctrl = []
+    ctrl_codes = [c for c, *_ in ctrl]
+    rows = read_series(con, [code, *ctrl_codes])
+    series: dict[str, dict[int, float]] = {}
+    for r in rows:
+        v = getattr(r, outcome)
+        if v is not None:
+            series.setdefault(r.hcpcs, {})[r.year] = float(v)
+    treated = series.get(code, {})
+    controls = {c: series[c] for c in ctrl_codes if c in series}
+    scope = "code"
+    pre_years, _ = windows(year, available=range(1900, 2200), pre=pre, post=post)
+    if (
+        exposure is not None
+        and exposure.family
+        and not any(y in treated for y in pre_years)
+    ):
+        fam_series = family_series(con, exposure.family, outcome=outcome)
+        if any(y in fam_series for y in pre_years):
+            treated, scope = fam_series, "family"
+    effect = did_effect(
+        treated,
+        controls,
+        year,
+        code=code,
+        kind=kind,
+        outcome=outcome,
+        pre=pre,
+        post=post,
+    )
+    effect = Effect(**{**effect.__dict__, "scope": scope})
+    return Panel(exposure=exposure, effect=effect, treated=treated, controls=controls)
+
+
+def family_series(
+    con: Any, family: str, *, outcome: str = "n_services"
+) -> dict[int, float]:
+    """The family's member codes (``pfs.code_family``) summed per year."""
+    from pfs.utilization import read_series
+
+    try:
+        codes = [
+            r[0]
+            for r in con.execute(
+                "SELECT DISTINCT code FROM pfs.code_family WHERE key = ?", [family]
+            ).fetchall()
+        ]
+    except Exception:  # noqa: BLE001 — no family table on this replica
+        return {}
+    out: dict[int, float] = {}
+    for r in read_series(con, codes):
+        v = getattr(r, outcome)
+        if v is not None:
+            out[r.year] = out.get(r.year, 0.0) + float(v)
+    return out

tests/notebooks/test_event_study_nb.py

@@ -0,0 +1,72 @@
+"""notebooks/event_study.py — structure, headless degradation, live fixture panels (#696)."""
+
+from __future__ import annotations
+
+import ast
+import importlib.util
+import re
+from pathlib import Path
+
+import pytest
+
+NB = Path(__file__).resolve().parents[2] / "notebooks" / "event_study.py"
+_ROOT = Path(__file__).resolve().parents[2]
+_HAS_DATA = (_ROOT / "data" / "replica" / "aco.ro.duckdb").exists() and (
+    _ROOT / "data" / "bib.sqlite"
+).exists()
+_live = pytest.mark.skipif(not _HAS_DATA, reason="needs the live replica and bib")
+
+
+def _load():
+    spec = importlib.util.spec_from_file_location("event_study_nb", NB)
+    mod = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(mod)
+    return mod
+
+
+def test_notebook_is_a_marimo_app():
+    assert _load().app.__class__.__name__ == "App"
+
+
+def test_cells_are_anonymous_and_banners_present():
+    src = NB.read_text()
+    names = [n.name for n in ast.parse(src).body if isinstance(n, ast.FunctionDef)]
+    assert names and set(names) == {"_"}
+    assert re.findall(r"# ── (\S+)\. ", src) == ["0", "1", "2", "3", "4", "5"]
+
+
+def test_no_causal_language():
+    src = NB.read_text().lower()
+    assert (
+        "caused" not in src and "causal estimate" in src
+    )  # the disclaimer, nothing else
+
+
+def test_headless_run_degrades_without_data(monkeypatch):
+    monkeypatch.setenv("STACK_DUCKDB_REPLICA", "1")
+    mod = _load()
+    import conf.connect as cc
+
+    monkeypatch.setattr(
+        cc,
+        "duckdb",
+        lambda *a, **k: (_ for _ in ()).throw(FileNotFoundError("no replica")),
+    )
+    monkeypatch.setattr(
+        cc, "bib", lambda *a, **k: (_ for _ in ()).throw(FileNotFoundError("no bib"))
+    )
+    outputs, _ = mod.app.run()
+    rendered = "\n".join(o._repr_html_() for o in outputs if hasattr(o, "_repr_html_"))
+    assert "Not built yet" in rendered and "Traceback" not in rendered
+
+
+@_live
+def test_ccm_2015_and_apcm_2025_panels_render_with_anchors():
+    mod = _load()
+    outputs, _ = mod.app.run()
+    rendered = "\n".join(o._repr_html_() for o in outputs if hasattr(o, "_repr_html_"))
+    assert "Traceback" not in rendered
+    # section 4 lists every fixture with its anchor: CCM 2015 (DE2VH9PD) and APCM 2025 (JJ6AM5HJ)
+    assert "99490 becomes-payable 2015" in rendered and "DE2VH9PD" in rendered
+    assert "G0556 becomes-payable 2025" in rendered and "JJ6AM5HJ" in rendered
+    assert "no post-period on file yet" in rendered

tests/pfs/test_event_study.py

@@ -0,0 +1,263 @@
+"""pfs.event_study — pre/post windows around exposures with never-treated controls (#696)."""
+
+from __future__ import annotations
+
+import datetime as dt
+import math
+
+import duckdb
+import pytest
+
+from pfs.codetables import ExposureRow, ensure_tables, write_exposures
+from pfs.event_study import (
+    Effect,
+    PrePost,
+    did_effect,
+    log_change,
+    pre_post,
+    utilization_panel,
+    windows,
+)
+from pfs.utilization import UtilizationRow
+from pfs.utilization import write as write_util
+
+
+class TestWindows:
+    def test_pre_before_post_from_exposure_year(self):
+        pre, post = windows(2017, available=range(2013, 2025))
+        assert pre == (2015, 2016) and post == (2017, 2018)
+
+    def test_truncated_at_the_data_edge(self):
+        pre, post = windows(2024, available=range(2013, 2025), pre=3, post=3)
+        assert pre == (2021, 2022, 2023) and post == (2024,)
+        assert windows(2013, available=range(2013, 2025)) == ((), (2013, 2014))
+
+
+class TestLogChange:
+    def test_ratio_of_means(self):
+        s = {2015: 100.0, 2016: 100.0, 2017: 200.0, 2018: 200.0}
+        assert log_change(s, (2015, 2016), (2017, 2018)) == pytest.approx(math.log(2))
+
+    def test_missing_side_is_none(self):
+        assert log_change({2017: 1.0}, (2015, 2016), (2017,)) is None
+        assert log_change({2015: 1.0}, (2015,), (2017,)) is None
+        assert (
+            log_change({2015: 0.0, 2017: 1.0}, (2015,), (2017,)) is None
+        )  # zero pre → undefined
+
+
+class TestDidEffect:
+    def _controls(self, n=20, growth=0.1):
+        return {
+            f"C{i:04d}": {
+                2015: 100.0,
+                2016: 100.0,
+                2017: 100.0 * (1 + growth + 0.01 * (i % 5)),
+                2018: 100.0 * (1 + growth + 0.01 * (i % 5)),
+            }
+            for i in range(n)
+        }
+
+    def test_effect_is_treated_change_minus_control_change(self):
+        treated = {2015: 100.0, 2016: 100.0, 2017: 300.0, 2018: 300.0}
+        e = did_effect(
+            treated,
+            self._controls(),
+            2017,
+            code="X",
+            kind="becomes-payable",
+            outcome="services",
+        )
+        assert isinstance(e, Effect)
+        assert e.pre_years == (2015, 2016) and e.post_years == (2017, 2018)
+        assert e.treated_change == pytest.approx(math.log(3))
+        assert 0.09 < e.control_change < 0.14
+        assert e.effect == pytest.approx(e.treated_change - e.control_change)
+        assert e.ci_low is not None and e.ci_low < e.effect < e.ci_high
+        assert e.n_controls == 20
+        assert e.effect_pct == pytest.approx(math.exp(e.effect) - 1)
+
+    def test_ci_is_deterministic_and_shrinks_with_more_controls(self):
+        treated = {2015: 100.0, 2016: 100.0, 2017: 150.0, 2018: 150.0}
+        a = did_effect(
+            treated, self._controls(10), 2017, code="X", kind="k", outcome="o"
+        )
+        b = did_effect(
+            treated, self._controls(10), 2017, code="X", kind="k", outcome="o"
+        )
+        c = did_effect(
+            treated, self._controls(200), 2017, code="X", kind="k", outcome="o"
+        )
+        assert (a.ci_low, a.ci_high) == (b.ci_low, b.ci_high)
+        assert (c.ci_high - c.ci_low) < (a.ci_high - a.ci_low)
+
+    def test_too_few_controls_gives_no_ci(self):
+        treated = {2015: 1.0, 2016: 1.0, 2017: 2.0}
+        e = did_effect(
+            treated,
+            {"C1": {2015: 1.0, 2017: 1.1}},
+            2017,
+            code="X",
+            kind="k",
+            outcome="o",
+        )
+        assert e.n_controls == 1 and e.ci_low is None and e.effect is not None
+        e0 = did_effect(treated, {}, 2017, code="X", kind="k", outcome="o")
+        assert e0.n_controls == 0 and e0.control_change is None and e0.effect is None
+
+    def test_treated_without_post_data(self):
+        e = did_effect(
+            {2015: 1.0, 2016: 1.0},
+            self._controls(),
+            2017,
+            code="X",
+            kind="k",
+            outcome="o",
+        )
+        assert e.treated_change is None and e.effect is None and e.n_controls == 20
+
+
+class TestPrePost:
+    def test_means_and_change(self):
+        p = pre_post({2015: 10.0, 2016: 12.0, 2017: 20.0, 2018: 24.0}, 2017)
+        assert isinstance(p, PrePost)
+        assert p.pre_mean == 11.0 and p.post_mean == 22.0
+        assert p.change == pytest.approx(math.log(2))
+
+    def test_empty(self):
+        p = pre_post({}, 2017)
+        assert p.pre_mean is None and p.change is None
+
+
+RVU_COLS = "hcpcs VARCHAR, mod VARCHAR, description VARCHAR, status_code VARCHAR, non_fac_total DOUBLE, year INTEGER"
+
+
+def _util(code, year, services):
+    return UtilizationRow(
+        year,
+        code,
+        "d",
+        "N",
+        "O",
+        10,
+        100,
+        float(services),
+        float(services),
+        1.0,
+        2.0,
+        3.0,
+        4.0,
+        "ds",
+        "K",
+    )
+
+
+def _exp(code, year, kind, status="A", band=1, item_key="ITEM", p_id=5):
+    return ExposureRow(
+        code,
+        "FAM",
+        dt.date(year, 1, 1),
+        year,
+        kind,
+        status,
+        "",
+        band,
+        1.0,
+        0.0,
+        item_key,
+        p_id,
+        9,
+        "created",
+        False,
+        "",
+    )
+
+
+@pytest.fixture
+def con():
+    c = duckdb.connect(":memory:")
+    ensure_tables(c)
+    c.execute(f"CREATE TABLE pfs.rvu ({RVU_COLS})")
+    rows = []
+    # treated code T appears payable 2017; controls C1..C6 payable across the span in the same band
+    for y in range(2015, 2021):
+        for i in range(1, 7):
+            rows.append((f"C000{i}", None, "c", "A", 1.0, y))
+        if y >= 2017:
+            rows.append(("T0001", None, "t", "A", 1.0, y))
+    c.executemany("INSERT INTO pfs.rvu VALUES (?,?,?,?,?,?)", rows)
+    write_exposures(c, [_exp("T0001", 2017, "becomes-payable")])
+    util = []
+    for y in range(2015, 2021):
+        for i in range(1, 7):
+            util.append(_util(f"C000{i}", y, 100 * (1 + 0.05 * (y - 2015))))
+        if y >= 2017:
+            util.append(_util("T0001", y, 500 + 100 * (y - 2017)))
+    for y in range(2015, 2021):
+        write_util(c, y, [r for r in util if r.year == y])
+    yield c
+    c.close()
+
+
+class TestUtilizationPanel:
+    def test_panel_uses_control_set_and_carries_anchor(self, con):
+        panel = utilization_panel(con, "T0001", 2017, "becomes-payable")
+        assert panel.effect.code == "T0001" and panel.effect.outcome == "n_services"
+        assert panel.effect.n_controls == 6
+        assert (
+            panel.exposure is not None
+            and panel.exposure.item_key == "ITEM"
+            and panel.exposure.p_id == 5
+        )
+        # treated has no pre-period utilization (it did not exist) → change undefined, stated not faked
+        assert panel.effect.treated_change is None and panel.effect.scope == "code"
+        assert set(panel.controls) == {f"C000{i}" for i in range(1, 7)}
+        assert panel.treated == {2017: 500.0, 2018: 600.0, 2019: 700.0, 2020: 800.0}
+
+    def test_panel_with_pre_period(self, con):
+        write_exposures(
+            con,
+            [_exp("T0001", 2017, "becomes-payable"), _exp("C0001", 2018, "revalued")],
+        )
+        panel = utilization_panel(con, "C0001", 2018, "revalued")
+        assert panel.effect.treated_change is not None
+        assert panel.effect.n_controls == 5  # the other controls, not itself
+        assert panel.effect.effect == pytest.approx(
+            0.0, abs=1e-9
+        )  # it grows like its controls
+
+    def test_new_code_falls_back_to_the_family_series(self, con):
+        # T0001 has no pre-period of its own; its family (FAM) also holds P0001,
+        # billed before 2017 — the family total stands in as the treated series.
+        con.executemany(
+            "INSERT INTO pfs.code_family VALUES (?,?,?,?,?,?,?,?,?)",
+            [
+                ("FAM", "Fam", "T0001", "successor", 2017, None, "", 0, ""),
+                ("FAM", "Fam", "P0001", "predecessor", 2015, 2016, "", 0, ""),
+            ],
+        )
+        for y in (2015, 2016):
+            write_util(
+                con,
+                y,
+                [
+                    _util(f"C000{i}", y, 100 * (1 + 0.05 * (y - 2015)))
+                    for i in range(1, 7)
+                ]
+                + [_util("P0001", y, 400.0)],
+            )
+        panel = utilization_panel(con, "T0001", 2017, "becomes-payable")
+        assert panel.effect.scope == "family"
+        assert panel.treated[2015] == 400.0 and panel.treated[2017] == 500.0
+        assert (
+            panel.effect.treated_change is not None and panel.effect.effect is not None
+        )
+        assert panel.effect.n_controls == 6
+
+    def test_unknown_exposure_still_returns_a_panel(self, con):
+        panel = utilization_panel(con, "ZZZZZ", 2018, "ends")
+        assert (
+            panel.exposure is None
+            and panel.effect.n_controls == 0
+            and panel.treated == {}
+        )