"""correlation-instability-rates-equities-dollar.py — how much the correlations between the 10-year Treasury yield,
U.S. equities, and the dollar move around, and how much of that movement a constant-correlation world would produce
on its own.

Input : FRED's public CSV endpoint (https://fred.stlouisfed.org/graph/fredgraph.csv?id=<ID>), no API key.
        Series pulled with --download:
          DGS10      Market Yield on U.S. Treasury Securities at 10-Year Constant Maturity, percent, daily (H.15)
          NASDAQCOM  NASDAQ Composite, index Feb 5 1971=100, daily close (Nasdaq, Inc.)
          DTWEXBGS   Nominal Broad U.S. Dollar Index, index Jan 2006=100, daily (H.10)
          GS10       the monthly 10-year constant maturity series, percent (H.15) — used only as a method check
          TWEXBGSMTH the monthly broad dollar index (G.5) — used only as a method check
        THE NASDAQ CLOSES ARE NEVER WRITTEN TO DISK. FRED tags NASDAQCOM "Copyrighted: Pre-Approval Required" and
        its series notes read "Copyright © 2016, NASDAQ OMX Group, Inc." (source: Nasdaq, Inc.;
        https://fred.stlouisfed.org/series/NASDAQCOM, read 2026-10-10). --download holds the closes in memory and
        writes only r_NASDAQCOM_pct, the daily return computed from them, which is this site's own output. FRED's
        S&P 500 series is not used either: its notes read "Reproduction of S&P 500 in any form is prohibited except
        with the prior written permission of S&P Dow Jones Indices LLC", and FRED carries only ten years of it.
        Without --download the script reads the two input files that --download writes:
          datasets/correlation-instability-rates-equities-dollar-raw.csv  DGS10, DTWEXBGS, GS10 and TWEXBGSMTH as
                                                                          FRED served them, "." as an empty cell
          datasets/correlation-instability-rates-equities-dollar.csv      the aligned daily changes, the Nasdaq
                                                                          return at full round-trip precision so
                                                                          that every figure reproduces exactly
        One printed count needs the Nasdaq calendar, which no shipped file carries: how many Nasdaq closes fell
        inside the sample window. An offline run quotes it from the 2026-09-18 pull; --download recounts it.
        NOTE ON THE REQUEST: fredgraph.csv stalls and then times out if the request carries no Accept header, which
        is what urllib sends by default. Sending "Accept: */*" returns the file immediately.
Sample: every date from START to END on which all three daily series have a value, and the change between each
        consecutive pair of those dates. END is the last date the dollar index had reached when the data were
        pulled for the article (2026-09-18) and the method check stops at CHECK_END, so a later pull reproduces the
        published figures; move both to extend the study.
Method: Daily changes are computed on the intersection of dates where all three series have a value, between
        consecutive dates in that intersection: the yield in percentage points, the two indexes in percent.
          rolling correlation = Pearson r over a trailing window of WINDOW observations, computed from running sums.
          Fisher-z interval  = tanh(atanh(r) +/- 1.96/sqrt(n-3)), the textbook interval for one window.
          noise test         = simulate SIMS paths of bivariate normal pairs with the SAMPLE correlation held
                               constant for the whole path, same length as the data, and compare the observed
                               spread of rolling correlations with the simulated distribution. A path is drawn as
                               x = z1, y = rho*z1 + sqrt(1-rho^2)*z2 (seeded; stdlib random only).
          horizons           = the same correlations recomputed on weekly (Friday-to-Friday) and monthly
                               (last-observation-to-last-observation) changes. Every statistic here is
                               scale-invariant, so the Nasdaq path for these is rebuilt from its daily returns,
                               starting at 1.0 on the first date of the sample.
        Method check against a published figure: H.15 states that monthly rates are averages of business days, so
        averaging the DGS10 observations in each month must reproduce the published monthly series GS10. The same
        test is run on DTWEXBGS against TWEXBGSMTH, where no such statement was found on the source pages.
        Nothing here is a forecast or investment advice; it describes what the series did over the sample.
Output: prints every number used in the article; writes
          datasets/correlation-instability-rates-equities-dollar-rolling.csv  rolling correlations, one row per date
          datasets/correlation-instability-rates-equities-dollar-annual.csv   correlations by calendar year
          datasets/correlation-instability-rates-equities-dollar-monthly.csv  the monthly-average method check
        and, with --download, the two input files above.
        With --figures, writes two PNG charts to content/images/correlation-instability-rates-equities-dollar/.
Run   : python code/correlation-instability-rates-equities-dollar.py [--download] [--figures]
Needs : Python 3.13. Standard library only for every number. matplotlib 3.10.9 only for --figures.
"""
from __future__ import annotations

