"""Revision checks (editorial-board review, 2026-10-17) from stored event files.

Everything here is recomputed from outputs/estimates/events*/ (stage F/F2
outputs); no raw minute data are needed. Day-block bootstrap, B=500, same seed.

  A. Pre-event levels: real vs placebo breadth/intensity at t-5 (the "no
     pre-event response" falsification quantified; by decade as well).
  B. Base rates: placebo breadth, eligible receivers, neighbor / non-neighbor
     levels by arm and decade (components of the network differential).
  C. Co-arrivals: share of real events that share their minute with other
     events; network differential and market-wide ΔΔ by co-arrival class.
  D. Pre-level-stratified post comparison (mean-reversion check): real minus
     placebo breadth at t+h within deciles of the pre-event breadth level.
  E. Alternative clustering: ticker-cluster bootstrap SEs.
  F. Signed response (real minus placebo), network differential by shock sign.

Outputs: outputs/estimates/revision_checks.json,
         outputs/tables/table9_design_checks.{tex,csv},
         outputs/tables/table10_network_components.{tex,csv}
         (+ copies in paper/tables/).
"""

from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import polars as pl

ROOT = Path(__file__).resolve().parents[1]
EV = ROOT / "outputs" / "estimates" / "events"
EV2 = ROOT / "outputs" / "estimates" / "events2"
TAB = ROOT / "outputs" / "tables"
PTAB = ROOT / "paper" / "tables"
B = 500
SEED = 20260713
DEC = {"1992–1999": (1992, 1999), "2000–2009": (2000, 2009),
       "2010–2019": (2010, 2019), "2020–2026": (2020, 2026)}


def load(folder: Path, prefix: str) -> pl.DataFrame:
    """Concatenate monthly event files, keeping the real→placebo row order and
    labelling each placebo with the minute of the real event it matches."""
    frames = []
    for p in sorted(folder.glob(f"{prefix}_*.parquet")):
        if p.stat().st_size <= 500:
            continue
        d = pl.read_parquet(p).with_row_index("i")
        prev = d.select((pl.col("i") + 1).alias("i"), pl.col("kind").alias("pk"),
                        pl.col("ticker").alias("pt"), pl.col("z").alias("pz"),
                        pl.col("date").alias("ev_date"), pl.col("mod").alias("ev_mod"))
        d = d.join(prev, on="i", how="left")
        ok = (pl.col("kind") == "real") | (
            (pl.col("pk") == "real") & (pl.col("pt") == pl.col("ticker"))
            & (pl.col("pz") == pl.col("z")))
        d = d.filter(ok).with_columns(
            pl.when(pl.col("kind") == "real").then(pl.col("date")).otherwise(pl.col("ev_date")).alias("ev_date"),
            pl.when(pl.col("kind") == "real").then(pl.col("mod")).otherwise(pl.col("ev_mod")).alias("ev_mod"),
        ).drop("i", "pk", "pt", "pz")
        frames.append(d)
    df = pl.concat(frames, how="diagonal")
    return df.with_columns(pl.col("date").str.slice(0, 4).cast(pl.Int32).alias("year"))


def dd_boot(df: pl.DataFrame, col: str, rng, cluster: str = "date"):
    """Real-minus-placebo difference of means with cluster bootstrap SE."""
    g = (df.filter(pl.col(col).is_not_null())
         .group_by(cluster, "kind")
         .agg(pl.col(col).sum().alias("s"), pl.len().alias("n"))
         .pivot(on="kind", index=cluster, values=["s", "n"]))
    need = ("s_real", "n_real", "s_placebo", "n_placebo")
    if any(c not in g.columns for c in need) or g.height < 10:
        return np.nan, np.nan, 0
    g = g.fill_null(0).sort(cluster)  # deterministic cluster order → reproducible draws
    sr, nr = g["s_real"].to_numpy(), g["n_real"].to_numpy()
    sp, npb = g["s_placebo"].to_numpy(), g["n_placebo"].to_numpy()
    est = sr.sum() / max(nr.sum(), 1) - sp.sum() / max(npb.sum(), 1)
    idx = rng.integers(0, len(sr), size=(B, len(sr)))
    bs = (sr[idx].sum(1) / np.maximum(nr[idx].sum(1), 1)
          - sp[idx].sum(1) / np.maximum(npb[idx].sum(1), 1))
    return float(est), float(bs.std()), int(nr.sum())


