"""New York's extra unemployment is people out of work for six months or more: the analysis
behind a SparkyData article.

    py src/build_article.py the-long-unemployed

The question: New York City's unemployment rate is above the nation's (4.8% against 4.2% in
September 2026). Is that unusual, and whose unemployment makes up the difference?

Datasets (cards in datasets/):
- ``bls_labour``: the State's LAUS estimates for the city, monthly from January 1976, seasonally
  adjusted and not, and the national rate (BLS LNS14000000 and LNU04000000, via FRED): the
  headline and the fifty-year comparison.
- ``cps_basic``: the CPS monthly files from June 2024 to September 2026 (no October 2025 file),
  every state: the story's window is the latest twelve, the check June 2024 to May 2025, the
  anchors' windows July 2025 to June 2026 and calendar 2025. Who is unemployed in the city against
  the nation, how long they have been looking, and whether the city's mix of workers or higher
  rates within the same groups make the gap.
- ``acs1_national``: the ACS one-year national files, ages 16 to 29, for young adults without a
  degree, city against nation, in 2019 and 2024.
The anchors are BLS's alternative measures of labour underutilization for the city (U-1, U-3,
July 2025 to June 2026) and the State Comptroller's 2025 group rates, both reproduced from the
same monthly files.

Nothing here says why the gap exists: the decompositions are descriptive, one characteristic at
a time (symmetric Kitagawa weights) and all four at once (direct standardisation on the nation's
rates). Deterministic: no sampling.
"""

from __future__ import annotations

import numpy as np
import pandas as pd

from article_kit import Article, lines, r, strip
from nycdata import acs_national, bls, cps

SLUG = "the-long-unemployed"
SEED = 20261014
MONTHS = ["January", "February", "March", "April", "May", "June", "July", "August", "September", "October",
          "November", "December"]
CPS_COLUMNS = ["GESTFIPS", "GTCO", "HRHHID", "HRHHID2", "HRYEAR4", "HRMONTH", "PEEDUCA", "PEHSPNON", "PEMLR", "PESEX",
               "PRCITSHP", "PRTAGE", "PRUNEDUR", "PRUNTYPE", "PTDTRACE", "PWCMPWGT", "PESCHENR"]
BLS_ALT = ("bls_alternative_measures_states_2026_06", "bls_alternative_measures_states_2026-06.html")
OSC = ("osc_nyc_labor_force_2025", "osc_nyc_labor_force_2025.pdf")
NYBTN = ("nyc_comptroller_nybtn_117", "nyc_comptroller_nybtn_117_2026-09.html")
# The story's documents (v0.9.3; style guide, section 1): the State's September statement (the hook),
# its release schedule (the dated watch), the EDC's reading of the same months, and the State
# Comptroller's attributions, each quoted from its saved copy.
DOL_NYC = ("nysdol_nyc_labor_statistics_2026_10", "nysdol_nyc_labor_statistics_2026-10-18.html")
DOL_SCHEDULE = ("nysdol_press_release_schedule_2026", "nysdol_press_release_schedule_2026-10-18.html")
EDC = ("nycedc_economic_snapshot_2026_09", "nycedc_economic_snapshot_2026-09.pdf")
DURATIONS = [("lt5", "Under 5 weeks", 0, 4), ("w5_14", "5 to 14 weeks", 5, 14), ("w15_26", "15 to 26 weeks", 15, 26),
             ("w27", "27 weeks or more", 27, 10_000)]
REASONS = {1: "lost", 2: "lost", 3: "lost", 4: "left", 5: "reentrant", 6: "new"}
REASON_LABELS = {"lost": "Lost a job or a temporary job ended", "left": "Left a job", "reentrant": "Returned to the labour force",
                 "new": "Looking for a first job"}


# ------------------------------------------------------------------ fifty years of LAUS


