"""fomc-event-study.py - an event study around FOMC announcement dates, built
entirely from free sources with no API key.

Two downloads, both public:

  Federal Reserve   the FOMC calendar and the historical-materials pages, from
                    which the ANNOUNCEMENT dates are recovered by the statement
                    URL each page links (monetaryYYYYMMDDa.htm and its older
                    spellings). This is the day the decision was published, not
                    the day the meeting started, which is the date an event
                    study needs.
  FRED              daily series through the public fredgraph CSV endpoint.
                    DGS2 and DGS10 are constant-maturity Treasury yields in
                    per cent, VIXCLS is the VIX close, NASDAQCOM is the Nasdaq
                    Composite close.

Yield and VIX changes are reported in basis points of level change. Nasdaq
changes are log returns in basis points, so 100 is one per cent.

For every series and every offset in the event window it reports the mean move,
its standard error, a t statistic, how often the move was positive, the mean
absolute move, and a permutation p value obtained by drawing the same number of
days at random from the same period 20,000 times with a fixed seed. The
permutation test is the one to read: daily returns are fat-tailed and the t
statistic assumes they are not.

It then prints the multiple-testing arithmetic, because this script runs one
test per series per offset and some of them will look significant by chance,
and a robustness block that re-runs the headline result on sub-periods, without
the crisis years, and with the largest moves removed.

Standard library only (urllib, csv, json, math, random, statistics, argparse).
Python 3.9 or newer. One run makes about 30 HTTP requests.

    python fomc-event-study.py
    python fomc-event-study.py --since 2010 --csv fomc-event-study.csv
"""
from __future__ import annotations

import argparse
import csv
import io
import json
import math
import random
import re
import statistics
import urllib.request

FED = "https://www.federalreserve.gov"
CALENDAR = f"{FED}/monetarypolicy/fomccalendars.htm"
HISTORICAL = [f"{FED}/monetarypolicy/fomchistorical{year}.htm"
              for year in range(2000, 2021)]

# Every spelling the Fed has used for a statement URL. All of them carry the
# announcement date, which is what is being extracted; the page each one sits
# on is only a means of finding it.
STATEMENT_PATTERNS = [
    r"/newsevents/pressreleases/monetary(\d{8})a\.htm",
    r"/newsevents/press/monetary/(\d{8})a\.htm",
    r"/boarddocs/press/monetary/\d{4}/(\d{8})/",
    r"/monetarypolicy/files/monetary(\d{8})a1\.pdf",
]

SERIES = {"DGS2": "level", "DGS10": "level", "VIXCLS": "level",
          "NASDAQCOM": "logret"}
OFFSETS = [-2, -1, 0, 1, 2]
HEADLINE = ("NASDAQCOM", 0)
DRAWS = 20000
SEED = 20260922


def get(url):
    request = urllib.request.Request(
        url, headers={"User-Agent": "prismdatalab-research/1.0", "Accept": "*/*"})
    with urllib.request.urlopen(request, timeout=90) as response:
        return response.read().decode("utf-8", "replace")


def announcement_dates(since):
    """Every FOMC announcement date the Fed's own pages link a statement for."""
    found = set()
    for page in [CALENDAR] + HISTORICAL:
        text = get(page)
        for pattern in STATEMENT_PATTERNS:
            found.update(re.findall(pattern, text))
    dates = sorted(f"{d[:4]}-{d[4:6]}-{d[6:]}" for d in found)
    return [d for d in dates if d[:4] >= since]


def fred(series_id):
    """A FRED daily series as {date: value}, missing values dropped.

    The fredgraph CSV endpoint stalls until timeout if the request carries no
    Accept header, so get() always sends one.
    """
    text = get(f"https://fred.stlouisfed.org/graph/fredgraph.csv?id={series_id}")
    rows = list(csv.reader(io.StringIO(text)))
    return {row[0]: float(row[1]) for row in rows[1:]
            if row and row[1] not in (".", "")}


