"""month-of-year-effects-false-discovery-rate.py

Month-of-year (seasonality) tests on five FRED series, corrected for multiple testing four ways, with a time-split
replication and two hypotheses fixed in advance. Prism Data Lab, 2026-10-09. Nothing here is investment advice.

What it does
  1. Reads five FRED series (no API key) and keeps the last observation dated in each calendar month.
  2. Turns each into monthly changes, March 1971 at the earliest and September 2026 at the latest: percent log
     changes for the Nasdaq Composite, WTI crude oil, retail gasoline and Henry Hub natural gas, and basis points
     for the 10-year Treasury yield.
  3. For every series and calendar month, compares that month's mean change with the mean of the other eleven
     months: a two-sided permutation p from 20,000 shuffles of the month labels (fixed seeds), and a Welch t with its
     own p value (Student t, Welch-Satterthwaite degrees of freedom).
  4. Corrects the 60 permutation p values (5 series x 12 months) with Bonferroni, Holm, Benjamini-Hochberg and
     Benjamini-Yekutieli, all at 0.05, and repeats the corrections on the Welch p values as a robustness check (added
     after the analysis plan was fixed; the article says so).
  5. Runs one omnibus test per series (are the twelve month means equal?) from the same shuffles.
  6. Splits each series at its middle month, screens the first half the same way, and checks every first-half
     discovery on the second half, one-sided in the direction it was found.
  7. Tests two hypotheses fixed in advance on the Nasdaq: January beats the other months, and November-April beats
     May-October.

Inputs (one of)
  --download        fetch the five series from https://fred.stlouisfed.org/graph/fredgraph.csv?id=<ID>
  --raw-dir DIR     read DIR/<ID>.csv, files saved from that endpoint earlier
Outputs
  The tables, printed. With --csv, three files in --out (default: the datasets folder beside this script's folder):
    month-of-year-effects-false-discovery-rate.csv          one row per series, sample and calendar month
    month-of-year-effects-false-discovery-rate-series.csv   one row per series and sample
    month-of-year-effects-false-discovery-rate-monthly.csv  the monthly changes of the four public-domain series
  FRED marks NASDAQCOM "Copyrighted: Pre-Approval Required", so this script never writes its monthly values to
  disk; the files carry only its per-month statistics. With --figures, the chart
  content/images/month-of-year-effects-false-discovery-rate/month-effects-by-series.png (or --fig-out).

Run
  python month-of-year-effects-false-discovery-rate.py --download
  python month-of-year-effects-false-discovery-rate.py --raw-dir fred --csv --figures
Standard library only for every number, Python 3.9 or later; matplotlib (3.10) only for --figures. A run draws 20,000
shuffles for each series in each of three samples.
"""
from __future__ import annotations

import argparse
import csv
import io
import math
import os
import random
import sys
import time
import urllib.request

# (FRED id, label, unit of the monthly change, transform, public domain on FRED)
SERIES = [
    ("NASDAQCOM", "Nasdaq Composite", "% (log)", "log", False),
    ("DGS10", "10-year Treasury yield", "bp", "diff", True),
    ("DCOILWTICO", "WTI crude oil", "% (log)", "log", True),
    ("GASREGW", "Retail gasoline", "% (log)", "log", True),
    ("DHHNGSP", "Henry Hub natural gas", "% (log)", "log", True),
]
LAST_MONTH = "2026-09"          # October 2026 was incomplete when the data were pulled
DRAWS = 20000
SEED = 20261009
ALPHA = 0.05
EPS = 1e-9                      # tolerance for ties between a shuffled statistic and the observed one
MONTHS = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
WINTER = (10, 11, 0, 1, 2, 3)   # November to April, as 0-based month indexes
USER_AGENT = "month-of-year-study/1.0 (put your own contact here)"
SLUG = "month-of-year-effects-false-discovery-rate"
STAGES = ("full", "first", "second")