def fifty_years_of_rates(A: Article) -> pd.DataFrame:
    nyc = bls.laus_nyc(adjusted=True).rename(columns={"unemployment_rate": "nyc"})
    us = bls.us_unemployment(adjusted=True).rename(columns={"unemployment_rate": "us"})
    m = nyc.merge(us, on=["year", "month"], how="inner").sort_values(["year", "month"]).reset_index(drop=True)
    m["gap"] = (m["nyc"] - m["us"]).round(1)
    last = m.iloc[-1]
    us_all = bls.us_unemployment(adjusted=True)
    us_last = us_all.iloc[-1]
    above = m["gap"] > 0
    annual = m.groupby("year").agg(nyc=("nyc", "mean"), us=("us", "mean"), gap=("gap", "mean"), months=("gap", "size")).reset_index()
    full = annual[annual["months"] >= 11]                                     # 2025 has 11 months (no October)
    narrow = m.iloc[:-1][m.iloc[:-1]["gap"] <= last["gap"]].iloc[-1]
    feb = m[(m["year"] == 2026) & (m["month"] == 2)].iloc[0]
    A.results["laus"] = {
        "first": f"{int(m.iloc[0]['year'])}-{int(m.iloc[0]['month']):02d}", "last": f"{int(last['year'])}-{int(last['month']):02d}",
        "latest": {"nyc": float(last["nyc"]), "us": float(last["us"]), "gap": float(last["gap"])},
        "us_latest": {"month": f"{int(us_last['year'])}-{int(us_last['month']):02d}", "rate": float(us_last["unemployment_rate"])},
        "months": int(len(m)), "months_above": int(above.sum()), "share_above": r(100 * above.mean(), 0),
        "median_gap": r(float(m["gap"].median()), 1),
        "latest_gap_percentile": r(100 * float((m["gap"] <= last["gap"]).mean()), 0),
        "last_as_narrow": f"{int(narrow['year'])}-{int(narrow['month']):02d}",
        "annual_gap": {str(int(y)): r(g, 1) for y, g in zip(full["year"], full["gap"])},
        "gap_2015_2019": {"lo": r(float(full[full["year"].between(2015, 2019)]["gap"].min()), 1),
                          "hi": r(float(full[full["year"].between(2015, 2019)]["gap"].max()), 1)},
        "years_gap_under_half_point": [int(y) for y, g in zip(full["year"], full["gap"]) if g < 0.5],
        "years_city_at_or_below": [int(y) for y, g in zip(full["year"], full["gap"]) if g <= 0.0],
        "peak_2026": {"nyc": float(feb["nyc"])},
    }
    sa = bls.laus_nyc(adjusted=True)
    f = sa[(sa["year"] == 2026) & (sa["month"] == 2)].iloc[0]
    l = sa.iloc[-1]
    A.results["laus"]["since_feb_2026"] = {"labor_force": int(l["labor_force"] - f["labor_force"]), "employed": int(l["employed"] - f["employed"]),
                                           "unemployed": int(l["unemployed"] - f["unemployed"])}
    Q = A.quotes()
    q = Q.quote(*NYBTN, "all of the improvement came from a dip in the labor force")
    qs = A.results.setdefault("quotes", {})
    qs["nybtn_labor_force"] = q
    qs["nybtn_august"] = Q.quote(*NYBTN, "NYC's unemployment rate declined from 5.0% to 4.8% in August")
    qs["dol_september"] = Q.quote(*DOL_NYC, "The city's seasonally adjusted unemployment rate was 4.8 percent in September, unchanged from August")
    qs["dol_state_september"] = Q.quote(*DOL_NYC, "New York State's rate was 4.3 percent in September 2026")
    qs["dol_schedule_september"] = Q.quote(*DOL_SCHEDULE, "September 2026 October 15, 2026")
    qs["dol_schedule_october"] = Q.quote(*DOL_SCHEDULE, "October 2026 November 19, 2026")
    qs["edc_sixth"] = Q.quote(*EDC, "the unemployment rate fell to 4.8%, its sixth consecutive monthly decline and well below "
                                    "the 5.7% rate at the beginning of the year")
    qs["edc_participation"] = Q.quote(*EDC, "labor force participation also continued to decline, reaching 62.1% in August")
    qs["edc_long_term"] = Q.quote(*EDC, "Long-Term Unemployment Rate (2026:Q2) 2.3%")
    qs["edc_cooling"] = Q.quote(*EDC, "suggesting a labor market that is cooling with relatively little churn")
    qs["osc_youth_black"] = Q.quote(*OSC, "Some demographic groups such as youth and Black individuals, though improved, continue to "
                                          "face elevated unemployment as they did prior to the pandemic")
    qs["osc_native_born"] = Q.quote(*OSC, "In 2025, a much higher share of native-born individuals was unemployed than before the "
                                          "pandemic, whereas the gap for foreign-born workers was significantly smaller")
    qs["osc_immigration"] = Q.quote(*OSC, "federal immigration policy restrictions before, during and after the pandemic have led to a "
                                          "decline in this population")
    qs["osc_college"] = Q.quote(*OSC, "college-educated individuals have yet to recover their pre-pandemic level of unemployment, "
                                      "perhaps due to limited entry")
    qs["osc_uneven"] = Q.quote(*OSC, "New York City continues to experience uneven outcomes across demographic groups, even as "
                                     "overall unemployment has improved since the onset of the pandemic")
    # The series' own landmarks, for the history (our figures, not a document's).
    pre = m[m["year"] < 2020]
    pk, pd_, lo = pre.loc[pre["nyc"].idxmax()], m.loc[m["nyc"].idxmax()], m.loc[m["nyc"].idxmin()]
    y76 = m[m["year"] == int(m.iloc[0]["year"])]
    declines = 0
    for k in range(len(m) - 1, 0, -1):
        if m.iloc[k]["nyc"] < m.iloc[k - 1]["nyc"]:
            declines += 1
        else:
            break
    jan = m[(m["year"] == int(last["year"])) & (m["month"] == 1)].iloc[0]
    A.results["laus"]["landmarks"] = {
        "first_year": {"year": int(y76.iloc[0]["year"]), "nyc": r(float(y76["nyc"].mean()), 1), "us": r(float(y76["us"].mean()), 1)},
        "pre_pandemic_peak": {"month": f"{int(pk['year'])}-{int(pk['month']):02d}", "nyc": float(pk["nyc"]), "us": float(pk["us"])},
        "pandemic_peak": {"month": f"{int(pd_['year'])}-{int(pd_['month']):02d}", "nyc": float(pd_["nyc"]), "us": float(pd_["us"])},
        "low": {"month": f"{int(lo['year'])}-{int(lo['month']):02d}", "nyc": float(lo["nyc"]), "us": float(lo["us"])},
        "consecutive_declines_to_latest": int(declines),
        "widest_gap": {"month": f"{int(m.loc[m['gap'].idxmax(), 'year'])}-{int(m.loc[m['gap'].idxmax(), 'month']):02d}",
                       "gap": float(m["gap"].max())},
        # the falls before the latest month (the EDC's "sixth consecutive monthly decline" to August)
        "declines_before_latest": int(_declines_ending_at(m, len(m) - 2)),
        "january_this_year": {"month": f"{int(jan['year'])}-01", "nyc": float(jan["nyc"])},
    }
    # The readings in the coverage, each figure read from its quote (never typed).
    import re as _re

    def qnum(key, pat):
        mm = _re.search(pat, qs[key]["quote"])
        if not mm:
            raise AssertionError(f"{key}: {pat!r} not in the quote")
        return float(mm.group(1))
    A.results["news"] = {
        "dol_september_rate": qnum("dol_september", r"was ([\d.]+) percent in September, unchanged"),
        "edc_august_rate": qnum("edc_sixth", r"fell to ([\d.]+)%, its sixth"),
        "edc_consecutive_declines": 6.0,
        "edc_start_of_year_rate": qnum("edc_sixth", r"below the ([\d.]+)% rate at the beginning"),
        "edc_participation": qnum("edc_participation", r"reaching ([\d.]+)% in August"),
        "edc_long_term_q2": qnum("edc_long_term", r"\(2026:Q2\) ([\d.]+)%"),
        "dol_state_september_rate": qnum("dol_state_september", r"rate was ([\d.]+) percent in September 2026"),
    }
    if L_last_is(m, 2026, 9):
        A.reproduce("The city's unemployment rate, September 2026, seasonally adjusted (NYS Department of Labor, %)",
                    published=A.results["news"]["dol_september_rate"], published_moe90=None, reproduced=float(last["nyc"]),
                    reproduced_moe90=None, kind="pct", tolerance=0.05)
    return m


def _declines_ending_at(m, k: int) -> int:
    """How many consecutive month-on-month falls in the city's rate end at row ``k``."""
    n = 0
    while k > 0 and m.iloc[k]["nyc"] < m.iloc[k - 1]["nyc"]:
        n += 1
        k -= 1
    return n


def L_last_is(m, year: int, month: int) -> bool:
    last = m.iloc[-1]
    return int(last["year"]) == year and int(last["month"]) == month


# ------------------------------------------------------------------ the CPS


def load_cps() -> tuple[pd.DataFrame, cps.Design, list, np.ndarray]:
    """Adults 16 and over in every state, January 2025 to the latest file; the story's window is
    the latest twelve files (``period``), the anchors' windows are earlier."""
    months = cps.months()
    d = cps.load(CPS_COLUMNS, months=months)
    d = d[(d["PRTAGE"] >= 16).to_numpy() & d["PEMLR"].between(1, 7).to_numpy()].reset_index(drop=True)
    latest = months[-12:]
    period = np.isin(cps.yyyymm(d), [y * 100 + m for y, m in latest])
    return d, cps.Design(d), latest, period


def rate(d, des, mask, y=None, weight="PWCMPWGT") -> tuple[float, float]:
    lf = (d["PEMLR"] <= 4).to_numpy()
    un = d["PEMLR"].isin([3, 4]).to_numpy() if y is None else y
    v, m = cps.ratio_moe(d, y=un.astype(float), x=lf.astype(float), weight=weight, mask=mask & lf, design=des)
    return 100 * v, 100 * m


def diff(d, des, mask_a, mask_b, y=None, weight="PWCMPWGT") -> tuple[float, float]:
    lf = (d["PEMLR"] <= 4).to_numpy().astype(float)
    un = (d["PEMLR"].isin([3, 4]).to_numpy() if y is None else y).astype(float)
    w = d[weight]
    ra, za = cps.ratio(un, lf, w, mask_a & (lf > 0))
    rb, zb = cps.ratio(un, lf, w, mask_b & (lf > 0))
    return 100 * (ra - rb), 100 * des.moe(za - zb)