import csv
import datetime
import math
import os
import random
import statistics
import sys
import time
import urllib.request

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
SLUG = "correlation-instability-rates-equities-dollar"
DS = os.path.join(ROOT, "datasets")
IMG_DIR = os.path.join(ROOT, "content", "images", SLUG)

DAILY = ["DGS10", "NASDAQCOM", "DTWEXBGS"]
EQUITY = "NASDAQCOM"        # held in memory by --download, never written (see the header)
MONTHLY_CHECK = [("DGS10", "GS10"), ("DTWEXBGS", "TWEXBGSMTH")]
SHIPPED = ["DGS10", "DTWEXBGS", "GS10", "TWEXBGSMTH"]   # the series the raw file carries
RAW = os.path.join(DS, SLUG + "-raw.csv")          # the four Federal Reserve series as pulled, wide
OUT_DAILY = os.path.join(DS, SLUG + ".csv")        # the aligned daily changes: written by --download, read by every run
OUT_ROLL = os.path.join(DS, SLUG + "-rolling.csv")
OUT_YEAR = os.path.join(DS, SLUG + "-annual.csv")
OUT_MONTH = os.path.join(DS, SLUG + "-monthly.csv")

START = "2006-01-02"        # DTWEXBGS begins here
END = "2026-09-11"          # the last date DTWEXBGS had reached on 2026-09-18, the day the article's data was pulled
CHECK_END = "2026-08"       # the last month GS10 and TWEXBGSMTH had reached on 2026-09-18
NASDAQ_CLOSES_2026_09_18 = 5206   # NASDAQCOM closes dated START to END in the 2026-09-18 pull, quoted offline
WINDOW = 126                # trailing observations, about six months of trading days
SIMS = 1500
SEED = 20260918
# the three pairs, in a fixed order used by every table and figure
PAIRS = [("d_DGS10", "r_NASDAQCOM", "10-year yield change vs equity return"),
         ("r_NASDAQCOM", "r_DTWEXBGS", "equity return vs dollar return"),
         ("d_DGS10", "r_DTWEXBGS", "10-year yield change vs dollar return")]
UA = {"User-Agent": "prismdatalab-research/1.0", "Accept": "*/*"}


# --------------------------------------------------------------- download ---
def fetch(series: str, tries: int = 4) -> list[tuple[str, str]]:
    url = f"https://fred.stlouisfed.org/graph/fredgraph.csv?id={series}"
    for i in range(tries):
        try:
            with urllib.request.urlopen(urllib.request.Request(url, headers=UA), timeout=120) as r:
                text = r.read().decode("utf-8")
            break
        except Exception as exc:  # noqa: BLE001
            if i == tries - 1:
                raise
            print(f"  retry in {5 * (i + 1)}s after {exc}", file=sys.stderr)
            time.sleep(5 * (i + 1))
    rows = [r for r in csv.reader(text.splitlines()) if r]
    head = rows[0]
    if len(head) != 2 or head[1].strip() != series:
        sys.exit(f"unexpected columns {head} for {series}")
    return [(r[0], "" if r[1] in (".", "") else r[1]) for r in rows[1:]]


def download() -> int:
    """Pull the five series and write the two input files. Returns the count of Nasdaq closes in the window."""
    pulled: dict[str, list[tuple[str, str]]] = {}
    for series in DAILY + [m for _, m in MONTHLY_CHECK]:
        obs = fetch(series)
        vals = [o for o in obs if o[1]]
        print(f"{series}: {len(obs)} rows, {len(vals)} values, {vals[0][0]} to {vals[-1][0]}")
        pulled[series] = obs
        time.sleep(1.5)     # FRED's robots.txt asks for a one-second crawl delay
    return write_inputs(pulled)


def latest_common_before(levels: dict[str, dict], day: str) -> str:
    """the last date before `day` on which both Federal Reserve daily series have a value"""
    return max(d for d in set(levels["DGS10"]) & set(levels["DTWEXBGS"]) if d < day)