# ----------------------------------------------------------------- data ---
def parse_fredgraph(text: str) -> list[tuple[str, str]]:
    rows = list(csv.reader(io.StringIO(text)))
    out = [(r[0], r[1].strip()) for r in rows[1:] if r]
    if [d for d, _ in out] != sorted(d for d, _ in out):
        raise ValueError("FRED rows are not in date order")
    return out


def load(series_id: str, raw_dir: str | None) -> list[tuple[str, str]]:
    if raw_dir:
        with open(os.path.join(raw_dir, series_id + ".csv"), encoding="utf-8") as f:
            return parse_fredgraph(f.read())
    url = "https://fred.stlouisfed.org/graph/fredgraph.csv?id=" + series_id
    # FRED's CSV endpoint stalls when a request has no Accept header (urllib sends none by default).
    req = urllib.request.Request(url, headers={"User-Agent": USER_AGENT, "Accept": "*/*"})
    with urllib.request.urlopen(req, timeout=90) as r:
        text = r.read().decode("utf-8")
    time.sleep(1.5)                 # FRED's robots.txt asks for a one-second crawl delay
    return parse_fredgraph(text)


def month_ends(rows) -> dict[str, tuple[str, float, str]]:
    """The last observation dated in each calendar month: {YYYY-MM: (date, value, value as FRED wrote it)}. Missing
    values are skipped."""
    last: dict[str, tuple[str, float, str]] = {}
    for d, v in rows:
        if v in ("", "."):
            continue
        last[d[:7]] = (d, float(v), v)
    return last


def next_month(ym: str) -> str:
    y, m = int(ym[:4]), int(ym[5:])
    return f"{y + (m == 12)}-{m % 12 + 1:02d}"


def monthly_changes(last: dict, transform: str) -> list[tuple[str, float]]:
    out = []
    months = sorted(last)
    for a, b in zip(months, months[1:]):
        if b > LAST_MONTH:
            break
        if next_month(a) != b:      # a calendar month with no observation breaks the chain
            continue
        va, vb = last[a][1], last[b][1]
        if transform == "log":
            if va <= 0 or vb <= 0:
                raise ValueError(f"non-positive month-end price at {a} or {b}")
            out.append((b, 100 * math.log(vb / va)))
        else:
            out.append((b, 100 * (vb - va)))
    return out


# ----------------------------------------------------------- statistics ---
def mean(xs):
    return sum(xs) / len(xs)


def sd(xs):
    m = mean(xs)
    return math.sqrt(sum((x - m) ** 2 for x in xs) / (len(xs) - 1))


def welch_t(xs, ys):
    return (mean(xs) - mean(ys)) / math.sqrt(sd(xs) ** 2 / len(xs) + sd(ys) ** 2 / len(ys))


def welch_df(xs, ys):
    a, b = sd(xs) ** 2 / len(xs), sd(ys) ** 2 / len(ys)
    return (a + b) ** 2 / (a * a / (len(xs) - 1) + b * b / (len(ys) - 1))


def _betacf(a, b, x):
    """Continued fraction for the regularized incomplete beta function (modified Lentz)."""
    tiny = 1e-300
    qab, qap, qam = a + b, a + 1, a - 1
    c, d = 1.0, 1 - qab * x / qap
    d = 1 / (d if abs(d) > tiny else tiny)
    h = d
    for m in range(1, 500):
        m2 = 2 * m
        for aa in (m * (b - m) * x / ((qam + m2) * (a + m2)), -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))):
            d = 1 + aa * d
            d = 1 / (d if abs(d) > tiny else tiny)
            c = 1 + aa / c
            c = c if abs(c) > tiny else tiny
            h *= d * c
        if abs(d * c - 1) < 1e-15:
            break
    return h


def t_p_two(t, df):
    """Two-sided p of Student's t with df degrees of freedom: I_x(df/2, 1/2) at x = df / (df + t^2)."""
    x, a, b = df / (df + t * t), df / 2, 0.5
    if x <= 0:
        return 0.0
    bt = math.exp(math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) + a * math.log(x) + b * math.log(1 - x))
    if x < (a + 1) / (a + b + 2):
        return bt * _betacf(a, b, x) / a
    return 1 - bt * _betacf(b, a, 1 - x) / b