def groupings(d) -> dict:
    """Label arrays for the decompositions (all CPS adults 16+); education for 25 and over."""
    age = d["PRTAGE"].astype(int).to_numpy()
    edu = d["PEEDUCA"].astype(int).to_numpy()
    race = np.select([d["PEHSPNON"] == 1, d["PTDTRACE"] == 1, d["PTDTRACE"] == 2, d["PTDTRACE"] == 4],
                     ["hispanic", "white", "black", "asian"], "other")
    return {
        "race_ethnicity": (race, None),
        "age": (np.select([age < 25, age < 35, age < 45, age < 55, age < 65], ["16_24", "25_34", "35_44", "45_54", "55_64"], "65_plus"), None),
        "education": (np.select([edu <= 38, edu == 39, edu <= 42, edu == 43], ["no_diploma", "high_school", "some_college", "bachelors"], "graduate"), age >= 25),
        "nativity": (np.select([d["PRCITSHP"] <= 3, d["PRCITSHP"] == 4], ["us_born", "naturalized"], "noncitizen"), None),
        "sex": (np.where(d["PESEX"] == 1, "men", "women"), None),
    }


GROUP_LABELS = {"hispanic": "Hispanic", "white": "White", "black": "Black", "asian": "Asian", "other": "Other",
                "16_24": "16 to 24", "25_34": "25 to 34", "35_44": "35 to 44", "45_54": "45 to 54", "55_64": "55 to 64", "65_plus": "65 and over",
                "no_diploma": "No high-school diploma (25+)", "high_school": "High school only (25+)", "some_college": "Some college (25+)",
                "bachelors": "Bachelor's degree (25+)", "graduate": "Graduate degree (25+)",
                "us_born": "Native-born", "naturalized": "Naturalized citizens", "noncitizen": "Not citizens", "men": "Men", "women": "Women"}


def kitagawa(d, des, nyc, labels, universe, period) -> dict:
    """The NYC-US gap = composition + rate effects, with symmetric (Kitagawa) weights; margins
    by linearisation over the design. The nation includes the city."""
    w = d["PWCMPWGT"].to_numpy(float)
    lf = (d["PEMLR"] <= 4).to_numpy(); un = d["PEMLR"].isin([3, 4]).to_numpy()
    base = lf & period & (np.ones(len(d), bool) if universe is None else universe)
    cats = sorted(np.unique(labels[base]))
    G = len(cats)
    gi = np.searchsorted(np.array(cats), labels)
    wu, wl = np.where(base & un, w, 0.0), np.where(base, w, 0.0)
    T = np.zeros(4 * G)
    for g in range(G):
        m = gi == g
        T[g], T[G + g], T[2 * G + g], T[3 * G + g] = wu[m & nyc].sum(), wl[m & nyc].sum(), wu[m].sum(), wl[m].sum()

    def parts(T):
        UN, LN, UU, LU = T[:G], T[G:2 * G], T[2 * G:3 * G], T[3 * G:]
        rN, rU = UN / LN, UU / LU
        sN, sU = LN / LN.sum(), LU / LU.sum()
        return np.array([UN.sum() / LN.sum() - UU.sum() / LU.sum(), ((sN - sU) * (rN + rU) / 2).sum(), ((rN - rU) * (sN + sU) / 2).sum()])

    f0 = parts(T)
    grad = np.zeros((3, 4 * G))
    for k in range(4 * G):
        h = max(T[k] * 1e-6, 1e-3)
        Tp = T.copy(); Tp[k] += h
        grad[:, k] = (parts(Tp) - f0) / h
    moe = [des.moe(grad[j][gi] * wu * nyc + grad[j][G + gi] * wl * nyc + grad[j][2 * G + gi] * wu + grad[j][3 * G + gi] * wl)
           for j in range(3)]
    rows = {}
    for g, c in enumerate(cats):
        mN, mU = base & nyc & (gi == g), base & (gi == g)
        rn, mrn = rate(d, des, mN)
        ru, mru = rate(d, des, mU)
        dd, md = diff(d, des, mN, mU)
        rows[c] = {"label": GROUP_LABELS.get(c, c), "nyc": r(rn, 1), "nyc_moe": r(mrn, 1), "us": r(ru, 1), "us_moe": r(mru, 1),
                   "diff": r(dd, 1), "lo": r(dd - md, 1), "hi": r(dd + md, 1),
                   "nyc_lf_share": r(100 * wl[mN].sum() / wl[base & nyc].sum(), 1), "us_lf_share": r(100 * wl[mU].sum() / wl[base].sum(), 1),
                   "records_lf": int(mN.sum()), "records_unemployed": int((mN & un).sum())}
    return {"gap": r(100 * f0[0], 2), "gap_moe": r(100 * moe[0], 2), "composition": r(100 * f0[1], 2), "composition_moe": r(100 * moe[1], 2),
            "rate": r(100 * f0[2], 2), "rate_moe": r(100 * moe[2], 2),
            "composition_lo": r(100 * (f0[1] - moe[1]), 2), "composition_hi": r(100 * (f0[1] + moe[1]), 2),
            "rate_lo": r(100 * (f0[2] - moe[2]), 2), "rate_hi": r(100 * (f0[2] + moe[2]), 2), "groups": rows}


def joint_mix(d, des, nyc, period, race_only: bool = False) -> dict:
    """Direct standardisation on race and ethnicity, education, age and citizenship at once: the
    city's rate if each of its cells had the nation's rate, minus the nation's rate, is what the
    city's mix of workers alone would do to the gap. Margin by linearisation over the design."""
    G = groupings(d)
    edu = np.where(d["PRTAGE"].to_numpy() >= 25, G["education"][0], "under_25")
    lab = np.char.add(np.char.add(np.char.add(G["race_ethnicity"][0].astype(str), "|"), np.char.add(edu.astype(str), "|")),
                      np.char.add(np.char.add(G["age"][0].astype(str), "|"), G["nativity"][0].astype(str)))
    if race_only:
        lab = G["race_ethnicity"][0].astype(str)
    w = d["PWCMPWGT"].to_numpy(float)
    lf = (d["PEMLR"] <= 4).to_numpy() & period
    un = d["PEMLR"].isin([3, 4]).to_numpy()
    cats, gi = np.unique(lab, return_inverse=True)
    C = len(cats)
    LN = np.bincount(gi, weights=np.where(lf & nyc, w, 0.0), minlength=C)
    LU = np.bincount(gi, weights=np.where(lf, w, 0.0), minlength=C)
    UU = np.bincount(gi, weights=np.where(lf & un, w, 0.0), minlength=C)
    UN = np.bincount(gi, weights=np.where(lf & un & nyc, w, 0.0), minlength=C)
    with np.errstate(divide="ignore", invalid="ignore"):
        rU = np.where(LU > 0, UU / LU, 0.0)
    std = (LN * rU).sum() / LN.sum()                  # the city's mix at the nation's rates
    us, city = UU.sum() / LU.sum(), UN.sum() / LN.sum()
    comp, within = std - us, city - std
    # influence of each record on comp = std - us (std depends on LN, UU, LU; us on UU, LU)
    lnN = LN.sum()
    z = np.zeros(len(d))
    m_n = lf & nyc
    z += np.where(m_n, w * (rU[gi] - std) / lnN, 0.0)                                        # through LN
    with np.errstate(divide="ignore", invalid="ignore"):
        dr_du = np.where(LU > 0, 1.0 / LU, 0.0); dr_dl = np.where(LU > 0, -UU / LU ** 2, 0.0)
    wn = LN / lnN
    z += np.where(lf & un, w * wn[gi] * dr_du[gi], 0.0) + np.where(lf, w * wn[gi] * dr_dl[gi], 0.0)   # through rU
    z -= np.where(lf & un, w / LU.sum(), 0.0) - np.where(lf, w * us / LU.sum(), 0.0)                 # minus the nation's rate
    moe = des.moe(z)
    return {"composition": r(100 * comp, 2), "lo": r(100 * (comp - moe), 2), "hi": r(100 * (comp + moe), 2),
            "within": r(100 * within, 2), "cells": int(C), "cells_city": int((LN > 0).sum())}