def level_boot(df: pl.DataFrame, col: str, kind: str, rng):
    g = (df.filter((pl.col("kind") == kind) & pl.col(col).is_not_null())
         .group_by("date").agg(pl.col(col).sum().alias("s"), pl.len().alias("n")).sort("date"))
    s, n = g["s"].to_numpy(), g["n"].to_numpy()
    est = s.sum() / max(n.sum(), 1)
    idx = rng.integers(0, len(s), size=(B, len(s)))
    bs = s[idx].sum(1) / np.maximum(n[idx].sum(1), 1)
    return float(est), float(bs.std()), int(n.sum())


def stratified_boot(df: pl.DataFrame, post: str, pre: str, nbins: int, rng):
    """Real minus placebo mean of `post` within quantile bins of `pre`,
    weighted by the real-arm bin counts; day bootstrap."""
    d = df.filter(pl.col(post).is_not_null() & pl.col(pre).is_not_null())
    qs = [d[pre].quantile(q / nbins) for q in range(1, nbins)]
    d = d.with_columns(pl.col(pre).cut(qs, labels=[str(k) for k in range(nbins)]).alias("bin"))
    g = (d.group_by("date", "bin", "kind")
         .agg(pl.col(post).sum().alias("s"), pl.len().alias("n")))
    dates = sorted(g["date"].unique().to_list())
    di = {x: k for k, x in enumerate(dates)}
    S = {k: np.zeros((len(dates), nbins)) for k in ("real", "placebo")}
    N = {k: np.zeros((len(dates), nbins)) for k in ("real", "placebo")}
    for r in g.iter_rows(named=True):
        S[r["kind"]][di[r["date"]], int(r["bin"])] += r["s"]
        N[r["kind"]][di[r["date"]], int(r["bin"])] += r["n"]

    def stat(ix):
        sr, nr = S["real"][ix].sum(0), N["real"][ix].sum(0)
        sp, npb = S["placebo"][ix].sum(0), N["placebo"][ix].sum(0)
        w = nr / nr.sum()
        return float(np.sum(w * (sr / np.maximum(nr, 1) - sp / np.maximum(npb, 1))))

    est = stat(np.arange(len(dates)))
    idx = rng.integers(0, len(dates), size=(B, len(dates)))
    bs = np.array([stat(idx[b]) for b in range(B)])
    return est, float(bs.std()), int(N["real"].sum())