def write_inputs(pulled: dict[str, list[tuple[str, str]]]) -> int:
    """Write the raw file (the four Federal Reserve series only) and the aligned daily changes. The Nasdaq closes
    are used here, in memory, for the return between consecutive sample dates, and are not written anywhere."""
    table: dict[str, dict[str, str]] = {}
    for series in SHIPPED:
        for d, v in pulled[series]:
            table.setdefault(d, {})[series] = v
    with open(RAW, "w", encoding="utf-8", newline="") as f:
        w = csv.writer(f)
        w.writerow(["date"] + SHIPPED)
        for d in sorted(table):
            w.writerow([d] + [table[d].get(c, "") for c in SHIPPED])
    print(f"  wrote {RAW} ({len(table)} dates x {len(SHIPPED)} series, "
          f"retrieved {datetime.date.today().isoformat()}; {EQUITY} is not written)")

    lv = {s: {d: float(v) for d, v in pulled[s] if v} for s in DAILY}
    dates = sorted(d for d in set(lv["DGS10"]) & set(lv[EQUITY]) & set(lv["DTWEXBGS"]) if START <= d <= END)
    # An offline run takes the first sample date to be the last date before the first change on which both Federal
    # Reserve series have a value, because the shipped files carry no Nasdaq calendar. Stop if that would be wrong.
    if latest_common_before(lv, dates[1]) != dates[0]:
        sys.exit(f"the Nasdaq had no close on {latest_common_before(lv, dates[1])}: move START so that an offline "
                 f"run can find the first sample date ({dates[0]})")
    g, n, x = lv["DGS10"], lv[EQUITY], lv["DTWEXBGS"]
    with open(OUT_DAILY, "w", encoding="utf-8", newline="") as f:
        w = csv.writer(f)
        w.writerow(["date", "DGS10_pct", "DTWEXBGS_index", "d_DGS10_pp", "r_NASDAQCOM_pct", "r_DTWEXBGS_pct"])
        for prev, cur in zip(dates, dates[1:]):
            w.writerow([cur, f"{g[cur]:.2f}", f"{x[cur]:.4f}", f"{g[cur] - g[prev]:.4f}",
                        repr(100 * (n[cur] / n[prev] - 1)), f"{100 * (x[cur] / x[prev] - 1):.6f}"])
    print(f"  wrote {OUT_DAILY} ({len(dates) - 1} daily changes; the {EQUITY} column is the return, not the index)")
    return sum(1 for d in n if START <= d <= dates[-1])


def load_raw() -> dict[str, dict[str, float]]:
    """{series: {date: value}}"""
    with open(RAW, encoding="utf-8", newline="") as f:
        rows = list(csv.reader(f))
    names = rows[0][1:]
    out: dict[str, dict[str, float]] = {n: {} for n in names}
    for r in rows[1:]:
        for n, v in zip(names, r[1:]):
            if v != "":
                out[n][r[0]] = float(v)
    return out


def load_equity_returns() -> tuple[list[str], list[float]]:
    """the sample's change dates and the Nasdaq daily return on each, from the aligned daily-changes file"""
    with open(OUT_DAILY, encoding="utf-8", newline="") as f:
        rows = list(csv.DictReader(f))
    return [r["date"] for r in rows], [float(r["r_NASDAQCOM_pct"]) for r in rows]


# ------------------------------------------------------------ statistics ---
def pearson(xs: list[float], ys: list[float]) -> float:
    n = len(xs)
    sx, sy = sum(xs), sum(ys)
    sxx = sum(x * x for x in xs)
    syy = sum(y * y for y in ys)
    sxy = sum(x * y for x, y in zip(xs, ys))
    num = n * sxy - sx * sy
    den = math.sqrt((n * sxx - sx * sx) * (n * syy - sy * sy))
    return num / den if den else float("nan")


def rolling_corr(xs: list[float], ys: list[float], window: int) -> list[float | None]:
    """trailing Pearson r ending at each index, from running sums (None until the window fills)"""
    out: list[float | None] = [None] * len(xs)
    sx = sy = sxx = syy = sxy = 0.0
    for i, (x, y) in enumerate(zip(xs, ys)):
        sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y
        if i >= window:
            a, b = xs[i - window], ys[i - window]
            sx -= a; sy -= b; sxx -= a * a; syy -= b * b; sxy -= a * b
        if i >= window - 1:
            n = window
            den = math.sqrt((n * sxx - sx * sx) * (n * syy - sy * sy))
            out[i] = (n * sxy - sx * sy) / den if den else None
    return out