#: The valid entries of every CPS variable the analysis reads, as the 2024, 2025 and 2026 record
#: layouts give them (data/reference/cps/cps_basic_layout_<year>.txt; identical in all three).
CPS_VALID = {"PEMLR": range(1, 8), "PRUNEDUR": range(0, 120), "PRUNTYPE": range(1, 7), "PEEDUCA": range(31, 47),
             "PRCITSHP": range(1, 6), "PTDTRACE": range(1, 27), "PEHSPNON": range(1, 3), "PESCHENR": range(1, 3)}
#: The universe in which each is asked (the rest of the records carry a not-in-universe code).
CPS_UNIVERSE = {"PRUNEDUR": "unemployed", "PRUNTYPE": "unemployed", "PESCHENR": "16_54"}


def code_checks(d) -> dict:
    """Invalid entries per variable, file by file, within each variable's universe: all must be zero
    (the outside review of 18 October 2026, R4, asked for the count)."""
    un = d["PEMLR"].isin([3, 4]).to_numpy()
    ages = d["PRTAGE"].between(16, 54).to_numpy()
    ym = cps.yyyymm(d)
    out = {}
    for v, valid in CPS_VALID.items():
        scope = {"unemployed": un, "16_54": ages}.get(CPS_UNIVERSE.get(v, ""), np.ones(len(d), bool))
        bad = scope & ~d[v].isin(list(valid)).to_numpy()
        out[v] = {"invalid": int(bad.sum()), "records": int(scope.sum()), "files": int(len(np.unique(ym[scope])))}
    assert all(c["invalid"] == 0 for c in out.values()), f"CPS codes outside the layouts' valid entries: {out}"
    w = d["PWCMPWGT"].to_numpy()
    assert (w > 0).all(), "a composite weight of zero or less among adults 16 and over"
    return out


def cps_block(A: Article, d, des, months, period) -> None:
    A.results["code_checks"] = code_checks(d)
    nyc = cps.in_nyc(d) & period
    every = period
    lf = (d["PEMLR"] <= 4).to_numpy()
    un = d["PEMLR"].isin([3, 4]).to_numpy()
    rn, mn = rate(d, des, nyc)
    ru, mu = rate(d, des, every)
    g, mg = diff(d, des, nyc, every)
    A.results["cps"] = {"first": f"{months[0][0]}-{months[0][1]:02d}", "last": f"{months[-1][0]}-{months[-1][1]:02d}", "files": len(months),
                        "nyc": r(rn, 1), "nyc_moe": r(mn, 1), "us": r(ru, 1), "us_moe": r(mu, 2), "gap": r(g, 1), "gap_moe": r(mg, 1),
                        "gap_lo": r(g - mg, 1), "gap_hi": r(g + mg, 1),
                        "records_lf": int((nyc & lf).sum()), "records_unemployed": int((nyc & un).sum()),
                        "households": int(np.unique(des.cluster[nyc & lf]).size)}
    # by how long people have been looking: unemployed in each band per 100 in the labour force
    dur = d["PRUNEDUR"].to_numpy()
    rows = {}
    for key, lab, lo, hi in DURATIONS:
        y = un & (dur >= lo) & (dur <= hi)
        a, ma = rate(d, des, nyc, y)
        b, mb = rate(d, des, every, y)
        dd, md = diff(d, des, nyc, every, y)
        rows[key] = {"label": lab, "nyc": r(a, 2), "nyc_moe": r(ma, 2), "us": r(b, 2), "us_moe": r(mb, 2), "diff": r(dd, 2),
                     "lo": r(dd - md, 2), "hi": r(dd + md, 2), "records": int((nyc & y).sum())}
    A.results["duration"] = rows
    lt = un & (dur >= 27)
    sN, msN = cps.ratio_moe(d, y=lt.astype(float), x=un.astype(float), weight="PWCMPWGT", mask=nyc & un, design=des)
    sU, msU = cps.ratio_moe(d, y=lt.astype(float), x=un.astype(float), weight="PWCMPWGT", mask=every & un, design=des)
    A.results["long_term_share"] = {"nyc": r(100 * sN, 0), "nyc_moe": r(100 * msN, 0), "us": r(100 * sU, 0)}
    # the part of the gap from searches under 27 weeks (the complement of the long-term part)
    short = un & (dur >= 0) & (dur < 27)
    dd, md = diff(d, des, nyc, every, short)
    A.results["under_27_weeks_of_gap"] = {"diff": r(dd, 2), "lo": r(dd - md, 2), "hi": r(dd + md, 2)}
    # the year before (June 2024 to May 2025, twelve files, before the city's 2025-26 rise in unemployment)
    ym = cps.yyyymm(d)
    prior = (ym >= 202406) & (ym <= 202505)
    pn, pu = cps.in_nyc(d) & prior, prior
    g0, mg0 = diff(d, des, pn, pu)
    l0, ml0 = diff(d, des, pn, pu, un & (dur >= 27))
    A.results["prior_year"] = {"first": "2024-06", "last": "2025-05", "gap": r(g0, 1), "gap_lo": r(g0 - mg0, 1), "gap_hi": r(g0 + mg0, 1),
                               "long_term_diff": r(l0, 2), "lo": r(l0 - ml0, 2), "hi": r(l0 + ml0, 2)}
    # why they are looking
    reason = np.array([REASONS.get(int(v), "") for v in d["PRUNTYPE"].to_numpy()])
    rr = {}
    for k in ("lost", "left", "reentrant", "new"):
        y = un & (reason == k)
        dd, md = diff(d, des, nyc, every, y)
        a, _ = rate(d, des, nyc, y)
        b, _ = rate(d, des, every, y)
        rr[k] = {"label": REASON_LABELS[k], "nyc": r(a, 2), "us": r(b, 2), "diff": r(dd, 2), "lo": r(dd - md, 2), "hi": r(dd + md, 2)}
    A.results["reasons"] = rr
    # the city's mix of workers against higher rates within the same groups
    A.results["decomposition"] = {k: kitagawa(d, des, nyc, lab, uni, period) for k, (lab, uni) in groupings(d).items()}
    A.results["joint_mix"] = joint_mix(d, des, nyc, period)
    A.results["race_mix"] = joint_mix(d, des, nyc, period, race_only=True)
    # the eleven files of the latest twelve calendar months (October 2025 to September 2026; no October file)
    ym = cps.yyyymm(d)
    last = months[-1]
    eleven = (ym >= (last[0] - 1) * 100 + last[1] + 1) & (ym <= last[0] * 100 + last[1])
    e_n, e_u = cps.in_nyc(d) & eleven, eleven
    g1, mg1 = diff(d, des, e_n, e_u)
    l1, ml1 = diff(d, des, e_n, e_u, d["PEMLR"].isin([3, 4]).to_numpy() & (d["PRUNEDUR"].to_numpy() >= 27))
    A.results["eleven_files"] = {"files": int(len(set(ym[eleven].tolist()))), "gap": r(g1, 1), "gap_lo": r(g1 - mg1, 1), "gap_hi": r(g1 + mg1, 1),
                                 "long_term_diff": r(l1, 2), "lo": r(l1 - ml1, 2), "hi": r(l1 + ml1, 2)}
    # how many work: the employment-population ratio, 16 and over
    emp = d["PEMLR"].isin([1, 2]).to_numpy()
    eN, meN = cps.ratio_moe(d, y=emp.astype(float), x=np.ones(len(d)), weight="PWCMPWGT", mask=nyc, design=des)
    eU, meU = cps.ratio_moe(d, y=emp.astype(float), x=np.ones(len(d)), weight="PWCMPWGT", mask=every, design=des)
    pN, mpN = cps.ratio_moe(d, y=lf.astype(float), x=np.ones(len(d)), weight="PWCMPWGT", mask=nyc, design=des)
    pU, _ = cps.ratio_moe(d, y=lf.astype(float), x=np.ones(len(d)), weight="PWCMPWGT", mask=every, design=des)
    A.results["employment_ratio"] = {"nyc": r(100 * eN, 1), "nyc_moe": r(100 * meN, 1), "us": r(100 * eU, 1),
                                     "participation_nyc": r(100 * pN, 1), "participation_nyc_moe": r(100 * mpN, 1), "participation_us": r(100 * pU, 1)}
    # young people out of school without a degree, 22 to 27 as in the census check (the enrolment
    # question's universe is ages 16 to 54 in the 2024-2026 record layouts)
    age = d["PRTAGE"].to_numpy(); edu = d["PEEDUCA"].to_numpy()
    young = (age >= 22) & (age <= 27) & (d["PESCHENR"] == 2).to_numpy() & (edu < 43)
    a, ma = rate(d, des, nyc & young)
    b, mb = rate(d, des, every & young)
    A.results["young_cps"] = {"nyc": r(a, 1), "nyc_moe": r(ma, 1), "us": r(b, 1), "records_lf": int((nyc & young & lf).sum()),
                              "records_unemployed": int((nyc & young & un).sum())}