def lag1(xs):
    m = mean(xs)
    return sum((a - m) * (b - m) for a, b in zip(xs[1:], xs)) / sum((x - m) ** 2 for x in xs)


def group_stats(sums, n, total, count):
    """Month-minus-rest differences, the Nov-Apr minus May-Oct difference, and sum(s_m^2 / n_m), whose ordering is
    that of the between-month sum of squares (the omnibus statistic), all from the twelve month sums."""
    d = [sums[m] / n[m] - (total - sums[m]) / (count - n[m]) for m in range(12)]
    sw, nw = sum(sums[m] for m in WINTER), sum(n[m] for m in WINTER)
    h = sw / nw - (total - sw) / (count - nw)
    b = sum(sums[m] * sums[m] / n[m] for m in range(12))
    return d, h, b


def permutation(values, month_idx, draws, seed):
    """Joint permutation test: shuffle the values, deal them back into month groups of the observed sizes, and count
    how often each shuffled statistic is at least as extreme as the observed one. p = (1 + count) / (draws + 1)."""
    n = [0] * 12
    for m in month_idx:
        n[m] += 1
    total, count = sum(values), len(values)
    obs = [0.0] * 12
    for x, m in zip(values, month_idx):
        obs[m] += x
    d0, h0, b0 = group_stats(obs, n, total, count)
    bounds, a = [], 0
    for m in range(12):
        bounds.append((a, a + n[m]))
        a += n[m]
    two, up, down, h_up, b_ge = [0] * 12, [0] * 12, [0] * 12, 0, 0
    rng = random.Random(seed)
    vals = list(values)
    for _ in range(draws):
        rng.shuffle(vals)
        d, h, b = group_stats([sum(vals[i:j]) for i, j in bounds], n, total, count)
        for m in range(12):
            if abs(d[m]) >= abs(d0[m]) - EPS:
                two[m] += 1
            if d[m] >= d0[m] - EPS:
                up[m] += 1
            if d[m] <= d0[m] + EPS:
                down[m] += 1
        if h >= h0 - EPS:
            h_up += 1
        if b >= b0 - EPS:
            b_ge += 1
    p = lambda c: (1 + c) / (draws + 1)  # noqa: E731
    return {"n": n, "d": d0, "h": h0, "p_two": [p(c) for c in two], "p_up": [p(c) for c in up],
            "p_down": [p(c) for c in down], "p_halloween_up": p(h_up), "p_omnibus": p(b_ge)}


def corrections(ps, alpha=ALPHA):
    """Index sets kept by each procedure, and the adjusted p values (BH, BY, Holm, Bonferroni)."""
    m = len(ps)
    order = sorted(range(m), key=lambda i: (ps[i], i))  # smallest p first
    hm = sum(1 / j for j in range(1, m + 1))             # 4.6799 for m = 60
    kept = {"raw": {i for i in range(m) if ps[i] < alpha},
            "bonferroni": {i for i in range(m) if ps[i] <= alpha / m}}
    holm = set()
    for rank, i in enumerate(order, 1):                  # step down: stop at the first miss
        if ps[i] > alpha / (m - rank + 1):
            break
        holm.add(i)
    kept["holm"] = holm
    for name, q in (("bh", alpha), ("by", alpha / hm)):
        k = 0
        for rank, i in enumerate(order, 1):              # step up: keep up to the last pass
            if ps[i] <= rank * q / m:
                k = rank
        kept[name] = set(order[:k])
    adj = {"bh": [0.0] * m, "by": [0.0] * m, "holm": [0.0] * m, "bonferroni": [min(1.0, m * p) for p in ps]}
    run_bh, run_by = 1.0, 1.0
    for rank in range(m, 0, -1):                      # step-up: running minimum from the largest p down
        i = order[rank - 1]
        run_bh = min(run_bh, m * ps[i] / rank)
        run_by = min(run_by, m * hm * ps[i] / rank)
        adj["bh"][i], adj["by"][i] = min(1.0, run_bh), min(1.0, run_by)
    run = 0.0
    for rank, i in enumerate(order, 1):               # step-down: running maximum from the smallest p up
        run = max(run, (m - rank + 1) * ps[i])
        adj["holm"][i] = min(1.0, run)
    return kept, adj, hm