def fisher_ci(r: float, n: int, z: float = 1.959964) -> tuple[float, float]:
    se = 1.0 / math.sqrt(n - 3)
    lo, hi = math.atanh(r) - z * se, math.atanh(r) + z * se
    return math.tanh(lo), math.tanh(hi)


def sign_crossings(series: list[float | None]) -> int:
    vals = [v for v in series if v is not None]
    return sum(1 for a, b in zip(vals, vals[1:]) if a * b < 0)


def simulate_range(rho: float, n: int, window: int, sims: int, seed: int) -> dict:
    """rolling-correlation spread in worlds where the true correlation never changes"""
    rng = random.Random(seed)
    k = math.sqrt(max(0.0, 1.0 - rho * rho))
    ranges, mins, maxs, crossings = [], [], [], []
    for _ in range(sims):
        xs, ys = [], []
        for _ in range(n):
            z1, z2 = rng.gauss(0.0, 1.0), rng.gauss(0.0, 1.0)
            xs.append(z1)
            ys.append(rho * z1 + k * z2)
        r = [v for v in rolling_corr(xs, ys, window) if v is not None]
        ranges.append(max(r) - min(r))
        mins.append(min(r))
        maxs.append(max(r))
        crossings.append(sum(1 for a, b in zip(r, r[1:]) if a * b < 0))
    ranges.sort(); mins.sort(); maxs.sort()
    q = lambda s, p: s[min(len(s) - 1, int(p * len(s)))]  # noqa: E731
    return {"range_med": statistics.median(ranges), "range_p95": q(ranges, 0.95),
            "range_max": ranges[-1], "min_p05": q(mins, 0.05), "max_p95": q(maxs, 0.95),
            "cross_med": statistics.median(crossings), "ranges": ranges}


# ------------------------------------------------------------------ main ---
def build_changes(raw, change_dates, equity_returns) -> tuple[list[str], dict[str, list[float]], list[str]]:
    """the yield and dollar changes between consecutive sample dates, beside the shipped Nasdaq returns"""
    sample = [latest_common_before(raw, change_dates[0])] + change_dates
    cols: dict[str, list[float]] = {"d_DGS10": [], "r_NASDAQCOM": list(equity_returns), "r_DTWEXBGS": []}
    for prev, cur in zip(sample, sample[1:]):
        cols["d_DGS10"].append(raw["DGS10"][cur] - raw["DGS10"][prev])
        cols["r_DTWEXBGS"].append(100 * (raw["DTWEXBGS"][cur] / raw["DTWEXBGS"][prev] - 1))
    return change_dates, cols, sample


def period_corr(dates, cols, keys, pick) -> tuple[float, int]:
    """correlation of changes between the picked period-end dates (weekly/monthly)"""
    raw_levels = pick
    xs = [raw_levels[k] for k in keys]
    return pearson(*xs), len(xs[0])