def reproductions(A: Article, d, des) -> None:
    """BLS's U-1 and U-3 for the city (July 2025 to June 2026, eleven months) and the State
    Comptroller's 2025 group rates, from the same monthly files."""
    Q = A.quotes()
    qs = A.results.setdefault("quotes", {})
    qs["bls_alt_nyc"] = Q.quote(*BLS_ALT, "New York City 3.1 2.5 5.5 6.0 6.9 10.2")
    qs["bls_alt_states"] = Q.quote(*BLS_ALT, "Total, all states 1.8 2.0 4.4 4.6 5.3 8.2")
    qs["bls_alt_next"] = Q.quote(*BLS_ALT, "covering the 12 months ending in September 2026 (excluding October 2025), is scheduled for October 23, 2026")
    ym = cps.yyyymm(d)
    nyc = cps.in_nyc(d)
    un = d["PEMLR"].isin([3, 4]).to_numpy()
    dur = d["PRUNEDUR"].to_numpy()
    win = (ym >= 202507) & (ym <= 202606)
    if win.any() and len(np.unique(ym[win])) == 11:
        u1, mu1 = rate(d, des, nyc & win, un & (dur >= 15))
        u3, mu3 = rate(d, des, nyc & win)
        u1us, mu1us = rate(d, des, win, un & (dur >= 15))
        A.reproduce("U-1, New York City, July 2025 to June 2026 (unemployed 15 weeks or longer, % of the labour force)", published=3.1,
                    published_moe90=None, reproduced=r(u1, 1), reproduced_moe90=r(mu1, 1), kind="pct")
        A.reproduce("U-3, New York City, July 2025 to June 2026 (the unemployment rate, %)", published=5.5, published_moe90=None,
                    reproduced=r(u3, 1), reproduced_moe90=r(mu3, 1), kind="pct")
        A.reproduce("U-1, all states, July 2025 to June 2026 (%)", published=1.8, published_moe90=None, reproduced=r(u1us, 1),
                    reproduced_moe90=r(mu1us, 2), kind="pct")
        A.results["bls_alt"] = {"u1_nyc": 3.1, "u1_states": 1.8, "u3_nyc": 5.5, "next_issue": "2026-10-23",
                                "ours_u1_nyc": r(u1, 1), "ours_u1_nyc_moe": r(mu1, 1)}
    cy = (ym >= 202501) & (ym <= 202512)
    if cy.any():
        qs["osc_black"] = Q.quote(*OSC, "Black unemployment stood over three times higher than for Whites at 8.9 percent in 2025")
        qs["osc_hispanic"] = Q.quote(*OSC, "Hispanic unemployment declined to 6.3 percent in 2025")
        qs["osc_youth"] = Q.quote(*OSC, "falling 1.3 points to 11.9 percent")
        black = (d["PTDTRACE"] == 2).to_numpy() & (d["PEHSPNON"] != 1).to_numpy()
        hisp = (d["PEHSPNON"] == 1).to_numpy()
        youth = d["PRTAGE"].between(16, 24).to_numpy()
        for lab, pub, mask in (("Black", 8.9, black), ("Hispanic", 6.3, hisp), ("Ages 16 to 24", 11.9, youth)):
            v, m = rate(d, des, nyc & cy & mask)
            A.reproduce(f"Unemployment rate, {lab}, New York City, 2025 (State Comptroller, %)", published=pub, published_moe90=None,
                        reproduced=r(v, 1), reproduced_moe90=r(m, 1), kind="pct")


# ------------------------------------------------------------------ young adults in the ACS