# --------------------------------------------------------------- figure ---
def figure(res, cells, kept, path):
    """Five small multiples, one per series, each on its own scale (the units differ). Bar colour is the strength of
    evidence on a one-hue ramp: p >= 0.05, p < 0.05 but not kept, kept by Benjamini-Hochberg. Needs matplotlib."""
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from matplotlib.patches import Patch
    from matplotlib.ticker import FuncFormatter, MaxNLocator

    paper, ink, muted, grid, axis = "#ffffff", "#101820", "#56616b", "#e3e6e8", "#9aa3ab"
    ramp = {"none": "#8fbcc2", "raw": "#4a8a94", "bh": "#0b5563"}  # one-hue ramp, validated: L 0.765 / 0.595 / 0.415
    minus = lambda s: s.replace("-", "−")  # noqa: E731
    plt.rcParams.update({"font.family": "DejaVu Sans", "font.size": 10, "text.color": ink, "xtick.color": muted,
                         "ytick.color": muted})
    fig, axes = plt.subplots(len(SERIES), 1, figsize=(8.8, 10.4), dpi=150, facecolor=paper)
    fig.subplots_adjust(left=0.075, right=0.985, top=0.885, bottom=0.03, hspace=0.75)
    fig.text(0.075, 0.966, "Mean change in each calendar month minus the mean of the other eleven", fontsize=12.5,
             ha="left", va="center")
    handles = [Patch(color=ramp["bh"], label="kept by Benjamini-Hochberg, q = 0.05"),
               Patch(color=ramp["raw"], label="p < 0.05 unadjusted, not kept"),
               Patch(color=ramp["none"], label="p ≥ 0.05")]
    fig.legend(handles=handles, loc="center left", bbox_to_anchor=(0.068, 0.93), ncol=3, frameon=False, fontsize=9.5,
               handlelength=1.1, handleheight=0.9, columnspacing=1.6)
    unit_text = {"% (log)": "percent (log change)", "bp": "basis points"}
    for ax, (sid, label, unit, *_) in zip(axes, SERIES):
        r = res[(sid, "full")]
        vals = r["d"]
        lev = []
        for j in range(12):
            i = cells.index((sid, j))
            lev.append("bh" if i in kept["bh"] else "raw" if i in kept["raw"] else "none")
        ax.set_facecolor(paper)
        ax.bar(range(12), vals, width=0.58, color=[ramp[x] for x in lev], linewidth=0, zorder=3)
        ax.axhline(0, color=axis, linewidth=1.0, zorder=2)
        ax.yaxis.grid(True, color=grid, linewidth=0.8, zorder=0)
        for s in ax.spines.values():
            s.set_visible(False)
        lo, hi = min(min(vals), 0), max(max(vals), 0)
        pad = 0.24 * (hi - lo)
        ax.set_ylim(lo - pad, hi + pad)
        ax.yaxis.set_major_locator(MaxNLocator(nbins=4, steps=[1, 2, 2.5, 5, 10]))
        ax.yaxis.set_major_formatter(FuncFormatter(lambda v, _: minus(f"{v:+g}") if v else "0"))
        ax.set_xticks(range(12), [m for m in MONTHS])
        ax.tick_params(length=0, labelsize=9, pad=3)
        ax.set_xlim(-0.6, 11.6)
        ax.set_title(f"{label}, {unit_text[unit]}", loc="left", fontsize=10.5, pad=5)
        ax.set_title(f"{r['first']} to {r['last']}, {r['count']} months", loc="right", fontsize=9, color=muted, pad=5)
        for j, v in enumerate(vals):
            if lev[j] == "bh":
                ax.annotate(minus(f"{v:+.1f}"), (j, v), textcoords="offset points", xytext=(0, 3 if v > 0 else -3),
                            ha="center", va="bottom" if v > 0 else "top", fontsize=9, color=ink, zorder=4)
    os.makedirs(os.path.dirname(path), exist_ok=True)
    fig.savefig(path, facecolor=paper)
    plt.close(fig)
    print(f"wrote {path} ({os.path.getsize(path) / 1024:.0f} KB)")