def main() -> None:
    rng = np.random.default_rng(SEED)
    e1 = load(EV, "events").with_columns(
        *[(pl.col(f"breadth_{h}") - pl.col("pre_breadth")).alias(f"db_{h}") for h in (5, 15, 30)],
        (pl.col("intensity_5") - pl.col("pre_intensity")).alias("di_5"))
    e2 = load(EV2, "events2").with_columns(
        *[(pl.col(f"nbh_{h}") - pl.col(f"non_{h}")).alias(f"nd_{h}") for h in (5, 15, 30)])
    real1 = e1.filter(pl.col("kind") == "real")
    print(f"events: {real1.height:,} real, {e1.height - real1.height:,} paired placebos")

    out: dict = {"params": {"bootstrap": B, "seed": SEED}}

    # ---------------- A. pre-event levels
    pre = {}
    for col, lab in [("pre_breadth", "breadth"), ("pre_intensity", "intensity")]:
        r = level_boot(e1, col, "real", rng)
        p = level_boot(e1, col, "placebo", rng)
        dlt = dd_boot(e1, col, rng)
        pre[lab] = {"real": r[0], "placebo": p[0], "diff": dlt[0], "se": dlt[1], "n": dlt[2]}
        print(f"pre {lab}: real {r[0]:.5f} placebo {p[0]:.5f} diff {dlt[0]*1e4:+.2f} ({dlt[1]*1e4:.2f}) x1e-4")
    pre["breadth_by_decade"] = []
    for dec, (a, b) in DEC.items():
        sub = e1.filter(pl.col("year").is_between(a, b))
        dlt = dd_boot(sub, "pre_breadth", rng)
        pre["breadth_by_decade"].append({"dec": dec, "diff": dlt[0], "se": dlt[1], "n": dlt[2]})
    out["pre_event"] = pre

    # ---------------- B. base rates and network components
    base = {
        "placebo_breadth_5": float(e1.filter(pl.col("kind") == "placebo")["breadth_5"].mean()),
        "placebo_breadth_30": float(e1.filter(pl.col("kind") == "placebo")["breadth_30"].mean()),
        "nelig_5_mean": float(real1["nelig_5"].mean()),
        "nelig_5_by_decade": {d: float(real1.filter(pl.col("year").is_between(a, b))["nelig_5"].mean())
                              for d, (a, b) in DEC.items()},
        "placebo_nbh_5": float(e2.filter(pl.col("kind") == "placebo")["nbh_5"].mean()),
        "placebo_non_5": float(e2.filter(pl.col("kind") == "placebo")["non_5"].mean()),
        "event_days": int(real1["date"].n_unique()),
        "match_rate": float((e1.height - real1.height) / real1.height),
        "match_rate_by_decade": {d: float(e1.filter(pl.col("year").is_between(a, b) & (pl.col("kind") == "placebo")).height
                                         / max(real1.filter(pl.col("year").is_between(a, b)).height, 1))
                                 for d, (a, b) in DEC.items()},
        "placebo_cum_15": float(e2.filter(pl.col("kind") == "placebo")["cum_15"].mean()),
        "events_per_stock_week": float(real1.height / (500 * real1["date"].n_unique() / 5)),
    }
    comp = []
    for scope, (a, b) in [("all", (1992, 2026))] + list(DEC.items()):
        sub = e2.filter(pl.col("year").is_between(a, b))
        row = {"scope": scope}
        for kind in ("real", "placebo"):
            for col in ("nbh_5", "non_5", "nd_5"):
                est, se, n = level_boot(sub, col, kind, rng)
                row[f"{kind}_{col}"] = est
                row[f"{kind}_{col}_se"] = se
        dd, se, n = dd_boot(sub, "nd_5", rng)
        row.update({"dd": dd, "dd_se": se, "n": n})
        comp.append(row)
        print(f"net components {scope:10} real nd {row['real_nd_5']*1e4:+.1f} placebo nd {row['placebo_nd_5']*1e4:+.1f} dd {dd*1e4:+.1f} ({se*1e4:.1f})")
    out["base_rates"] = base
    out["network_components"] = comp

    # ---------------- C. co-arrivals
    cnt = real1.group_by("date", "mod").len().rename({"len": "n_co"})
    e1c = e1.join(cnt, left_on=["ev_date", "ev_mod"], right_on=["date", "mod"], how="left")
    e2c = e2.join(cnt, left_on=["ev_date", "ev_mod"], right_on=["date", "mod"], how="left")
    dist = (real1.join(cnt, on=["date", "mod"], how="left")
            .with_columns(pl.when(pl.col("n_co") == 1).then(pl.lit("1"))
                          .when(pl.col("n_co") == 2).then(pl.lit("2"))
                          .when(pl.col("n_co") <= 5).then(pl.lit("3–5"))
                          .otherwise(pl.lit("6+")).alias("cls"))
            .group_by("cls").len().sort("cls"))
    tot = dist["len"].sum()
    coarr = {"share_by_class": {r["cls"]: r["len"] / tot for r in dist.iter_rows(named=True)},
             "mean_events_per_event_minute": float(real1.height / cnt.height),
             "by_class": []}
    classes = {"1": (1, 1), "2": (2, 2), "3–5": (3, 5), "6+": (6, 10**6)}
    for cls, (lo, hi) in classes.items():
        s1 = e1c.filter(pl.col("n_co").is_between(lo, hi))
        s2 = e2c.filter(pl.col("n_co").is_between(lo, hi))
        row = {"cls": cls, "share": coarr["share_by_class"].get(cls, 0.0)}
        for h in (5, 30):
            est, se, n = dd_boot(s1, f"db_{h}", rng)
            row[f"db{h}"], row[f"db{h}_se"], row["n"] = est, se, n
        for h in (5, 15, 30):
            est, se, n2 = dd_boot(s2, f"nd_{h}", rng)
            row[f"nd{h}"], row[f"nd{h}_se"] = est, se
        coarr["by_class"].append(row)
        print(f"co-arrival {cls:4} share {row['share']:.3f} ΔΔb5 {row['db5']*1e4:+.2f} ({row['db5_se']*1e4:.2f}) nd5 {row['nd5']*1e4:+.1f} ({row['nd5_se']*1e4:.1f})")
    out["coarrival"] = coarr

    # ---------------- D. pre-level-stratified post comparison
    strat = []
    for h in (5, 15, 30):
        raw = dd_boot(e1, f"breadth_{h}", rng)          # post level, real − placebo
        dd = dd_boot(e1, f"db_{h}", rng)                # ΔΔ (baseline estimand)
        st = stratified_boot(e1, f"breadth_{h}", "pre_breadth", 10, rng)
        strat.append({"h": h, "post_diff": raw[0], "post_se": raw[1], "dd": dd[0], "dd_se": dd[1],
                      "strat": st[0], "strat_se": st[1], "n": st[2]})
        print(f"h={h}: post diff {raw[0]*1e4:+.2f} ({raw[1]*1e4:.2f}); ΔΔ {dd[0]*1e4:+.2f} ({dd[1]*1e4:.2f}); pre-stratified {st[0]*1e4:+.2f} ({st[1]*1e4:.2f})")
    out["stratified"] = strat

    # ---------------- E. ticker-cluster bootstrap
    clus = {}
    for lab, df, col in [("db5", e1, "db_5"), ("db30", e1, "db_30"), ("nd5", e2, "nd_5")]:
        d_est, d_se, _ = dd_boot(df, col, rng, cluster="date")
        t_est, t_se, _ = dd_boot(df, col, rng, cluster="ticker")
        clus[lab] = {"est": d_est, "se_day": d_se, "se_ticker": t_se}
        print(f"{lab}: est {d_est*1e4:+.2f} day-SE {d_se*1e4:.2f} ticker-SE {t_se*1e4:.2f}")
    out["cluster"] = clus

    # ---------------- F. signed response and network differential by sign
    signed = []
    for h in (1, 5, 15, 30):
        est, se, n = dd_boot(e1, f"signed_{h}", rng)
        signed.append({"h": h, "diff": est, "se": se, "n": n})
        print(f"signed h={h}: {est*1e4:+.2f} ({se*1e4:.2f}) x1e-4")
    out["signed"] = signed
    bysign = []
    for sgn, lab in [(-1, "negative"), (1, "positive")]:
        est, se, n = dd_boot(e2.filter(pl.col("sign") == sgn), "nd_5", rng)
        bysign.append({"sign": lab, "dd": est, "se": se, "n": n})
        print(f"network diff, {lab} shocks: {est*1e4:+.1f} ({se*1e4:.1f})")
    out["network_by_sign"] = bysign

    (ROOT / "outputs" / "estimates" / "revision_checks.json").write_text(json.dumps(out, indent=1))
    (ROOT / "outputs" / "diagnostics" / "revision_checks.json").write_text(json.dumps(out, indent=1))

    # ---------------- Table 9: design diagnostics
    f4 = lambda v, s: f"{v*1e4:,.2f} ({s*1e4:,.2f})"
    rows9 = []
    L = ["% Auto-generated by scripts/17_revision_checks.py — do not edit numbers.",
         "\\begin{table}[htbp]\\centering\\begin{adjustbox}{max width=\\textwidth}\\begin{threeparttable}",
         "\\caption{Design diagnostics: pre-event levels, co-arrivals, mean reversion, and clustering}",
         "\\label{tab:design_checks}",
         "\\small\\begin{tabular}{lrrr}", "\\toprule"]
    L += ["\\multicolumn{4}{l}{\\emph{Panel A: Pre-event levels at $t-5$ ($\\times10^{-4}$)}} \\\\",
          " & Real & Placebo & Real $-$ placebo \\\\ \\midrule"]
    for lab, key in [("Breadth", "breadth"), ("Intensity", "intensity")]:
        r = pre[key]
        L.append(f"{lab} & {r['real']*1e4:,.0f} & {r['placebo']*1e4:,.0f} & {f4(r['diff'], r['se'])} \\\\")
        rows9.append({"panel": "A", "row": lab, "real": r["real"], "placebo": r["placebo"], "diff": r["diff"], "se": r["se"]})
    L += ["\\addlinespace", "\\multicolumn{4}{l}{\\emph{Panel B: By number of source events in the same minute}} \\\\",
          "Events in minute (share of events) & $\\Delta\\Delta$Breadth$_{5}$ & $\\Delta\\Delta$Breadth$_{30}$ & Network diff$_{5}$ \\\\ \\midrule"]
    for r in coarr["by_class"]:
        L.append(f"{r['cls']} ({r['share']*100:.1f}\\%) & {f4(r['db5'], r['db5_se'])} & {f4(r['db30'], r['db30_se'])} & {f4(r['nd5'], r['nd5_se'])} \\\\")
        rows9.append({"panel": "B", "row": r["cls"], "share": r["share"], "db5": r["db5"], "db5_se": r["db5_se"],
                      "db30": r["db30"], "db30_se": r["db30_se"], "nd5": r["nd5"], "nd5_se": r["nd5_se"]})
    L += ["\\addlinespace", "\\multicolumn{4}{l}{\\emph{Panel C: Market-wide breadth, alternative estimands ($\\times10^{-4}$)}} \\\\",
          "Horizon & Post level, real $-$ placebo & $\\Delta\\Delta$ (baseline) & Pre-level stratified \\\\ \\midrule"]
    for r in strat:
        L.append(f"$h={r['h']}$ & {f4(r['post_diff'], r['post_se'])} & {f4(r['dd'], r['dd_se'])} & {f4(r['strat'], r['strat_se'])} \\\\")
        rows9.append({"panel": "C", "row": r["h"], "post": r["post_diff"], "post_se": r["post_se"], "dd": r["dd"],
                      "dd_se": r["dd_se"], "strat": r["strat"], "strat_se": r["strat_se"]})
    L += ["\\addlinespace", "\\multicolumn{4}{l}{\\emph{Panel D: Bootstrap standard errors by clustering unit ($\\times10^{-4}$)}} \\\\",
          "Statistic & Estimate & Day clusters & Ticker clusters \\\\ \\midrule"]
    for lab, key in [("$\\Delta\\Delta$Breadth$_{5}$", "db5"), ("$\\Delta\\Delta$Breadth$_{30}$", "db30"), ("Network diff$_{5}$", "nd5")]:
        r = clus[key]
        L.append(f"{lab} & {r['est']*1e4:,.2f} & ({r['se_day']*1e4:,.2f}) & ({r['se_ticker']*1e4:,.2f}) \\\\")
        rows9.append({"panel": "D", "row": key, "est": r["est"], "se_day": r["se_day"], "se_ticker": r["se_ticker"]})
    L += ["\\addlinespace", "\\multicolumn{4}{l}{\\emph{Panel E: Signed response, real $-$ placebo ($\\times10^{-4}$)}} \\\\",
          "Horizon & $h=1$ & $h=5$ & $h=30$ \\\\ \\midrule"]
    sg = {r["h"]: r for r in signed}
    L.append("sign$(e_{i,t})\\times$ mean $z_{j,t+h}$ & " + " & ".join(f4(sg[h]["diff"], sg[h]["se"]) for h in (1, 5, 30)) + " \\\\")
    L += ["\\bottomrule\\end{tabular}",
          "\\begin{tablenotes}[flushleft]\\footnotesize",
          "\\item Notes: All statistics from the baseline event files (isolated $|z|\\ge6$ source",
          "shocks, lagged top-500 universe, 1992--2026; matched within-stock placebos).",
          "Day-block bootstrap standard errors (500 replications) in parentheses unless stated.",
          "Panel A: levels of the outcomes five minutes before the event minute. Panel B: events",
          "classified by the number of source events detected in the same minute (the isolation",
          "filters cap the share of simultaneous $|z|>4$ names at 2\\%, not the number of $|z|\\ge6$",
          "events); each placebo inherits the class of its real event. Panel C: `post level' is",
          "the real-minus-placebo difference in breadth at $t+h$ without pre-differencing; `pre-level",
          "stratified' reweights the same difference within deciles of the pre-event breadth level.",
          "Panel D: standard errors from resampling trading days versus resampling source tickers.",
          "Panel E: cross-sectional mean standardized return of receivers signed by the source shock.",
          "\\end{tablenotes}\\end{threeparttable}\\end{adjustbox}\\end{table}"]
    (TAB / "table9_design_checks.tex").write_text("\n".join(L))
    (PTAB / "table9_design_checks.tex").write_text("\n".join(L))
    pl.DataFrame(rows9, strict=False).write_csv(TAB / "table9_design_checks.csv")

    # ---------------- Table 10: network components by decade
    L = ["% Auto-generated by scripts/17_revision_checks.py — do not edit numbers.",
         "\\begin{table}[htbp]\\centering\\begin{adjustbox}{max width=\\textwidth}\\begin{threeparttable}",
         "\\caption{Components of the network differential at $h=5$, by decade}",
         "\\label{tab:network_components}",
         "\\small\\begin{tabular}{lrrrrrrr}", "\\toprule",
         " & \\multicolumn{3}{c}{Real events} & \\multicolumn{3}{c}{Placebo events} & \\\\",
         "\\cmidrule(lr){2-4}\\cmidrule(lr){5-7}",
         "Sample & Neighbors & Others & Diff. & Neighbors & Others & Diff. & Real $-$ placebo \\\\ \\midrule"]
    rows10 = []
    for r in comp:
        lab = "Full sample" if r["scope"] == "all" else r["scope"]
        L.append(f"{lab} & {r['real_nbh_5']*1e4:,.0f} & {r['real_non_5']*1e4:,.0f} & {r['real_nd_5']*1e4:,.1f} ({r['real_nd_5_se']*1e4:,.1f})"
                 f" & {r['placebo_nbh_5']*1e4:,.0f} & {r['placebo_non_5']*1e4:,.0f} & {r['placebo_nd_5']*1e4:,.1f} ({r['placebo_nd_5_se']*1e4:,.1f})"
                 f" & {r['dd']*1e4:,.1f} ({r['dd_se']*1e4:,.1f}) \\\\")
        rows10.append(r)
    L += ["\\bottomrule\\end{tabular}",
          "\\begin{tablenotes}[flushleft]\\footnotesize",
          "\\item Notes: Share of receivers with $|z|>2$ five minutes after the event minute, among",
          "the source's twenty lagged-correlation neighbors and among all other receivers, for real",
          "events and for their matched placebos ($\\times10^{-4}$). `Diff.' is neighbors minus",
          "others within an arm; the last column is the network differential of",
          "Table~\\ref{tab:network_mechanism} (real minus placebo). Day-block bootstrap SEs (500",
          "replications) in parentheses.",
          "\\end{tablenotes}\\end{threeparttable}\\end{adjustbox}\\end{table}"]
    (TAB / "table10_network_components.tex").write_text("\n".join(L))
    (PTAB / "table10_network_components.tex").write_text("\n".join(L))
    pl.DataFrame(rows10).write_csv(TAB / "table10_network_components.csv")
    print("revision checks complete")


if __name__ == "__main__":
    main()