"""real-versus-nominal-deflating-correctly.py - deflating a nominal series, checked against the BLS's own real
earnings series, and six mistakes measured on public data.

Input : FRED's public CSV endpoint, no API key (https://fred.stlouisfed.org/graph/fredgraph.csv?id=<ID>):
          monthly   CES0500000003 AHETPI CPIAUCSL CPIAUCNS CWSR0000SA0 PCEPI CPILFESL PCEPILFE FEDFUNDS
          quarterly GDP GDPC1 GDPDEF PCEC PCECC96 GPDI GPDIC1 GCE GCEC1 NETEXP NETEXC
        and the BLS Public Data API v1, no key (https://api.bls.gov/publicAPI/v1/timeseries/data/), for the
        BLS's own constant-dollar series: CES0500000013 (average hourly earnings of all employees in 1982-1984
        dollars, deflated by the CPI-U) and CES0500000032 (production and nonsupervisory employees, deflated by
        the CPI-W). With --download the script fetches all of them and rewrites the two CSVs below; without it,
        it reads the CSVs, so every number can be reproduced offline.
Method: real_t = nominal_t / (P_t / P_ref), where P is a price index and P_ref its value in the month whose
        prices you want to express everything in (P_ref = 100 gives the index's own base period, which for the
        CPI is 1982-1984). A real change between two months is (1 + nominal change) / (1 + inflation) - 1.
        Every statistic is computed from unrounded values; tables print them rounded. The BLS check rounds
        to the cent, half up, because that is how the BLS publishes its real series.
Output: prints every table in the article and writes
          datasets/real-versus-nominal-deflating-correctly.csv      monthly, 1954-07 to the latest month
          datasets/real-versus-nominal-deflating-correctly-gdp.csv  quarterly, 1947Q1 to the latest quarter
        With --figures, writes one PNG chart to content/images/real-versus-nominal-deflating-correctly/.
Run   : python code/real-versus-nominal-deflating-correctly.py [--download] [--figures]
Needs : Python 3.10+. Standard library only for every number. matplotlib (3.10) only for --figures.
        The API needs no key; the script still names its requester and a contact in the User-Agent (CONTACT),
        so change that line to your own before running --download.
Not investment advice; a description of public statistics.
"""
from __future__ import annotations

import csv
import datetime
import json
import os
import statistics
import sys
import time
import urllib.request
from decimal import ROUND_HALF_UP, Decimal

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
SLUG = "real-versus-nominal-deflating-correctly"
MONTHLY_CSV = os.path.join(ROOT, "datasets", SLUG + ".csv")
GDP_CSV = os.path.join(ROOT, "datasets", SLUG + "-gdp.csv")
IMG_DIR = os.path.join(ROOT, "content", "images", SLUG)
CONTACT = "Prism Data Lab research script (contact@prismdatalab.com)"

# dataset column -> source series (FRED unless marked BLS)
MONTHLY = {"ahe_all": "CES0500000003", "ahe_pns": "AHETPI", "cpi_u": "CPIAUCSL", "cpi_u_nsa": "CPIAUCNS",
           "cpi_w": "CWSR0000SA0", "pce_price": "PCEPI", "core_cpi": "CPILFESL", "core_pce": "PCEPILFE",
           "fedfunds": "FEDFUNDS", "bls_real_ahe_all": "BLS:CES0500000013", "bls_real_ahe_pns": "BLS:CES0500000032"}
QUARTERLY = {"gdp": "GDP", "gdp_real": "GDPC1", "gdp_deflator": "GDPDEF", "pce": "PCEC", "pce_real": "PCECC96",
             "investment": "GPDI", "investment_real": "GPDIC1", "government": "GCE", "government_real": "GCEC1",
             "net_exports": "NETEXP", "net_exports_real": "NETEXC"}
MONTHLY_FROM = "1954-07-01"  # the first FEDFUNDS month

# The August 2026 Real Earnings release (USDL-26-1497, 11 September 2026), Table A-1, as published that day.
# Average hourly earnings for August 2026 was the first estimate; the 2 October jobs report revised it to 37.76.
RELEASE = {"ahe_2025_08": 36.62, "ahe_2026_08_first": 37.75, "cpi_2025_08": 323.291, "cpi_2026_08": 334.131,
           "real_2025_08": 11.33, "real_2026_08": 11.30}