# ----------------------------------------------------------------- run ---
def fmt_p(p):
    return f"{p:.5f}"


def main():
    ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
    src = ap.add_mutually_exclusive_group(required=True)
    src.add_argument("--download", action="store_true", help="fetch the five series from FRED")
    src.add_argument("--raw-dir", help="read <DIR>/<ID>.csv saved from FRED's fredgraph.csv endpoint")
    ap.add_argument("--csv", action="store_true", help="write the three CSV files")
    ap.add_argument("--out", default=os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "datasets"))
    ap.add_argument("--draws", type=int, default=DRAWS)
    ap.add_argument("--figures", action="store_true", help="draw the chart (needs matplotlib)")
    ap.add_argument("--fig-out", default=os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
                                                      "content", "images", SLUG, "month-effects-by-series.png"))
    a = ap.parse_args()
    sys.stdout.reconfigure(encoding="utf-8")
    started = time.time()

    data = {}
    for k, (sid, label, unit, transform, public) in enumerate(SERIES):
        rows = load(sid, None if a.download else a.raw_dir)
        last = month_ends(rows)
        ch = monthly_changes(last, transform)
        obs = [(d, v) for d, v in rows if v not in ("", ".")]
        data[sid] = {"rows": rows, "last": last, "changes": ch, "k": k}
        print(f"{sid:<11} {len(obs):>6,} observations {obs[0][0]} to {obs[-1][0]}; "
              f"{len(ch)} monthly changes {ch[0][0]} to {ch[-1][0]}")

    # ------------------------------------------------ permutation tests ---
    res = {}
    for stage_no, stage in enumerate(STAGES):
        for sid, label, unit, transform, public in SERIES:
            ch = data[sid]["changes"]
            half = len(ch) // 2
            part = {"full": ch, "first": ch[:half], "second": ch[half:]}[stage]
            vals = [v for _, v in part]
            idx = [int(d[5:]) - 1 for d, _ in part]
            r = permutation(vals, idx, a.draws, SEED + 100 * stage_no + data[sid]["k"])
            by_m = [[v for v, m in zip(vals, idx) if m == j] for j in range(12)]
            rest = [[v for v, m in zip(vals, idx) if m != j] for j in range(12)]
            t = [welch_t(by_m[j], rest[j]) for j in range(12)]
            dfs = [welch_df(by_m[j], rest[j]) for j in range(12)]
            r.update(first=part[0][0], last=part[-1][0], count=len(vals), vals=vals, by_m=by_m,
                     mean_m=[mean(g) for g in by_m], mean_o=[mean(g) for g in rest], sd_m=[sd(g) for g in by_m],
                     t=t, df=dfs, p_welch=[t_p_two(t[j], dfs[j]) for j in range(12)],
                     lag1=lag1(vals), sd_all=sd(vals))
            res[(sid, stage)] = r

    cells = [(sid, j) for sid, *_ in SERIES for j in range(12)]
    fam = {}
    for stage in STAGES:
        ps = [res[(sid, stage)]["p_two"][j] for sid, j in cells]
        fam[stage] = corrections(ps)

    labels = {sid: label for sid, label, *_ in SERIES}
    units = {sid: unit for sid, _, unit, *_ in SERIES}

    # ------------------------------------------------------------ tables ---
    print("\nSamples (monthly changes; the split is each series' middle month)")
    for sid, *_ in SERIES:
        f, s1, s2 = res[(sid, "full")], res[(sid, "first")], res[(sid, "second")]
        print(f"  {labels[sid]:<23} {f['first']} to {f['last']}  n={f['count']}  first half {s1['first']} to "
              f"{s1['last']} ({s1['count']}), second {s2['first']} to {s2['last']} ({s2['count']}); "
              f"sd {f['sd_all']:.2f} {units[sid]}, lag-1 autocorrelation {f['lag1']:+.3f}")

    kept, adj, hm = fam["full"]
    m = len(cells)
    print(f"\nFull sample: {m} tests. Unadjusted p < 0.05: {len(kept['raw'])}; Bonferroni (p <= {ALPHA / m:.6f}): "
          f"{len(kept['bonferroni'])}; Holm: {len(kept['holm'])}; Benjamini-Hochberg q=0.05: {len(kept['bh'])}; "
          f"Benjamini-Yekutieli q=0.05 (harmonic sum {hm:.4f}): {len(kept['by'])}")
    exp_null = m * ALPHA
    print(f"  expected unadjusted discoveries if every null were true: {exp_null:.1f}")

    print("\nEvery test with unadjusted p < 0.05, smallest p first")
    print("  series                  month     n   month mean  other mean    diff        t   perm p   BH adj  Holm adj  kept by")
    for i in sorted(kept["raw"], key=lambda i: (res[(cells[i][0], 'full')]['p_two'][cells[i][1]], i)):
        sid, j = cells[i]
        r = res[(sid, "full")]
        by = [name for name in ("bonferroni", "holm", "bh", "by") if i in kept[name]]
        print(f"  {labels[sid]:<23} {MONTHS[j]:<5} {r['n'][j]:>5} {r['mean_m'][j]:>+11.2f} {r['mean_o'][j]:>+11.2f} "
              f"{r['d'][j]:>+8.2f} {r['t'][j]:>+8.2f}  {fmt_p(r['p_two'][j])}  {adj['bh'][i]:.4f}  {adj['holm'][i]:.4f}   "
              f"{', '.join(by) if by else 'none'}")

    ps_full = [res[(sid, "full")]["p_two"][j] for sid, j in cells]
    order = sorted(range(m), key=lambda i: (ps_full[i], i))
    print(f"\nBenjamini-Hochberg step-up, q = {ALPHA}: compare the i-th smallest p with i x {ALPHA} / {m}")
    k_bh = max([rank for rank, i in enumerate(order, 1) if ps_full[i] <= rank * ALPHA / m], default=0)
    for rank, i in enumerate(order[:10], 1):
        sid, j = cells[i]
        thr = rank * ALPHA / m
        print(f"  {rank:>2}  {labels[sid]:<23} {MONTHS[j]}  p {ps_full[i]:.5f}  threshold {thr:.5f}  "
              f"{'p <= threshold' if ps_full[i] <= thr else 'p > threshold'}")
    print(f"  largest i with p(i) <= i x q / m: {k_bh}, so the {k_bh} smallest are kept")
    for p_show in (ps_full[order[2]],):
        print(f"  Monte Carlo standard error of a permutation p of {p_show:.5f} from {a.draws:,} shuffles: "
              f"{math.sqrt(p_show * (1 - p_show) / a.draws):.5f}")

    pw = [res[(sid, "full")]["p_welch"][j] for sid, j in cells]
    kw, adjw, _ = corrections(pw)
    print("\nRobustness (added after the plan): the same corrections on Welch t p values, which use each month's own "
          "variance")
    print(f"  unadjusted p < 0.05: {len(kw['raw'])}  Bonferroni: {len(kw['bonferroni'])}  Holm: {len(kw['holm'])}  "
          f"Benjamini-Hochberg: {len(kw['bh'])}  Benjamini-Yekutieli: {len(kw['by'])}")
    for name in ("raw", "bonferroni", "holm", "bh", "by"):
        same = "the same cells as the permutation p values" if kw[name] == kept[name] else "different cells"
        names = ", ".join(f"{labels[cells[i][0]]} {MONTHS[cells[i][1]]}" for i in sorted(kw[name], key=lambda i: (pw[i], i)))
        print(f"  {name:<10} {len(kw[name]):>2}: {names or 'none'} ({same})")
    for sid, j in (("DHHNGSP", 1), ("GASREGW", 2)):
        r = res[(sid, "full")]
        print(f"  {labels[sid]} {MONTHS[j]}: Welch t {r['t'][j]:+.2f} on {r['df'][j]:.1f} df, p {r['p_welch'][j]:.5f}; "
              f"permutation p {r['p_two'][j]:.5f}; month sd {r['sd_m'][j]:.1f}")
    g = res[("DHHNGSP", "full")]["sd_m"]
    print(f"  Henry Hub natural gas, sd of the monthly change by month: " + ", ".join(f"{MONTHS[j]} {g[j]:.1f}" for j in range(12)))

    print("\nAll twelve months, full sample: mean change in the month minus the other months (two-sided permutation p)")
    for sid, *_ in SERIES:
        r = res[(sid, "full")]
        cellsfmt = "  ".join(f"{MONTHS[j]} {r['d'][j]:+.2f} ({r['p_two'][j]:.3f})" for j in range(12))
        print(f"  {labels[sid]:<23} {cellsfmt}")

    print("\nOne omnibus test per series (are the twelve month means equal?), full sample")
    om = [res[(sid, "full")]["p_omnibus"] for sid, *_ in SERIES]
    okept, oadj, _ = corrections(om)
    for k, (sid, *_) in enumerate(SERIES):
        print(f"  {labels[sid]:<23} permutation p {fmt_p(om[k])}  BH-adjusted {oadj['bh'][k]:.4f}  "
              f"{'kept' if k in okept['bh'] else 'not kept'} at q=0.05")

    # ------------------------------------------------------ replication ---
    k1, adj1, _ = fam["first"]
    print("\nFirst half screened, second half checked (one-sided p < 0.05 in the direction found)")
    rep_rows = {}
    for name in ("raw", "bonferroni", "holm", "bh", "by"):
        found = sorted(k1[name], key=lambda i: (res[(cells[i][0], 'first')]['p_two'][cells[i][1]], i))
        ok = []
        for i in found:
            sid, j = cells[i]
            up = res[(sid, "first")]["d"][j] > 0
            p2 = res[(sid, "second")]["p_up" if up else "p_down"][j]
            ok.append((i, p2, p2 < ALPHA, up))
        rep_rows[name] = ok
        print(f"  {name:<10} first-half discoveries {len(found):>2}  replicated {sum(1 for x in ok if x[2]):>2}")
    print("  first-half discoveries (unadjusted p < 0.05) and their second-half test")
    for i, p2, good, up in rep_rows["raw"]:
        sid, j = cells[i]
        r1, r2 = res[(sid, "first")], res[(sid, "second")]
        tags = [nm for nm in ("bonferroni", "holm", "bh", "by") if i in k1[nm]]
        print(f"    {labels[sid]:<23} {MONTHS[j]:<4} first {r1['d'][j]:+.2f} p {fmt_p(r1['p_two'][j])}  second "
              f"{r2['d'][j]:+.2f} one-sided p {fmt_p(p2)}  {'REPLICATED' if good else 'not replicated'}"
              f"{'  [' + ', '.join(tags) + ']' if tags else ''}")

    # ------------------------------------------------- fixed in advance ---
    print("\nTwo hypotheses fixed in advance, Nasdaq Composite (one-sided permutation p; Holm across the two at 0.05)")
    for stage in STAGES:
        r = res[("NASDAQCOM", stage)]
        pair = sorted([("January", r["p_up"][0]), ("November-April", r["p_halloween_up"])], key=lambda x: x[1])
        holm_pair = []
        for rank, (nm, p) in enumerate(pair, 1):        # Holm for two tests: the smaller p against 0.025, then 0.05
            if p > ALPHA / (2 - rank + 1):
                break
            holm_pair.append(nm)
        print(f"  {stage:<6} {r['first']} to {r['last']}: January {r['mean_m'][0]:+.2f} vs other months "
              f"{r['mean_o'][0]:+.2f}, diff {r['d'][0]:+.2f}, p {fmt_p(r['p_up'][0])}; Nov-Apr minus May-Oct "
              f"{r['h']:+.2f}, p {fmt_p(r['p_halloween_up'])}; Holm keeps {', '.join(holm_pair) or 'neither'}")

    # --------------------------------------------------------------- csv ---
    if a.csv:
        os.makedirs(a.out, exist_ok=True)
        with open(os.path.join(a.out, SLUG + ".csv"), "w", encoding="utf-8", newline="") as f:
            w = csv.writer(f, lineterminator="\n")
            w.writerow(["series", "label", "unit", "sample", "first_month", "last_month", "month", "month_name", "n",
                        "mean_month", "mean_other", "diff", "sd_month", "welch_t", "welch_df", "p_welch", "p_two", "p_up",
                        "p_down", "p_bh", "p_by", "p_holm", "p_bonferroni"])
            for stage in STAGES:
                kept_s, adj_s, _ = fam[stage]
                for i, (sid, j) in enumerate(cells):
                    r = res[(sid, stage)]
                    w.writerow([sid, labels[sid], units[sid], stage, r["first"], r["last"], j + 1, MONTHS[j], r["n"][j],
                                f"{r['mean_m'][j]:.6f}", f"{r['mean_o'][j]:.6f}", f"{r['d'][j]:.6f}", f"{r['sd_m'][j]:.6f}",
                                f"{r['t'][j]:.6f}", f"{r['df'][j]:.4f}", f"{r['p_welch'][j]:.10f}",
                                f"{r['p_two'][j]:.10f}", f"{r['p_up'][j]:.10f}", f"{r['p_down'][j]:.10f}",
                                f"{adj_s['bh'][i]:.10f}", f"{adj_s['by'][i]:.10f}", f"{adj_s['holm'][i]:.10f}",
                                f"{adj_s['bonferroni'][i]:.10f}"])
        with open(os.path.join(a.out, SLUG + "-series.csv"), "w", encoding="utf-8", newline="") as f:
            w = csv.writer(f, lineterminator="\n")
            w.writerow(["series", "label", "unit", "sample", "first_month", "last_month", "n", "mean", "sd", "lag1_autocorr",
                        "p_omnibus", "nov_apr_minus_may_oct", "p_nov_apr_up", "draws", "seed"])
            for stage_no, stage in enumerate(STAGES):
                for sid, *_ in SERIES:
                    r = res[(sid, stage)]
                    w.writerow([sid, labels[sid], units[sid], stage, r["first"], r["last"], r["count"],
                                f"{mean(r['vals']):.6f}", f"{r['sd_all']:.6f}", f"{r['lag1']:.6f}", f"{r['p_omnibus']:.10f}",
                                f"{r['h']:.6f}", f"{r['p_halloween_up']:.10f}", a.draws, SEED + 100 * stage_no + data[sid]["k"]])
        public = [s for s in SERIES if s[4]]
        chd = {s[0]: dict(data[s[0]]["changes"]) for s in public}
        months = sorted(set().union(*[set(chd[s[0]]) for s in public]))
        with open(os.path.join(a.out, SLUG + "-monthly.csv"), "w", encoding="utf-8", newline="") as f:
            w = csv.writer(f, lineterminator="\n")
            head = ["month"]
            for sid, *_ in public:
                head += [f"{sid}_date", f"{sid}_level", f"{sid}_change"]
            w.writerow(head)
            for mo in months:
                row = [mo]
                for sid, *_ in public:
                    if mo in chd[sid]:
                        dt, _, raw = data[sid]["last"][mo]
                        row += [dt, raw, f"{chd[sid][mo]:.6f}"]
                    else:
                        row += ["", "", ""]
                w.writerow(row)
        print(f"\nwrote {SLUG}.csv, {SLUG}-series.csv and {SLUG}-monthly.csv to {a.out}")
    if a.figures:
        figure(res, cells, fam["full"][0], a.fig_out)
    print(f"done in {time.time() - started:.0f} s")


if __name__ == "__main__":
    main()