def main() -> None:
    live_count = download() if "--download" in sys.argv else None
    raw = load_raw()
    dates, cols, sample = build_changes(raw, *load_equity_returns())
    n_common = len(sample)
    n = len(dates)
    print(f"\naligned sample: {n} daily changes, {dates[0]} to {dates[-1]} "
          f"({n_common} dates carried all three series)")
    for s in DAILY:
        if s == EQUITY:
            k = live_count if live_count is not None else NASDAQ_CLOSES_2026_09_18
            note = "" if live_count is not None else " (the 2026-09-18 pull's count; --download recounts it)"
        else:
            k, note = len([d for d in raw[s] if START <= d <= dates[-1]]), ""
        print(f"  {s}: {k} observations in the window; {k - n_common} not in the common set{note}")

    # ---- full sample, and the rolling series
    print(f"\nfull-sample correlation of daily changes (n={n}), 95% Fisher interval:")
    rolls: dict[str, list[float | None]] = {}
    full: dict[str, float] = {}
    for a, b, label in PAIRS:
        r = pearson(cols[a], cols[b])
        lo, hi = fisher_ci(r, n)
        full[label] = r
        print(f"  {label:<38} r = {r:+.3f}  [{lo:+.3f}, {hi:+.3f}]")
        rolls[label] = rolling_corr(cols[a], cols[b], WINDOW)

    print(f"\nrolling {WINDOW}-observation correlation:")
    for _, _, label in PAIRS:
        vals = [(d, v) for d, v in zip(dates, rolls[label]) if v is not None]
        vs = [v for _, v in vals]
        lo_d, lo_v = min(vals, key=lambda t: t[1])
        hi_d, hi_v = max(vals, key=lambda t: t[1])
        neg = sum(1 for v in vs if v < 0)
        print(f"  {label}")
        print(f"    {len(vs)} windows | mean {statistics.fmean(vs):+.3f} | sd {statistics.pstdev(vs):.3f} | "
              f"min {lo_v:+.3f} ({lo_d}) | max {hi_v:+.3f} ({hi_d}) | range {hi_v - lo_v:.3f}")
        print(f"    negative in {neg} of {len(vs)} windows ({100 * neg / len(vs):.1f}%) | "
              f"sign crossings {sign_crossings(rolls[label])}")

    print("\n  the latest window, with its own 95% interval:")
    for _, _, label in PAIRS:
        vals = [(d, v) for d, v in zip(dates, rolls[label]) if v is not None]
        d_last, v_last = vals[-1]
        lo, hi = fisher_ci(v_last, WINDOW)
        print(f"    {label:<38} {v_last:+.3f} on {d_last}  [{lo:+.3f}, {hi:+.3f}]")
    label0 = PAIRS[0][2]
    runs, cur = [], None
    for d, v in [(d, v) for d, v in zip(dates, rolls[label0]) if v is not None]:
        if v < 0:
            cur = [d, d, 1] if cur is None else [cur[0], d, cur[2] + 1]
        elif cur is not None:
            runs.append(cur)
            cur = None
    if cur is not None:
        runs.append(cur)
    runs.sort(key=lambda r: -r[2])
    print(f"  longest unbroken stretches with a negative {label0} ({len(runs)} stretches in all):")
    for a, b, k in runs[:3]:
        print(f"    {k} windows, {a} to {b}")

    with open(OUT_ROLL, "w", encoding="utf-8", newline="") as f:
        w = csv.writer(f)
        w.writerow(["date", "corr_yield_equity", "corr_equity_dollar", "corr_yield_dollar", "window_obs"])
        for i, d in enumerate(dates):
            vals = [rolls[label][i] for _, _, label in PAIRS]
            if any(v is None for v in vals):
                continue
            w.writerow([d] + [f"{v:.6f}" for v in vals] + [WINDOW])

    # ---- one window's sampling error
    lo, hi = fisher_ci(0.0, WINDOW)
    print(f"\nsampling error in a single {WINDOW}-observation window: a sample r of 0.00 has a 95% interval of "
          f"[{lo:+.3f}, {hi:+.3f}]")
    for r0 in (0.20, 0.40):
        lo, hi = fisher_ci(r0, WINDOW)
        print(f"  a sample r of {r0:+.2f} has a 95% interval of [{lo:+.3f}, {hi:+.3f}]")

    # ---- the constant-correlation simulation
    print(f"\nconstant-correlation simulation ({SIMS} paths of {n} pairs, seed {SEED}):")
    for _, _, label in PAIRS:
        vs = [v for v in rolls[label] if v is not None]
        obs_range = max(vs) - min(vs)
        sim = simulate_range(full[label], n, WINDOW, SIMS, SEED)
        beat = sum(1 for x in sim["ranges"] if x >= obs_range)
        print(f"  {label}")
        print(f"    true r fixed at {full[label]:+.3f}: simulated rolling range median {sim['range_med']:.3f}, "
              f"95th pct {sim['range_p95']:.3f}, max {sim['range_max']:.3f}")
        print(f"    simulated 5th pct of the minimum {sim['min_p05']:+.3f}, 95th pct of the maximum "
              f"{sim['max_p95']:+.3f}, median sign crossings {sim['cross_med']:.0f}")
        print(f"    observed range {obs_range:.3f}: {beat} of {SIMS} constant-correlation paths reached it "
              f"({100 * beat / SIMS:.2f}%)")

    # ---- by calendar year
    print("\ncorrelation of daily changes by calendar year:")
    print("  year |    n | yield-equity | equity-dollar | yield-dollar")
    years = sorted({d[:4] for d in dates})
    year_rows = []
    for y in years:
        idx = [i for i, d in enumerate(dates) if d.startswith(y)]
        rs = [pearson([cols[a][i] for i in idx], [cols[b][i] for i in idx]) for a, b, _ in PAIRS]
        year_rows.append((y, len(idx), rs))
        print(f"  {y} | {len(idx):>4} | {rs[0]:+12.3f} | {rs[1]:+13.3f} | {rs[2]:+12.3f}")
    with open(OUT_YEAR, "w", encoding="utf-8", newline="") as f:
        w = csv.writer(f)
        w.writerow(["year", "n_days", "corr_yield_equity", "corr_equity_dollar", "corr_yield_dollar"])
        for y, cnt, rs in year_rows:
            w.writerow([y, cnt] + [f"{r:.6f}" for r in rs])

    # ---- horizon dependence
    print("\nsame period, different horizons (correlation of changes):")
    path = [1.0]                                    # the Nasdaq rebased to 1.0 on the first sample date
    for r in cols["r_NASDAQCOM"]:
        path.append(path[-1] * (1 + r / 100))
    equity_level = dict(zip(sample, path))
    weekly_dates, monthly_dates = [], []
    for i, d in enumerate(dates):
        dt = datetime.date.fromisoformat(d)
        nxt = datetime.date.fromisoformat(dates[i + 1]) if i + 1 < len(dates) else None
        if nxt is None or nxt.isocalendar()[:2] != dt.isocalendar()[:2]:
            weekly_dates.append(i)
        if nxt is None or nxt.month != dt.month:
            monthly_dates.append(i)
    for label, idxs in (("daily", list(range(len(dates)))), ("weekly", weekly_dates), ("monthly", monthly_dates)):
        if label == "daily":
            series = {k: cols[k] for k in cols}
        else:
            series = {"d_DGS10": [], "r_NASDAQCOM": [], "r_DTWEXBGS": []}
            for a, b in zip(idxs, idxs[1:]):
                d0, d1 = dates[a], dates[b]
                series["d_DGS10"].append(raw["DGS10"][d1] - raw["DGS10"][d0])
                series["r_NASDAQCOM"].append(100 * (equity_level[d1] / equity_level[d0] - 1))
                series["r_DTWEXBGS"].append(100 * (raw["DTWEXBGS"][d1] / raw["DTWEXBGS"][d0] - 1))
        rs = [pearson(series[a], series[b]) for a, b, _ in PAIRS]
        cnt = len(series["d_DGS10"])
        print(f"  {label:<8} n={cnt:>5} | yield-equity {rs[0]:+.3f} | equity-dollar {rs[1]:+.3f} | "
              f"yield-dollar {rs[2]:+.3f}")

    # ---- method check against published monthly series
    # Both the daily and the monthly series are published to a fixed number of decimals, so the whole comparison is
    # done in integer units of that last decimal: no float rounding can decide the answer. The mean of cnt values
    # totalling `total` units, rounded half up, is (2*total + cnt) // (2*cnt).
    print("\nmethod check: averaging daily business-day values against the published monthly series")
    month_rows = []
    for daily_id, monthly_id in MONTHLY_CHECK:
        dec = 2 if daily_id == "DGS10" else 4
        q = 10 ** dec
        buckets: dict[str, list[int]] = {}
        for d, v in raw[daily_id].items():
            buckets.setdefault(d[:7], []).append(round(v * q))
        exact, halfway, rows = 0, [], []
        for m in sorted(buckets):
            key = m + "-01"
            if m > CHECK_END or key not in raw[monthly_id]:
                continue
            units = buckets[m]
            total, cnt = sum(units), len(units)
            mine_units = (2 * total + cnt) // (2 * cnt)
            pub_units = round(raw[monthly_id][key] * q)
            exact += 1 if mine_units == pub_units else 0
            if (2 * total) % (2 * cnt) == cnt:          # the average lands exactly on half a unit
                halfway.append(m)
            rows.append((monthly_id, m, cnt, total / (cnt * q), pub_units / q,
                         (mine_units - pub_units) / q, mine_units - pub_units))
        print(f"  {daily_id} averaged by month vs {monthly_id}: {len(rows)} months compared "
              f"({rows[0][1]} to {rows[-1][1]})")
        print(f"    rounding half up, {exact} of {len(rows)} reproduce the published value exactly at "
              f"{dec} decimals ({100 * exact / len(rows):.1f}%)")
        off = [(r[1], r[6]) for r in rows if r[6] != 0]
        if off:
            print(f"    the {len(off)} that differ are each one unit in the last decimal: {off}")
        print(f"    {len(halfway)} monthly averages land exactly on half a unit, where Python's round() rounds to "
              f"even and the published series rounds away from zero: {halfway}")
        month_rows += rows
    with open(OUT_MONTH, "w", encoding="utf-8", newline="") as f:
        w = csv.writer(f)
        w.writerow(["published_series", "month", "n_business_days", "computed_average", "published_value",
                    "difference", "difference_last_decimal_units"])
        for sid, m, cnt, mine, pub, diff, units in month_rows:
            dec = 2 if sid == "GS10" else 4
            w.writerow([sid, m, cnt, f"{mine:.6f}", f"{pub:.{dec}f}", f"{diff:.6f}", units])

    if "--figures" in sys.argv:
        figures(dates, rolls, full, n)