# ----------------------------------------------------------------- download ---
def fred(series_id: str) -> dict[str, str]:
    url = f"https://fred.stlouisfed.org/graph/fredgraph.csv?id={series_id}"
    # fredgraph.csv stalls when a request carries no Accept header, which is urllib's default.
    req = urllib.request.Request(url, headers={"User-Agent": "prismdatalab-research/1.0", "Accept": "*/*"})
    with urllib.request.urlopen(req, timeout=60) as r:
        rows = list(csv.reader(r.read().decode("utf-8").splitlines()))
    assert rows[0][1] == series_id, rows[0]
    return {d: ("" if v == "." else v) for d, v in rows[1:]}


def bls(series_ids: list[str], first_year: int, last_year: int) -> dict[str, dict[str, str]]:
    out: dict[str, dict[str, str]] = {s: {} for s in series_ids}
    for a in range(first_year, last_year + 1, 10):  # API v1 serves at most ten years a request
        body = json.dumps({"seriesid": series_ids, "startyear": str(a), "endyear": str(min(a + 9, last_year))})
        for attempt in range(3):
            req = urllib.request.Request("https://api.bls.gov/publicAPI/v1/timeseries/data/", data=body.encode(),
                                         headers={"Content-Type": "application/json", "User-Agent": CONTACT})
            try:
                with urllib.request.urlopen(req, timeout=90) as r:
                    payload = json.loads(r.read())
                break
            except OSError:
                if attempt == 2:
                    raise
                time.sleep(10)
        assert payload["status"] == "REQUEST_SUCCEEDED", payload.get("message")
        for s in payload["Results"]["series"]:
            for row in s["data"]:
                if row["period"].startswith("M") and row["period"] != "M13":
                    d = f"{row['year']}-{row['period'][1:]}-01"
                    out[s["seriesID"]][d] = "" if row["value"] in ("-", "") else row["value"]
        time.sleep(2)
    return out


def download() -> None:
    monthly = {}
    for col, sid in MONTHLY.items():
        if not sid.startswith("BLS:"):
            monthly[col] = fred(sid)
            time.sleep(0.5)
    real = bls(["CES0500000013", "CES0500000032"], 2006, datetime.date.today().year)
    monthly["bls_real_ahe_all"], monthly["bls_real_ahe_pns"] = real["CES0500000013"], real["CES0500000032"]
    dates = sorted(d for d in set().union(*[set(v) for v in monthly.values()]) if d >= MONTHLY_FROM)
    with open(MONTHLY_CSV, "w", encoding="utf-8", newline="") as f:
        w = csv.writer(f, lineterminator="\n")
        w.writerow(["date"] + list(MONTHLY))
        for d in dates:
            w.writerow([d] + [monthly[c].get(d, "") for c in MONTHLY])
    quarterly = {}
    for col, sid in QUARTERLY.items():
        quarterly[col] = fred(sid)
        time.sleep(0.5)
    dates = sorted(set().union(*[set(v) for v in quarterly.values()]))
    with open(GDP_CSV, "w", encoding="utf-8", newline="") as f:
        w = csv.writer(f, lineterminator="\n")
        w.writerow(["date"] + list(QUARTERLY))
        for d in dates:
            w.writerow([d] + [quarterly[c].get(d, "") for c in QUARTERLY])
    print(f"downloaded {len(MONTHLY)} monthly and {len(QUARTERLY)} quarterly series "
          f"({datetime.date.today().isoformat()})")


# ------------------------------------------------------------------- helpers ---
def load(path: str) -> dict[str, dict[str, float | None]]:
    with open(path, encoding="utf-8", newline="") as f:
        return {r["date"]: {k: (float(v) if v else None) for k, v in r.items() if k != "date"}
                for r in csv.DictReader(f)}


def pct(a: float, b: float) -> float:
    """percent change from a to b"""
    return 100 * (b / a - 1)


def real_pct(nominal_pct: float, inflation_pct: float) -> float:
    """deflate a change: (1 + g) / (1 + pi) - 1, in percent"""
    return 100 * ((1 + nominal_pct / 100) / (1 + inflation_pct / 100) - 1)


def cents(x: Decimal) -> Decimal:
    return x.quantize(Decimal("0.01"), rounding=ROUND_HALF_UP)