def young_acs(A: Article) -> None:
    """Ages 22 to 27, not enrolled, by degree: the unemployment rate, city against nation."""
    out, raw = {}, {}
    for year in (2019, 2024):
        s = acs_national.load(year)
        d = s.data
        R = s.estimator()
        nyc = s.in_nyc()
        age, esr = d["AGEP"].to_numpy(), d["ESR"].to_numpy()
        lf = np.isin(esr, [1, 2, 3]); un = esr == 3
        base = (age >= 22) & (age <= 27) & (d["SCH"].to_numpy() == 1) & lf
        grad = d["SCHL"].to_numpy() >= 21
        blk, raw[year] = {}, {}
        for k, m in (("no_degree", ~grad), ("degree", grad)):
            for where, wm in (("nyc", nyc), ("us", np.ones(len(d), bool))):
                e = R.share(base & m & wm, un)
                raw[year][f"{k}_{where}"] = (100 * e.value, 100 * e.moe)
                blk[f"{k}_{where}"] = r(100 * e.value, 1)
                blk[f"{k}_{where}_moe"] = r(100 * e.moe, 1)
                blk[f"{k}_{where}_records"] = int((base & m & wm).sum())
        # the city-nation difference for those without a degree, with a margin from the replicates
        a, b = base & ~grad & nyc, base & ~grad
        full = R.share(a, un).value - R.share(b, un).value
        rw = np.asarray(R.rw, dtype=float)
        ra = (rw[a] * un[a][:, None]).sum(0) / rw[a].sum(0)
        rb = (rw[b] * un[b][:, None]).sum(0) / rw[b].sum(0)
        moe = R.z * np.sqrt(R.factor * (((ra - rb) - full) ** 2).sum())
        raw[year]["no_degree_diff"] = (100 * full, 100 * moe)
        blk["no_degree_diff"] = r(100 * full, 1)
        blk["no_degree_lo"], blk["no_degree_hi"] = r(100 * (full - moe), 1), r(100 * (full + moe), 1)
        gq = acs_national.in_group_quarters(s)
        nat = d["NATIVITY"].to_numpy()
        for k, m in (("no_degree_household", ~gq), ("no_degree_native", nat == 1), ("no_degree_foreign", nat == 2)):
            for where, wm in (("nyc", nyc), ("us", np.ones(len(d), bool))):
                e = R.share(base & ~grad & m & wm, un)
                blk[f"{k}_{where}"] = r(100 * e.value, 1)
                blk[f"{k}_{where}_moe"] = r(100 * e.moe, 1)
                blk[f"{k}_{where}_records"] = int((base & ~grad & m & wm).sum())
        e = R.share(base & ~grad & nyc, gq)
        blk["no_degree_nyc_share_in_gq"] = r(100 * e.value, 1)
        out[str(year)] = blk
    # Changes from 2019 to 2024, the two years independent samples, from the unrounded estimates and
    # margins; only the results are rounded (outside review of 18 October 2026, R7).
    def change(key: str) -> dict:
        (v19, m19), (v24, m24) = raw[2019][key], raw[2024][key]
        ch, mm = v24 - v19, float(np.hypot(m19, m24))
        return {"pts": r(ch, 1), "lo": r(ch - mm, 1), "hi": r(ch + mm, 1)}
    out["widening"] = change("no_degree_diff")         # did the city's excess over the nation widen?
    out["degree_change"] = change("degree_nyc")        # young graduates in the city
    A.results["young_acs"] = out


# ------------------------------------------------------------------ main