# ------------------------------------------------------------------ figures ---
def figures(dates, rolls, full, n) -> None:
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    os.makedirs(IMG_DIR, exist_ok=True)
    paper, ink, soft, grid = "#ffffff", "#1b2430", "#5e646b", "#e3e2de"
    c1, c2, c3 = "#2a78d6", "#eb6834", "#0f766e"
    plt.rcParams.update({"font.family": "DejaVu Sans", "font.size": 10, "axes.edgecolor": grid,
                         "axes.labelcolor": ink, "xtick.color": soft, "ytick.color": soft, "text.color": ink,
                         "axes.spines.top": False, "axes.spines.right": False})
    x = [datetime.date.fromisoformat(d) for d in dates]

    fig, ax = plt.subplots(figsize=(11, 4.8), dpi=120, facecolor=paper)
    ax.set_facecolor(paper)
    ax.grid(axis="y", color=grid, linewidth=0.8)
    ax.axhline(0, color=ink, linewidth=1.0)
    for (_, _, label), col in zip(PAIRS, (c1, c2, c3)):
        ys = rolls[label]
        xs = [d for d, v in zip(x, ys) if v is not None]
        vs = [v for v in ys if v is not None]
        ax.plot(xs, vs, color=col, linewidth=1.6, label=label)
    ax.set_ylabel(f"trailing {WINDOW}-day correlation")
    ax.set_ylim(-1, 1)
    ax.set_title(f"Rolling {WINDOW}-observation correlation of daily changes, {dates[0]} to {dates[-1]}",
                 loc="left", fontsize=11)
    ax.legend(frameon=False, loc="lower left", fontsize=9)
    fig.tight_layout()
    fig.savefig(os.path.join(IMG_DIR, "rolling-correlations.png"), facecolor=paper)
    plt.close(fig)

    # the stock-bond pair against what a constant correlation would produce
    label = PAIRS[0][2]
    sim = simulate_range(full[label], n, WINDOW, 600, SEED + 1)
    fig, ax = plt.subplots(figsize=(11, 4.4), dpi=120, facecolor=paper)
    ax.set_facecolor(paper)
    ax.grid(axis="y", color=grid, linewidth=0.8)
    ax.axhline(0, color=ink, linewidth=1.0)
    ys = rolls[label]
    xs = [d for d, v in zip(x, ys) if v is not None]
    vs = [v for v in ys if v is not None]
    ax.axhspan(sim["min_p05"], sim["max_p95"], color=c1, alpha=0.13,
               label="where a constant correlation keeps 90% of its extremes")
    ax.axhline(full[label], color=soft, linewidth=1.0, linestyle="--",
               label=f"full-sample correlation {full[label]:+.3f}")
    ax.plot(xs, vs, color=c2, linewidth=1.7, label=f"trailing {WINDOW}-day correlation")
    ax.set_ylabel("correlation")
    ax.set_ylim(-1, 1)
    ax.set_title("10-year yield change against equity return: the observed swing versus sampling noise",
                 loc="left", fontsize=11)
    ax.legend(frameon=False, loc="lower left", fontsize=9)
    fig.tight_layout()
    fig.savefig(os.path.join(IMG_DIR, "stock-bond-vs-noise.png"), facecolor=paper)
    plt.close(fig)
    print(f"\nwrote 2 figures to {IMG_DIR}")


if __name__ == "__main__":
    main()