def window_label(d: str) -> str:
    return datetime.date.fromisoformat(d).strftime("%b %Y")


# ------------------------------------------------------------------ analysis ---
def analysis(M: dict, Q: dict) -> dict:
    out: dict = {}
    months = sorted(M)
    both = [d for d in months if M[d]["ahe_all"] is not None and M[d]["cpi_u"] is not None]
    end = both[-1]  # the latest month with both pay and the CPI
    out["end"] = end
    print(f"Latest month with average hourly earnings and the CPI-U: {end}")

    # 1. the check: our deflation against the BLS's published real series, to the cent
    print("\n1. BLS real average hourly earnings, recomputed (pay / CPI x 100, rounded half up to the cent)")
    for real_col, pay_col, cpi_col, name in (("bls_real_ahe_all", "ahe_all", "cpi_u", "all employees / CPI-U"),
                                             ("bls_real_ahe_pns", "ahe_pns", "cpi_w", "production and nonsupervisory / CPI-W")):
        compared = exact = blank = 0
        nsa_exact = 0
        for d in months:
            row = M[d]
            if d < "2006-01-01" or row[pay_col] is None:
                continue
            if row[real_col] is None:
                if row[cpi_col] is None and d <= end:
                    blank += 1
                continue
            ours = cents(Decimal(str(row[pay_col])) / Decimal(str(row[cpi_col])) * 100)
            compared += 1
            exact += ours == Decimal(str(row[real_col])).quantize(Decimal("0.01"))
            if real_col == "bls_real_ahe_all" and row["cpi_u_nsa"] is not None:
                nsa = cents(Decimal(str(row[pay_col])) / Decimal(str(row["cpi_u_nsa"])) * 100)
                nsa_exact += nsa == Decimal(str(row[real_col])).quantize(Decimal("0.01"))
        print(f"   {name:40s} months compared {compared}, identical to the cent {exact}, "
              f"months the BLS leaves blank (no CPI) {blank}")
        out[real_col] = (compared, exact, blank)
        if real_col == "bls_real_ahe_all":
            out["nsa_exact"] = nsa_exact
            print(f"   (with the NOT seasonally adjusted CPI-U instead, {nsa_exact} of {compared} months match)")

    # the release's own over-the-year figure
    R = RELEASE
    rel = real_pct(pct(R["ahe_2025_08"], R["ahe_2026_08_first"]), pct(R["cpi_2025_08"], R["cpi_2026_08"]))
    a, b = "2025-08-01", end
    now = real_pct(pct(M[a]["ahe_all"], M[b]["ahe_all"]), pct(M[a]["cpi_u"], M[b]["cpi_u"]))
    rounded = pct(R["real_2025_08"], R["real_2026_08"])
    out.update(release_real=rel, today_real=now, rounded_real=rounded)
    print(f"   Aug 2025 -> Aug 2026 real change: release inputs {rel:.4f}%, today's data {now:.4f}%, "
          f"from the cent-rounded real levels {rounded:.4f}%")

    # 2. units: 1982-84 dollars, August 2026 dollars, 2017 dollars
    s0 = "2006-03-01"
    r84_end = M[end]["ahe_all"] / M[end]["cpi_u"] * 100
    r84_start = M[s0]["ahe_all"] / M[s0]["cpi_u"] * 100
    in_end_dollars = M[s0]["ahe_all"] * M[end]["cpi_u"] / M[s0]["cpi_u"]
    pce_units = M[end]["ahe_all"] / M[end]["pce_price"] * 100
    out.update(r84_end=r84_end, r84_start=r84_start, in_end_dollars=in_end_dollars, pce_units=pce_units,
               g_84=pct(r84_start, r84_end), g_end=pct(in_end_dollars, M[end]["ahe_all"]),
               g_rounded=pct(round(r84_start, 2), round(r84_end, 2)))
    print(f"\n2. Units. {window_label(end)} pay ${M[end]['ahe_all']:.2f} = ${r84_end:.4f} in 1982-84 dollars "
          f"(CPI-U) = ${pce_units:.4f} in 2017 dollars (PCE price index)")
    print(f"   Mar 2006 pay ${M[s0]['ahe_all']:.2f} = ${r84_start:.4f} in 1982-84 dollars = "
          f"${in_end_dollars:.4f} in {window_label(end)} dollars")
    print(f"   real growth: {out['g_84']:.4f}% in 1982-84 dollars, {out['g_end']:.4f}% in {window_label(end)} dollars, "
          f"{out['g_rounded']:.4f}% from the cent-rounded 1982-84 levels")

    # 3. subtracting instead of dividing
    print("\n3. Subtracting inflation instead of dividing (CPI-U)")
    print(f"   {'window':22s} {'nominal':>9s} {'CPI-U':>9s} {'subtract':>9s} {'divide':>9s} {'1+infl':>7s}")
    out["subtract"] = {}
    for start in ("2025-08-01", "2020-02-01", "2006-03-01"):
        g, p = pct(M[start]["ahe_all"], M[end]["ahe_all"]), pct(M[start]["cpi_u"], M[end]["cpi_u"])
        out["subtract"][start] = (g, p, g - p, real_pct(g, p))
        print(f"   {window_label(start) + ' -> ' + window_label(end):22s} {g:9.3f} {p:9.3f} {g - p:9.3f} "
              f"{real_pct(g, p):9.3f} {1 + p / 100:7.4f}")
    # the same identity for an interest rate: the ex-post real fed funds rate
    fisher = []
    for d in months:
        prev = f"{int(d[:4]) - 1}{d[4:]}"
        if prev in M and M[d]["fedfunds"] is not None and M[d]["cpi_u"] and M[prev]["cpi_u"]:
            i, p = M[d]["fedfunds"], pct(M[prev]["cpi_u"], M[d]["cpi_u"])
            fisher.append((d, i, p, i - p, real_pct(i, p)))
    worst = max(fisher, key=lambda t: abs(t[3] - t[4]))
    last = fisher[-1]
    out.update(fisher_worst=worst, fisher_last=last, fisher_n=len(fisher))
    print(f"   real fed funds rate, {len(fisher)} months: largest gap {worst[0]}: rate {worst[1]:.2f}, CPI-U 12-month "
          f"{worst[2]:.3f}%, subtract {worst[3]:.3f}, divide {worst[4]:.3f}")
    print(f"   latest {last[0]}: rate {last[1]:.2f}, CPI-U 12-month {last[2]:.3f}%, subtract {last[3]:.3f}, "
          f"divide {last[4]:.3f}")

    # 4. the deflator is a choice
    print("\n4. Real change in average hourly earnings to the latest month, by deflator (percent)")
    defl = (("cpi_u", "CPI-U"), ("cpi_w", "CPI-W"), ("pce_price", "PCE price index"))
    print(f"   {'from':10s} {'nominal':>9s}" + "".join(f" {n:>16s}" for _, n in defl))
    out["deflators"] = {}
    for start in ("2025-08-01", "2021-01-01", "2020-02-01", "2006-03-01"):
        g = pct(M[start]["ahe_all"], M[end]["ahe_all"])
        vals = {c: real_pct(g, pct(M[start][c], M[end][c])) for c, _ in defl}
        out["deflators"][start] = (g, vals, {c: pct(M[start][c], M[end][c]) for c, _ in defl})
        print(f"   {window_label(start):10s} {g:9.3f}" + "".join(f" {vals[c]:16.3f}" for c, _ in defl))
    # every start month, not just the four above
    starts = [d for d in months if "2006-03-01" <= d < end and int(end[:4]) * 12 + int(end[5:7])
              - int(d[:4]) * 12 - int(d[5:7]) >= 12]
    by_start = []
    for d in starts:
        g = pct(M[d]["ahe_all"], M[end]["ahe_all"])
        cu = real_pct(g, pct(M[d]["cpi_u"], M[end]["cpi_u"])) if M[d]["cpi_u"] else None
        pc = real_pct(g, pct(M[d]["pce_price"], M[end]["pce_price"]))
        cc = real_pct(g, pct(M[d]["core_cpi"], M[end]["core_cpi"])) if M[d]["core_cpi"] else None
        by_start.append((d, cu, pc, cc))
    valid = [t for t in by_start if t[1] is not None]
    assert all(t[1] != 0 and t[2] != 0 for t in valid)
    disagree = [t[0] for t in valid if (t[1] < 0) != (t[2] < 0)]
    out.update(by_start=by_start, valid_starts=len(valid), disagree=disagree,
               cpi_negative=[t[0] for t in valid if t[1] < 0], pce_negative=[t[0] for t in valid if t[2] < 0])
    years = (int(end[:4]) * 12 + int(end[5:7]) - (2006 * 12 + 3)) / 12
    cu_ann = 100 * ((M[end]["cpi_u"] / M["2006-03-01"]["cpi_u"]) ** (1 / years) - 1)
    pc_ann = 100 * ((M[end]["pce_price"] / M["2006-03-01"]["pce_price"]) ** (1 / years) - 1)
    out.update(cu_ann=cu_ann, pc_ann=pc_ann, years=years)
    print(f"   start months a year or more before {window_label(end)}, from Mar 2006: {len(valid)}; the CPI-U and "
          f"the PCE price index disagree on the sign of the real change for {len(disagree)}: {', '.join(disagree)}")
    print(f"   CPI-U-deflated change negative for {len(out['cpi_negative'])} start months, PCE-deflated for "
          f"{len(out['pce_negative'])}")
    print(f"   average annual inflation Mar 2006 -> {window_label(end)} ({years:.4f} years): CPI-U {cu_ann:.4f}%, "
          f"PCE {pc_ann:.4f}%, gap {cu_ann - pc_ann:.4f} points a year")

    # 5. core is not a deflator
    print("\n5. The latest year, headline against core (percent)")
    a = "2025-08-01"
    g = pct(M[a]["ahe_all"], M[end]["ahe_all"])
    out["core"] = {}
    for c, n in (("cpi_u", "CPI-U"), ("core_cpi", "CPI-U less food and energy"), ("pce_price", "PCE price index"),
                 ("core_pce", "PCE less food and energy")):
        p = pct(M[a][c], M[end][c])
        out["core"][c] = (p, real_pct(g, p))
        print(f"   {n:28s} inflation {p:7.3f}   real pay {real_pct(g, p):7.3f}")
    core_disagree = [t[0] for t in valid if t[3] is not None and (t[1] < 0) != (t[3] < 0)]
    out["core_disagree"] = core_disagree
    print(f"   start months where the CPI-U and core CPI disagree on the sign: {len(core_disagree)} "
          f"({core_disagree[0]} to {core_disagree[-1]})")

    # 6. seasonally adjusted pay over a not seasonally adjusted price index
    print("\n6. Seasonally adjusted pay deflated by the NOT seasonally adjusted CPI-U")
    diffs, diffs12 = [], []
    for i, d in enumerate(months):
        if i == 0:
            continue
        p = months[i - 1]
        r, q = M[d], M[p]
        if None in (r["ahe_all"], q["ahe_all"], r["cpi_u"], q["cpi_u"], r["cpi_u_nsa"], q["cpi_u_nsa"]):
            continue
        sa = pct(q["ahe_all"] / q["cpi_u"], r["ahe_all"] / r["cpi_u"])
        nsa = pct(q["ahe_all"] / q["cpi_u_nsa"], r["ahe_all"] / r["cpi_u_nsa"])
        diffs.append((d, nsa - sa))
    for d in months:
        prev = f"{int(d[:4]) - 1}{d[4:]}"
        if prev not in M:
            continue
        r, q = M[d], M[prev]
        if None in (r["ahe_all"], q["ahe_all"], r["cpi_u"], q["cpi_u"], r["cpi_u_nsa"], q["cpi_u_nsa"]):
            continue
        diffs12.append((d, pct(q["ahe_all"] / q["cpi_u_nsa"], r["ahe_all"] / r["cpi_u_nsa"])
                        - pct(q["ahe_all"] / q["cpi_u"], r["ahe_all"] / r["cpi_u"])))
    by_cal = {m: statistics.fmean(x for d, x in diffs if int(d[5:7]) == m) for m in range(1, 13)}
    big = max(diffs, key=lambda t: abs(t[1]))
    out.update(sa_n=len(diffs), sa_sd=statistics.stdev(x for _, x in diffs), sa_big=big, sa_by_cal=by_cal,
               sa12_n=len(diffs12), sa12_max=max(abs(x) for _, x in diffs12),
               sa12_mean=statistics.fmean(abs(x) for _, x in diffs12))
    print(f"   monthly real change, NSA-deflated minus SA-deflated: {len(diffs)} months, sd {out['sa_sd']:.3f} points, "
          f"largest {big[1]:+.3f} in {big[0]}")
    print("   average by calendar month: " + ", ".join(f"{datetime.date(2000, m, 1):%b} {v:+.3f}" for m, v in by_cal.items()))
    print(f"   12-month real change: {len(diffs12)} months, largest gap {out['sa12_max']:.3f}, "
          f"mean absolute gap {out['sa12_mean']:.3f} points")

    # 7. chained dollars do not add
    print("\n7. Real GDP against the sum of its chained-dollar components (C + I + G + NX), billions")
    qs = [d for d in sorted(Q) if Q[d]["net_exports_real"] is not None]
    resid = {d: Q[d]["gdp_real"] - (Q[d]["pce_real"] + Q[d]["investment_real"] + Q[d]["government_real"]
                                    + Q[d]["net_exports_real"]) for d in qs}
    nominal_gap = max(abs(Q[d]["gdp"] - (Q[d]["pce"] + Q[d]["investment"] + Q[d]["government"] + Q[d]["net_exports"]))
                      for d in Q)
    years_q: dict[int, list[str]] = {}
    for d in qs:
        years_q.setdefault(int(d[:4]), []).append(d)

    def annual(col, y):
        return statistics.fmean(Q[d][col] for d in years_q[y])

    out["annual_resid"] = {}
    for y in (1970, 1990, 2017, max(y for y, v in years_q.items() if len(v) == 4)):
        gdp_r = annual("gdp_real", y)
        res = statistics.fmean(resid[d] for d in years_q[y])
        out["annual_resid"][y] = (gdp_r, res, 100 * res / gdp_r)
        print(f"   {y}: real GDP {gdp_r:,.3f}, residual {res:+,.3f} ({100 * res / gdp_r:+.3f}% of real GDP)")
    ref_q = [resid[d] for d in years_q[2017]]
    out["ref_year_max"] = max(abs(x) for x in ref_q)
    print("   2017 by quarter: " + ", ".join(f"{x:+.3f}" for x in ref_q))
    lq = qs[-1]
    worst_q = max(qs, key=lambda d: abs(resid[d] / Q[d]["gdp_real"]))
    out.update(last_q=lq, last_resid=resid[lq], last_sum=Q[lq]["gdp_real"] - resid[lq], worst_q=worst_q,
               worst_q_pct=100 * resid[worst_q] / Q[worst_q]["gdp_real"], nominal_gap=nominal_gap, n_q=len(qs))
    print(f"   latest quarter {lq}: real GDP {Q[lq]['gdp_real']:,.3f}, components {out['last_sum']:,.3f}, "
          f"residual {resid[lq]:+,.3f}; largest residual share {out['worst_q_pct']:+.3f}% in {worst_q}; "
          f"current-dollar components miss nominal GDP by at most {nominal_gap:.3f}")
    print("   shares of GDP (percent), chained 2017 dollars against current dollars:")
    out["shares"] = {}
    for y in (1970, max(y for y, v in years_q.items() if len(v) == 4)):
        row = {}
        for nom, real in (("pce", "pce_real"), ("investment", "investment_real"), ("government", "government_real"),
                          ("net_exports", "net_exports_real")):
            row[nom] = (100 * annual(real, y) / annual("gdp_real", y), 100 * annual(nom, y) / annual("gdp", y))
        out["shares"][y] = row
        print(f"   {y}: " + "; ".join(f"{k} {v[0]:.2f} vs {v[1]:.2f}" for k, v in row.items()))
    # the implicit price deflator is nominal over real
    same = sum(1 for d in Q if Q[d]["gdp_deflator"] is not None and
               Decimal(str(100 * Q[d]["gdp"] / Q[d]["gdp_real"])).quantize(Decimal("0.001"), rounding=ROUND_HALF_UP)
               == Decimal(str(Q[d]["gdp_deflator"])).quantize(Decimal("0.001")))
    out.update(deflator_same=same, n_quarters=len(Q))
    print(f"   GDPDEF equals 100 x GDP / GDPC1 to three decimals in {same} of {len(Q)} quarters")

    # 8. the 2020 composition jump
    f, ap = M["2020-02-01"], M["2020-04-01"]
    out.update(comp_nominal=pct(f["ahe_all"], ap["ahe_all"]), comp_cpi=pct(f["cpi_u"], ap["cpi_u"]))
    out["comp_real"] = real_pct(out["comp_nominal"], out["comp_cpi"])
    print(f"\n8. Feb 2020 -> Apr 2020: pay ${f['ahe_all']:.2f} -> ${ap['ahe_all']:.2f} ({out['comp_nominal']:+.3f}%), "
          f"CPI-U {out['comp_cpi']:+.3f}%, 'real' {out['comp_real']:+.3f}%")
    gaps = [d for d in months if d <= end and M[d]["cpi_u"] is None and M[d]["ahe_all"] is not None]
    print("   months with pay but no CPI-U: " + ", ".join(
        f"{d} (PCE price index {M[d]['pce_price']})" for d in gaps))
    return out