def changes(values, kind):
    """Daily change in basis points, keyed by the later of the two days."""
    dates = sorted(values)
    out = {}
    for i in range(1, len(dates)):
        before, after = values[dates[i - 1]], values[dates[i]]
        if kind == "level":
            out[dates[i]] = (after - before) * 100
        elif before > 0:
            out[dates[i]] = 10000 * math.log(after / before)
    return out


def permutation_p(sample_mean, pool, n, rng):
    """Two-sided p for a mean of n draws from pool being this far from zero."""
    hits = 0
    for _ in range(DRAWS):
        drawn = sum(rng.choice(pool) for _ in range(n)) / n
        if abs(drawn) >= abs(sample_mean):
            hits += 1
    return hits / DRAWS


def describe(values):
    mean = statistics.mean(values)
    stderr = statistics.stdev(values) / math.sqrt(len(values))
    return {"n": len(values), "mean": mean, "se": stderr, "t": mean / stderr,
            "up": sum(1 for v in values if v > 0),
            "abs": statistics.mean(abs(v) for v in values),
            "median": statistics.median(values)}


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--since", default="2003",
                        help="first year to include (default 2003, the first "
                             "year the Fed's pages cover completely)")
    parser.add_argument("--csv", metavar="FILE", help="write the result table")
    args = parser.parse_args()

    events = announcement_dates(args.since)
    per_year = {}
    for date in events:
        per_year[date[:4]] = per_year.get(date[:4], 0) + 1
    print(f"FOMC announcement dates from {args.since}: {len(events)}")
    print("  " + "  ".join(f"{y}:{n}" for y, n in sorted(per_year.items())))
    print()

    rng = random.Random(SEED)
    rows = []
    for series_id, kind in SERIES.items():
        moves = changes(fred(series_id), kind)
        trading = [d for d in sorted(moves) if d >= f"{args.since}-01-01"]
        position = {d: i for i, d in enumerate(trading)}
        on_event = set(events)
        quiet = [moves[d] for d in trading if d not in on_event]
        pool = [moves[d] for d in trading]
        unit = "bp" if kind == "level" else "bp log"
        print(f"{series_id}  ({unit}, {len(trading)} trading days from "
              f"{trading[0]}, {len(quiet)} of them not announcement days)")
        print(f"  {'offset':>7} {'n':>4} {'mean':>9} {'se':>7} {'t':>6} "
              f"{'up':>10} {'mean abs':>9} {'perm p':>8}")
        for offset in OFFSETS:
            picked = [moves[trading[position[d] + offset]] for d in events
                      if d in position and 0 <= position[d] + offset < len(trading)]
            stat = describe(picked)
            stat["p"] = permutation_p(stat["mean"], pool, stat["n"], rng)
            stat.update(series=series_id, offset=offset)
            rows.append(stat)
            share = 100 * stat["up"] / stat["n"]
            print(f"  {offset:>+7} {stat['n']:>4} {stat['mean']:>+9.2f} "
                  f"{stat['se']:>7.2f} {stat['t']:>+6.2f} "
                  f"{stat['up']:>4}/{stat['n']:<3} ({share:4.1f}%) "
                  f"{stat['abs']:>9.2f} {stat['p']:>8.4f}")
        rest = describe(quiet)
        print(f"  {'other':>7} {rest['n']:>4} {rest['mean']:>+9.2f} "
              f"{rest['se']:>7.2f} {'':>6} {rest['up']:>4}/{rest['n']:<3} "
              f"({100 * rest['up'] / rest['n']:4.1f}%) {rest['abs']:>9.2f}")
        print(f"  absolute-move ratio on the announcement day: "
              f"{[r for r in rows if r['series'] == series_id and r['offset'] == 0][0]['abs'] / rest['abs']:.2f}x")
        print()

    tests = len(rows)
    print(f"MULTIPLE TESTING: {tests} tests were run "
          f"({len(SERIES)} series x {len(OFFSETS)} offsets).")
    print(f"  At the 5% level, {0.05 * tests:.1f} of them are expected to look "
          f"significant with nothing going on.")
    print(f"  A Bonferroni threshold for a 5% family-wide rate is "
          f"{0.05 / tests:.5f}.")
    survivors = [r for r in rows if r["p"] < 0.05 / tests]
    nominal = [r for r in rows if r["p"] < 0.05]
    print(f"  Nominally significant: "
          + ", ".join(f"{r['series']} t{r['offset']:+d} p={r['p']:.4f}"
                      for r in sorted(nominal, key=lambda r: r["p"])))
    print(f"  Surviving Bonferroni: "
          + (", ".join(f"{r['series']} t{r['offset']:+d} p={r['p']:.4f}"
                       for r in survivors) or "none"))
    print()

    series_id, offset = HEADLINE
    moves = changes(fred(series_id), SERIES[series_id])
    trading = [d for d in sorted(moves) if d >= f"{args.since}-01-01"]
    position = {d: i for i, d in enumerate(trading)}
    picked = {d: moves[trading[position[d] + offset]] for d in events
              if d in position and 0 <= position[d] + offset < len(trading)}
    print(f"ROBUSTNESS on {series_id} at t{offset:+d}")
    print(f"  {'subset':<30} {'n':>4} {'mean':>9} {'t':>6} {'up':>9}")

    def line(label, keys):
        values = [picked[k] for k in keys]
        stat = describe(values)
        print(f"  {label:<30} {stat['n']:>4} {stat['mean']:>+9.2f} "
              f"{stat['t']:>+6.2f} {stat['up']:>4}/{stat['n']:<4}")
        return stat

    everything = line("all", list(picked))
    thirds = sorted(picked)
    cut = len(thirds) // 3
    for label, keys in (("first third by date", thirds[:cut]),
                        ("second third", thirds[cut:2 * cut]),
                        ("final third", thirds[2 * cut:])):
        line(f"{label} ({keys[0][:4]}-{keys[-1][:4]})", keys)
    line("excluding 2008 and 2009", [k for k in picked if k[:4] not in ("2008", "2009")])
    line("excluding 2020", [k for k in picked if k[:4] != "2020"])
    line("excluding 2008, 2009, 2020",
         [k for k in picked if k[:4] not in ("2008", "2009", "2020")])
    biggest = sorted(picked, key=lambda k: -abs(picked[k]))
    for drop in (5, 10):
        line(f"dropping the {drop} largest moves",
             [k for k in picked if k not in biggest[:drop]])
    print(f"  the {5} largest: "
          + ", ".join(f"{k} {picked[k]:+.0f}" for k in biggest[:5]))

    years = (int(max(picked)[:4]) - int(args.since)) + int(max(picked)[5:7]) / 12
    rate = len(picked) / years
    print(f"\n  {len(picked)} events over {years:.2f} years is {rate:.2f} a year, "
          f"so a mean of {everything['mean']:.2f} bp is "
          f"{everything['mean'] * rate / 100:.2f}% a year before costs.")

    if args.csv:
        with open(args.csv, "w", encoding="utf-8", newline="") as handle:
            writer = csv.writer(handle)
            writer.writerow(["series", "offset", "n", "mean_bp", "stderr_bp",
                             "t_stat", "days_up", "share_up_pct",
                             "mean_abs_bp", "median_bp", "permutation_p"])
            for r in rows:
                writer.writerow([r["series"], r["offset"], r["n"],
                                 f"{r['mean']:.4f}", f"{r['se']:.4f}",
                                 f"{r['t']:.4f}", r["up"],
                                 f"{100 * r['up'] / r['n']:.2f}",
                                 f"{r['abs']:.4f}", f"{r['median']:.4f}",
                                 f"{r['p']:.5f}"])
        print(f"\nwrote {args.csv}")


if __name__ == "__main__":
    main()