def main(argv: list[str] | None = None) -> int:
    """Run the analysis and write every file (no options)."""
    A = Article(SLUG, __file__, seed=SEED)
    m = fifty_years_of_rates(A)
    d, des, months, period = load_cps()
    cps_block(A, d, des, months, period)
    reproductions(A, d, des)
    young_acs(A)
    L = A.results["laus"]
    A.latest(L["last"], "State LAUS for New York City, seasonally adjusted (September 2026); CPS monthly files to September 2026")
    A.long_view(int(L["first"][:4]), int(L["last"][:4]), "State LAUS for New York City and the BLS national rate, monthly, seasonally adjusted")

    # ---- printed figures (the prose's strings)
    P = A.printed
    C, D0, Y = A.results["cps"], A.results["duration"], A.results["young_acs"]
    P("laus.latest.nyc", f"{L['latest']['nyc']:.1f}%")
    P("laus.latest.us", f"{L['latest']['us']:.1f}%")
    P("laus.share_above", f"{L['share_above']:.0f}%")
    P("laus.median_gap", f"{L['median_gap']:.1f}")
    P("laus.gap_2015_2019.lo", f"{L['gap_2015_2019']['lo']:.1f}")
    P("laus.gap_2015_2019.hi", f"{L['gap_2015_2019']['hi']:.1f}")
    assert L["annual_gap"]["2024"] == L["annual_gap"]["2025"], "the page says the gap was the same in 2024 and 2025"
    P("laus.annual_gap.2024", f"{L['annual_gap']['2024']:.1f}")
    P("laus.latest.gap", f"{L['latest']['gap']:.1f}")
    ly, lm = int(L["last_as_narrow"][:4]), int(L["last_as_narrow"][5:])
    P("laus.last_as_narrow", f"{MONTHS[lm - 1]} {ly}")
    P("laus.last", f"{MONTHS[int(L['last'][5:]) - 1]} {L['last'][:4]}")
    P("laus.peak_2026.nyc", f"{L['peak_2026']['nyc']:.1f}%")
    for k in ("labor_force", "employed", "unemployed"):
        P(f"laus.since_feb_2026.{k}", f"{abs(L['since_feb_2026'][k]):,}")
    P("cps.nyc", f"{C['nyc']:.1f}%")
    P("cps.us", f"{C['us']:.1f}%")
    P("cps.gap", f"{C['gap']:.1f}")
    P("cps.gap_lo", f"{C['gap_lo']:.1f}")
    P("cps.gap_hi", f"{C['gap_hi']:.1f}")
    A.results["cps"]["records_lf_rounded"] = int(round(C["records_lf"], -2))
    P("cps.records_lf_rounded", f"{A.results['cps']['records_lf_rounded']:,}")
    P("cps.records_unemployed", f"{C['records_unemployed']:,}")
    P("duration.w27.nyc", f"{D0['w27']['nyc']:.1f}")
    P("duration.w27.us", f"{D0['w27']['us']:.1f}")
    P("duration.w27.diff", f"{D0['w27']['diff']:.2f}")
    P("duration.w27.lo", f"{D0['w27']['lo']:.2f}")
    P("duration.w27.hi", f"{D0['w27']['hi']:.2f}")
    P("long_term_share.nyc", f"{A.results['long_term_share']['nyc']:.0f}%")
    P("long_term_share.us", f"{A.results['long_term_share']['us']:.0f}%")
    B = A.results["bls_alt"]
    P("bls_alt.u1_nyc", f"{B['u1_nyc']:.1f}%")
    P("bls_alt.u1_states", f"{B['u1_states']:.1f}%")
    P("bls_alt.u3_nyc", f"{B['u3_nyc']:.1f}%")
    rep = {x["statistic"]: x for x in A.results["reproduction"]}
    osc = {k: rep[f"Unemployment rate, {k}, New York City, 2025 (State Comptroller, %)"] for k in ("Black", "Hispanic", "Ages 16 to 24")}
    A.results["osc"] = {k.split()[0].lower(): {"published": v["published"], "ours": v["reproduced"]} for k, v in osc.items()}
    for k in ("black", "hispanic", "ages"):
        P(f"osc.{k}.published", f"{A.results['osc'][k]['published']:.1f}%")
        P(f"osc.{k}.ours", f"{A.results['osc'][k]['ours']:.1f}%")
    G0 = A.results["decomposition"]
    E11 = A.results["eleven_files"]
    assert E11["files"] == 11, "the page says eleven files"
    P("eleven_files.gap", f"{E11['gap']:.1f}")
    P("eleven_files.long_term_diff", f"{E11['long_term_diff']:.1f}")
    P("race_mix.composition", f"{A.results['race_mix']['composition']:.1f}")
    J = A.results["joint_mix"]
    P("joint_mix.composition", f"{abs(J['composition']):.1f}")
    assert J["lo"] < J["hi"] and J["hi"] < C["gap_lo"], "the page says the mix explains none of the gap"
    PY = A.results["prior_year"]
    P("prior_year.gap", f"{PY['gap']:.1f}")
    assert PY["lo"] > 0 and PY["long_term_diff"] >= 0.8 * PY["gap"], "the page says the long searches were about all of the earlier gap"
    S27 = A.results["under_27_weeks_of_gap"]
    for k in ("diff", "lo", "hi"):
        P(f"under_27_weeks_of_gap.{k}", f"{S27[k]:+.2f}".replace("+", "") if k != "diff" else f"{S27[k]:.2f}")
    YC = A.results["young_cps"]
    P("young_cps.nyc", f"{YC['nyc']:.1f}%")
    P("young_cps.us", f"{YC['us']:.1f}%")
    P("young_cps.nyc_moe", f"{YC['nyc_moe']:.1f}")
    P("decomposition.nativity.groups.us_born.nyc", f"{G0['nativity']['groups']['us_born']['nyc']:.1f}%")
    P("decomposition.nativity.groups.us_born.us", f"{G0['nativity']['groups']['us_born']['us']:.1f}%")
    P("decomposition.race_ethnicity.groups.hispanic.nyc", f"{G0['race_ethnicity']['groups']['hispanic']['nyc']:.1f}%")
    P("decomposition.race_ethnicity.groups.hispanic.us", f"{G0['race_ethnicity']['groups']['hispanic']['us']:.1f}%")
    P("decomposition.education.groups.graduate.nyc", f"{G0['education']['groups']['graduate']['nyc']:.1f}%")
    P("decomposition.education.groups.graduate.us", f"{G0['education']['groups']['graduate']['us']:.1f}%")
    for y, k in (("2024", "no_degree_nyc"), ("2024", "no_degree_us"), ("2019", "no_degree_nyc"), ("2019", "no_degree_us"),
                 ("2024", "degree_nyc"), ("2024", "degree_us"), ("2019", "degree_nyc"), ("2024", "no_degree_household_nyc"),
                 ("2024", "no_degree_native_nyc")):
        P(f"young_acs.{y}.{k}", f"{Y[y][k]:.1f}%")
    E = A.results["employment_ratio"]
    P("employment_ratio.nyc", f"{E['nyc']:.1f}%")
    P("employment_ratio.us", f"{E['us']:.1f}%")
    assert round(E["participation_nyc"]) == round(E["participation_us"]), "the page gives one participation figure for both"
    P("employment_ratio.participation_nyc", f"{E['participation_nyc']:.0f}%")
    P("reasons.lost.diff", f"{A.results['reasons']['lost']['diff']:.1f}")
    P("reasons.reentrant.diff", f"{A.results['reasons']['reentrant']['diff']:.1f}")
    # the history's landmarks and the readings checked (v0.9.3)
    LM = L["landmarks"]
    P("laus.landmarks.first_year.nyc", f"{LM['first_year']['nyc']:.1f}%")
    P("laus.landmarks.pre_pandemic_peak.nyc", f"{LM['pre_pandemic_peak']['nyc']:.1f}%")
    P("laus.landmarks.pandemic_peak.nyc", f"{LM['pandemic_peak']['nyc']:.1f}%")
    P("laus.landmarks.pandemic_peak.us", f"{LM['pandemic_peak']['us']:.1f}%")
    P("laus.landmarks.low.nyc", f"{LM['low']['nyc']:.1f}%")
    P("laus.landmarks.january_this_year.nyc", f"{LM['january_this_year']['nyc']:.1f}%")
    P("laus.landmarks.widest_gap.gap", f"{LM['widest_gap']['gap']:.1f}")
    assert LM["widest_gap"]["month"] == LM["pandemic_peak"]["month"], "the page says the widest gap came with the pandemic peak"
    assert LM["declines_before_latest"] == 6, "the page says six falls from March to August"
    N = A.results["news"]
    P("news.dol_september_rate", f"{N['dol_september_rate']:.1f}%")
    P("news.edc_august_rate", f"{N['edc_august_rate']:.1f}%")
    P("news.edc_start_of_year_rate", f"{N['edc_start_of_year_rate']:.1f}%")
    P("news.edc_participation", f"{N['edc_participation']:.1f}%")
    P("news.edc_long_term_q2", f"{N['edc_long_term_q2']:.1f}%")
    # the beliefs: the young, non-citizens and graduates, with their ranges
    for a, b in (("age", "16_24"), ("nativity", "noncitizen"), ("nativity", "naturalized")):
        g = G0[a]["groups"][b]
        P(f"decomposition.{a}.groups.{b}.nyc", f"{g['nyc']:.1f}%")
        P(f"decomposition.{a}.groups.{b}.us", f"{g['us']:.1f}%")
    A.claim("young_acs.degree_change", "up")
    P("young_acs.degree_change.lo", f"{Y['degree_change']['lo']:.1f}")
    P("young_acs.degree_change.hi", f"{Y['degree_change']['hi']:.1f}")
    P("young_acs.degree_change.pts", f"{A.results['young_acs']['degree_change']['pts']:.1f}")
    P("young_acs.2019.degree_nyc", f"{Y['2019']['degree_nyc']:.1f}%")
    for k, stated in (("w27", "up"),):
        A.claim(f"duration.{k}", stated)
    A.results["cps"]["lo"], A.results["cps"]["hi"] = C["gap_lo"], C["gap_hi"]
    A.claim("cps", "up")
    A.results["young_acs"]["2024"]["lo"], A.results["young_acs"]["2024"]["hi"] = Y["2024"]["no_degree_lo"], Y["2024"]["no_degree_hi"]
    A.claim("young_acs.2024", "up")
    A.claim("young_acs.widening", "up")

    # ---- charts
    m = m.copy()
    m["key"] = m["year"].astype(str) + "-" + m["month"].map("{:02d}".format)
    # Fifty years of months is too dense to mark (and too long for a line's tooltip): the chart
    # draws each whole year's average of the monthly rates, and the current year's months so far as
    # points of their own, since a part year is never set against whole years (house rule; outside
    # review of 18 October 2026, R3). The text's monthly figures stay monthly.
    yr = m.groupby("year").agg(nyc=("nyc", "mean"), us=("us", "mean"), months=("nyc", "size"),
                               last=("month", "max")).reset_index()
    cur = int(yr["year"].iloc[-1])
    whole, part = yr[yr["year"] < cur], yr[yr["year"] == cur].iloc[0]
    assert int(part["last"]) < 12, "the current year is complete: draw it on the line"
    xs = [[str(int(y)), str(int(y)) + (" (no October)" if int(y) == 2025 else "")] for y in whole["year"]]
    xs.append([str(cur), f"{cur} (January to {MONTHS[int(part['last']) - 1]})"])
    A.results["laus"]["chart_years"] = {"first": int(yr["year"].iloc[0]), "last_whole": int(whole["year"].iloc[-1]),
                                        "current": cur, "current_months": int(part["months"]),
                                        "current_nyc": r(float(part["nyc"]), 1), "current_us": r(float(part["us"]), 1),
                                        "years_city_at_or_below": [int(y) for y, g in zip(whole["year"], whole["nyc"] - whole["us"]) if round(g, 1) <= 0]}
    so_far = f"January to {MONTHS[int(part['last']) - 1]} {cur}"
    A.chart("fifty_years", lines(
        title="New York's unemployment rate has run above the nation's in most years since 1976",
        subtitle=f"Unemployment rate, seasonally adjusted, yearly average of the monthly rates, New York City and the United "
                 f"States, 1976 to {cur - 1}, and {so_far} on its own",
        x=xs, series=[["nyc", "New York City", "New York City"], ["us", "United States", "United States"],
                      ["nyc_ytd", f"New York City, {so_far}", f"New York City, {cur} so far"],
                      ["us_ytd", f"United States, {so_far}", f"United States, {cur} so far"]],
        cells={"": {"nyc": {str(int(y)): round(float(v), 1) for y, v in zip(whole["year"], whole["nyc"])},
                    "us": {str(int(y)): round(float(v), 1) for y, v in zip(whole["year"], whole["us"])},
                    "nyc_ytd": {str(cur): round(float(part["nyc"]), 1)},
                    "us_ytd": {str(cur): round(float(part["us"]), 1)}}},
        format="pct", axis="Unemployment rate", labels="legend",
        x_ticks=[str(int(y)) for y in yr["year"] if int(y) % 10 == 0],
        x_tick_labels={str(int(y)): str(int(y)) for y in yr["year"]},
        note="The city's rate is the State Labor Department's model-based estimate (LAUS), revised each spring; the "
             "nation's is the BLS survey rate. Each whole year is the average of its seasonally adjusted monthly rates "
             f"(2025 without October, when the survey was not conducted); {so_far} is averaged on its own and drawn "
             "apart, since a part year isn't set against whole years. The monthly rates are in data.csv.",
        provenance="New York State Department of Labor, LAUS; U.S. Bureau of Labor Statistics"))
    D = A.results["duration"]
    A.chart("by_duration", strip(
        title="Most of the city's extra unemployment is people out of work for half a year or more",
        subtitle=f"Unemployed people per 100 in the labour force, by how long they have been looking, the twelve latest "
                 f"monthly files ({MONTHS[months[0][1] - 1]} {months[0][0]} to {MONTHS[months[-1][1] - 1]} {months[-1][0]}, no October)",
        x=[["us", "United States", 1], ["nyc", "New York City", 0]],
        series=[[k, D[k]["label"], D[k]["label"]] for k, *_ in DURATIONS],
        cells={"": {k: {"nyc": D[k]["nyc"], "nyc_moe": D[k]["nyc_moe"], "us": D[k]["us"], "us_moe": D[k]["us_moe"]} for k, *_ in DURATIONS}},
        format="dec2", axis="Unemployed per 100 in the labour force", order="series", value_x="nyc",
        value_label="the city's figure", moe_label="90% margin",
        note="The CPS monthly files pooled over twelve months (no October 2025 file), composite weight; the nation includes "
             "the city. Margins are clustered on the household and approximate.",
        provenance="U.S. Census Bureau and BLS, Current Population Survey microdata"))
    G = A.results["decomposition"]
    rows = [("nativity", "us_born"), ("nativity", "naturalized"), ("nativity", "noncitizen"), ("race_ethnicity", "white"),
            ("race_ethnicity", "black"), ("race_ethnicity", "hispanic"), ("race_ethnicity", "asian"), ("education", "no_diploma"),
            ("education", "high_school"), ("education", "some_college"), ("education", "bachelors"), ("education", "graduate"),
            ("age", "16_24"), ("age", "25_34"), ("age", "35_44"), ("age", "45_54"), ("age", "55_64"), ("age", "65_plus")]
    A.chart("within_groups", strip(
        title="Within most groups, the city's estimated rate is the higher one",
        subtitle="Unemployment rate by group, New York City and the United States, the twelve latest monthly files (September 2025 to September 2026, no October)",
        x=[["us", "United States", 1], ["nyc", "New York City", 0]],
        series=[[f"{a}_{b}", G[a]["groups"][b]["label"], G[a]["groups"][b]["label"]] for a, b in rows],
        cells={"": {f"{a}_{b}": {"nyc": G[a]["groups"][b]["nyc"], "nyc_moe": G[a]["groups"][b]["nyc_moe"],
                                  "us": G[a]["groups"][b]["us"], "us_moe": G[a]["groups"][b]["us_moe"]} for a, b in rows}},
        format="pct1", axis="Unemployment rate", order="series", value_x="nyc", value_label="the city's rate", moe_label="90% margin",
        note="Education is for people 25 and over. Race groups other than Hispanic are non-Hispanic. Margins are clustered "
             "on the household and approximate; groups with few unemployed people in the city sample have wide margins.",
        provenance="U.S. Census Bureau and BLS, Current Population Survey microdata"))
    Y = A.results["young_acs"]
    A.chart("young_without_degrees", strip(
        title="In 2024, young New Yorkers without a degree were out of work about twice as often as their peers nationally",
        subtitle="Unemployment rate of people aged 22 to 27 who are not in school, by degree, 2019 and 2024",
        x=[["us2024", "United States, 2024", 1], ["nyc2019", "New York City, 2019", 2], ["nyc2024", "New York City, 2024", 0]],
        series=[["no_degree", "Without a bachelor's degree", "No degree"], ["degree", "With a bachelor's degree", "Degree"]],
        cells={"": {k: {"us2024": Y["2024"][f"{k}_us"], "us2024_moe": Y["2024"][f"{k}_us_moe"],
                        "nyc2019": Y["2019"][f"{k}_nyc"], "nyc2019_moe": Y["2019"][f"{k}_nyc_moe"],
                        "nyc2024": Y["2024"][f"{k}_nyc"], "nyc2024_moe": Y["2024"][f"{k}_nyc_moe"]} for k in ("no_degree", "degree")}},
        format="pct1", axis="Unemployment rate", order="series", value_x="nyc2024", value_label="the city's 2024 rate",
        moe_label="90% margin",
        note="The ACS one-year national files; the city is taken from the same files. Margins from the replicate weights.",
        provenance="U.S. Census Bureau, American Community Survey microdata"))

    window = f"{A.results['cps']['first']} to {A.results['cps']['last']} (12 files)"
    rows_out = [{"table": "laus_monthly", "period": k, "group": "", "nyc": n, "us": u, "nyc_moe": None, "us_moe": None}
                for k, n, u in zip(m["key"], m["nyc"], m["us"])]
    rows_out += [{"table": "duration_per_100_in_labour_force", "period": window, "group": v["label"], "nyc": v["nyc"], "us": v["us"],
                  "nyc_moe": v["nyc_moe"], "us_moe": v["us_moe"]} for v in D.values()]
    rows_out += [{"table": "reason_per_100_in_labour_force", "period": window, "group": v["label"], "nyc": v["nyc"], "us": v["us"],
                  "nyc_moe": None, "us_moe": None} for v in A.results["reasons"].values()]
    rows_out += [{"table": f"unemployment_rate_by_{a}", "period": window, "group": v["label"], "nyc": v["nyc"], "us": v["us"],
                  "nyc_moe": v["nyc_moe"], "us_moe": v["us_moe"]} for a, blk in G.items() for v in blk["groups"].values()]
    rows_out += [{"table": "unemployment_rate_22_27_out_of_school", "period": y, "group": lab, "nyc": Y[y][f"{k}_nyc"], "us": Y[y][f"{k}_us"],
                  "nyc_moe": Y[y][f"{k}_nyc_moe"], "us_moe": Y[y][f"{k}_us_moe"]}
                 for y in ("2019", "2024") for k, lab in (("no_degree", "without a bachelor's degree"), ("degree", "with a bachelor's degree"),
                                                         ("no_degree_household", "without a degree, not in group quarters"),
                                                         ("no_degree_native", "without a degree, native-born"), ("no_degree_foreign", "without a degree, foreign-born"))]
    A.write(data=pd.DataFrame(rows_out))
    print(f"{SLUG}: wrote {A.out}")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())