# ------------------------------------------------------------------- figure ---
def figure(out: dict) -> None:
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.dates as mdates
    import matplotlib.pyplot as plt

    os.makedirs(IMG_DIR, exist_ok=True)
    paper, ink, soft, grid, base, shade = "#ffffff", "#1b2430", "#5e646b", "#e3e2de", "#9a9893", "#efeeea"
    c1, c2 = "#2a78d6", "#eb6834"  # fixed order: CPI-U, PCE price index
    plt.rcParams.update({"font.family": "DejaVu Sans", "font.size": 11, "axes.edgecolor": grid,
                         "axes.labelcolor": ink, "xtick.color": soft, "ytick.color": soft, "text.color": ink,
                         "axes.spines.top": False, "axes.spines.right": False})
    rows = [t for t in out["by_start"] if t[1] is not None]
    x = [datetime.date.fromisoformat(t[0]) for t in rows]
    fig, ax = plt.subplots(figsize=(8.8, 5.2), dpi=150, facecolor=paper)
    ax.set_facecolor(paper)
    for d in out["disagree"]:
        d0 = datetime.date.fromisoformat(d)
        d1 = datetime.date(d0.year + (d0.month == 12), d0.month % 12 + 1, 1)
        ax.axvspan(d0, d1, color=shade, linewidth=0, zorder=0)
    ax.grid(axis="y", color=grid, linewidth=0.8, zorder=1)
    ax.axhline(0, color=base, linewidth=1.0, zorder=2)
    ax.plot(x, [t[1] for t in rows], color=c1, linewidth=2, zorder=3, label="deflated by the CPI-U")
    ax.plot(x, [t[2] for t in rows], color=c2, linewidth=2, zorder=3, label="deflated by the PCE price index")
    for value in (rows[0][2], rows[0][1]):  # start values, left of the first point so no line runs through them
        ax.annotate(f"{value:.1f}%", (x[0], value), textcoords="offset points", xytext=(-7, 0), fontsize=10,
                    color=ink, ha="right", va="center")
    ax.annotate("shaded: start months where the two\ndeflators disagree on the sign", (datetime.date(2020, 11, 1), -2.2),
                xytext=(datetime.date(2014, 6, 1), -3.3), fontsize=9.5, color=soft, ha="center", va="center",
                arrowprops={"arrowstyle": "-", "color": soft, "linewidth": 0.8})
    ax.set_ylim(-4.6, 21.5)
    ax.set_xlim(datetime.date(2004, 6, 1), datetime.date(2026, 3, 1))
    ax.xaxis.set_major_locator(mdates.YearLocator(4))
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
    ax.set_ylabel("percent")
    ax.set_xlabel("start month")
    ax.set_title(f"Real change in average hourly earnings from each start month to {window_label(out['end'])}",
                 loc="left", fontsize=11.5)
    ax.legend(frameon=False, loc="upper right")
    fig.tight_layout()
    path = os.path.join(IMG_DIR, "real-pay-by-start-month.png")
    fig.savefig(path, facecolor=paper)
    plt.close(fig)
    print(f"\nwrote {path}")


if __name__ == "__main__":
    sys.stdout.reconfigure(encoding="utf-8")
    if "--download" in sys.argv:
        download()
    results = analysis(load(MONTHLY_CSV), load(GDP_CSV))
    if "--figures" in sys.argv:
        figure(results)