"""Who gets the rent freeze: the Rent Guidelines Board's 2026 freeze set against the most
recent survey of the households it covers, and against every renewal guideline since 1968 —
every number and figure in the article, with the working.

    py src/build_article.py who-gets-the-rent-freeze      # this analysis, then viz.py, then the manifest

This file is the article's notebook. The article keeps to what was found; the how and the
why live here, next to the code that computes them.

======================================================================================
1. The question, and the data that answer it
======================================================================================

On 25 June 2026 the Rent Guidelines Board set 0% for one- and two-year renewals of
rent-stabilized leases starting 1 October 2026 to 30 September 2027 (Order #58). Who gets
the freeze, and what is it worth to them?

* **The orders** (administrative). Order #58 (the freeze) and Order #57 (3% one-year, 4.5%
  two-year, the year before), each read from its saved page; the board's chart of every
  order since 1968 (#1 to #58), each value used asserted against the saved PDF's text.
* **The households** (survey). The 2023 New York City Housing and Vacancy Survey public use
  files (HPD and the Census Bureau; fieldwork January to June 2023; incomes are for 2022,
  rents for 2023). Regulation status is the survey's CSR variable, which HPD builds from
  state, NYCHA and HPD records with the respondent's answers. It is the latest survey: the
  next one has not been published.
* **Published anchors**: HPD's 2023 Selected Initial Findings and the board's 2024 Income
  and Affordability Study, each figure checked against the saved PDF before it is compared.
* **The long view** (survey). Every Housing and Vacancy Survey from 1993 to 2023, read by
  nycdata.hvs_income (card nychvs_income): the share of stabilized and market-rate renter
  households with incomes of $100,000 or more in 2022 dollars, their medians, and six-figure
  households' share of the stabilized rent, which is how a freeze of one rate on every lease
  would split its dollars. Incomes are for the year before each survey and are converted with
  the New York metro CPI-U for all items. The series starts in 1993, because the 1991 file
  leaves a quarter of renter incomes unreported where later files allocate them. It breaks
  twice: at 2017, when the survey began using state records of about 62,000 apartments
  registered as permanently exempt from stabilization (HPD's 2017 report: compare the
  stabilized population across it "with caution"), and at the 2021 redesign (new income
  questions, more records behind stabilization status). 1993-2014, 2017 and 2021-2023 are
  separate stretches, and a change is named only within one.
  The reader's checks (Census Series IA Table 9 for 1999-2014 and HPD's medians from 1993 to
  2023) must all pass before the numbers are used.

Loaders and estimators: src/nycdata/hvs.py for 2023 (the idea's validity scan reads the same
module) and src/nycdata/hvs_income.py for the long view, both built on the replicate-weight
estimator in shared/python/sparkylib/replicates.py; all three are published beside this file.
Raw files: data/raw/housing/nychvs_<year>/ (sources nychvs_*_puf_*); they are not shipped,
being the Census Bureau's files, and the script downloads nothing at run time.

======================================================================================
2. Definitions
======================================================================================

* **Rent stabilized**: CSR 32. Market rental: CSR 80. Public housing: CSR 5. Renter
  households: CSR 5, 32, 80, 90 and 97 (the survey's renter tenure).
* **Income**: the household's 2022 income as the survey records it (HHINC_REC1), in the
  dollars reported. Bands: under $25k, $25-50k, $50-100k, $100-200k, $200k and over.
  Citywide fifths are the 2022 income quintiles of every household in the city, renters
  and owners.
* **Low income**: at or below 80% of HUD's FY2023 income limit for the household's size
  (the survey's HUDILFY23, income as a percentage of the 100% limit); zero or negative
  incomes count as low income.
* **Rent**: contract rent (RENT_AMOUNT), the rent due each month "regardless of whether this
  amount was paid by the tenant(s) or others on behalf of the tenant" (codebook p. 126), the
  base a renewal increase applies to. Gross rent adds utilities. So a saving computed on it is
  a reduction in contract rent, not necessarily money the household keeps.
* **Rent burden**: contract rent over 30% (over 50%: severe) of income, among households
  outside means-tested housing (RENTBURDEN_CAT 1-3), the board's convention.

======================================================================================
3. The freeze's value (computed)
======================================================================================

What a household does not pay because of the freeze, against a stated counterfactual:
Order #57 repeated (3% for a one-year lease, 4.5% for a two-year lease, by the household's
LEASE_LENGTH; a tenancy with no active lease or none reported gets the one-year rate). The
annual saving at renewal is 12 x contract rent x that rate. At 2023 rents it is measured
from the survey; "rolled forward" multiplies it by the three one-year renewals since the
survey (Orders #55 to #57: 3%, 2.75%, 3%), an assumption that every tenancy renewed at the
one-year guideline and saw no other increase (no vacancy, improvement or MCI increases).
Scaling both rates together scales every household alike and leaves the shares by income unchanged; one
rate for both lease lengths moves them a little, because two-year leases are less common at six figures.

Timing: a one-year lease renews every order year, so every one-year tenant renews in the
freeze year; a two-year lease renews every other year, phased by LEASE_START (unknown
phases count one half). "Once every lease has renewed" sums every household; "in the
freeze year" sums those whose renewal falls between October 2026 and September 2027.

======================================================================================
4. Uncertainty
======================================================================================

Every estimate uses the final weight FW; every margin is a 90% margin from the 80
successive-difference replicate weights (HPD's 2023 Guide to Estimating Variances, eq.
2.1: Var = 4/80 x sum of squared replicate deviations; MOE = 1.645 x SE), medians
recomputed on each replicate. A difference between groups is a sign claim only when its
replicate-based interval excludes zero. Cells under 100 records are flagged and Staten
Island (27 stabilized records) is not shown. The computed figures carry the survey's
sampling error; the counterfactual and the rolled-forward rents are assumptions, swept
below, not estimates.

Outputs (all regenerable; nothing is edited by hand), in output/articles/who_gets_the_rent_freeze/
and served from site/public/articles/who-gets-the-rent-freeze/:
    results.json        every number quoted in the prose
    charts.json         the cubes behind the bar charts (Figures 2, 5, 6 and 7)
    fig*.svg            their static twins, for readers without scripts
    viz.json            the template charts (Figures 1, 3 and 4), written by viz.py
    data.csv            the 2023 renter households, with the variables used
    long_view.csv       the thirty-year series, one row per survey and group
    hvs.py, hvs_income.py, replicates.py
                        the loaders and the estimator this file imports
    build.json          version and checksums of every served file and library module (the build)
The build also writes who-gets-the-rent-freeze-code.zip (this package with the library modules
it imports) to output/ only, until code bundles are served.
"""

from __future__ import annotations

import argparse
import json
import re
from pathlib import Path

import numpy as np
import pandas as pd

from _paths import REFERENCE
from nycdata import hvs as H
from nycdata import bls as BLS
from nycdata import hvs_income as HI
from nycdata import recode
import article_kit as K
from sparkylib import replicates
from sparkylib.quotes import QuoteBook

SLUG = "who-gets-the-rent-freeze"
PATHS = K.article_paths(SLUG)          # output/articles/who_gets_the_rent_freeze/ and the site folder
OUT = PATHS.out
REPORTS = REFERENCE / "reports"

ORDER57 = {"one_year": 0.03, "two_year": 0.045}   # leases 1 Oct 2025 - 30 Sep 2026
ORDER58 = {"one_year": 0.0, "two_year": 0.0}      # leases 1 Oct 2026 - 30 Sep 2027: the freeze
FREEZE_ORDER_YEAR = 2026
SWEEP = [0.02, 0.025, 0.03, 0.035, 0.04, 0.045]
INC_BANDS = [("lt25k", "Under $25k", -np.inf, 25_000), ("25k_50k", "$25k–50k", 25_000, 50_000),
             ("50k_100k", "$50k–100k", 50_000, 100_000), ("100k_200k", "$100k–200k", 100_000, 200_000),
             ("200k_plus", "$200k and over", 200_000, np.inf)]
THIN = 100

# Every renewal guideline since 1968, from the board's chart of Orders #1 to #58
# (reference/reports/rgb_apartment_orders_chart.pdf). One row per guideline year: the year
# its leases start, the order, the date range as the chart prints it (with the occurrence
# to use where two rows share one), and the one- and two-year renewal guidelines in percent.
# Where an order splits a rate, the value used is stated: heat-split orders (1980, 1981,
# 2004-2006, 2008, 2009) take the rate where the owner provides heat; 1968 covered two
# years (to June 1970); 1980 covered fifteen months; early orders added a stabilizer or
# fuel adjustments on top; 2012 was "2% or $20, whichever is greater"; 2021's one-year
# lease rose 1.5% halfway through (0% for its first six months); 2020's two-year lease was
# 0% for its first year and 1% for its second; 2023's two-year lease was 2.75% and then 3.2%.
ORDERS = [
    (1968, "1", "7/1/68 to 6/30/70", 1, 10.0, 10.0), (1970, "2", "7/1/70 to 6/30/71", 1, 6.0, 8.0),
    (1971, "3", "7/1/71 to 6/30/72", 1, 7.0, 9.0), (1972, "4", "7/1/72 to 6/30/73", 1, 6.0, 8.0),
    (1973, "5", "7/1/73 to 6/30/74", 1, 6.5, 8.5), (1974, "6", "7/1/74 to 6/30/75", 1, 8.5, 10.5),
    (1975, "7", "7/1/75 to 6/30/76", 1, 7.5, 9.5), (1976, "8", "7/1/76 to 6/30/77", 1, 6.5, 8.0),
    (1977, "9", "7/1/77 to 6/30/78", 1, 6.5, 8.5), (1978, "10", "7/1/78 to 6/30/79", 2, 4.5, 6.5),
    (1979, "11", "7/1/79 to 6/30/80", 1, 8.5, 12.0), (1980, "12", "7/1/80 to 6/30/81", 1, 11.0, 14.0),
    (1981, "13", "10/1/81 to 9/30/82", 1, 10.0, 13.0), (1982, "14", "10/1/82 to 9/30/83", 1, 4.0, 7.0),
    (1983, "15", "10/1/83 to 9/30/84", 1, 4.0, 7.0), (1984, "16", "10/1/84 to 9/30/85", 1, 6.0, 9.0),
    (1985, "17", "10/1/85 to 9/30/86", 1, 4.0, 6.5), (1986, "18", "10/1/86 to 9/30/87", 1, 6.0, 9.0),
    (1987, "19", "10/1/87 to 9/30/88", 1, 3.0, 6.5), (1988, "20", "10/1/88 to 9/30/89", 1, 6.0, 9.0),
    (1989, "21", "10/1/89 to 9/30/90", 1, 5.5, 9.0), (1990, "22", "10/1/90 to 9/30/91", 1, 4.5, 7.0),
    (1991, "23", "10/1/91 to 9/30/92", 1, 4.0, 6.5), (1992, "24", "10/1/92 to 9/30/93", 1, 3.0, 5.0),
    (1993, "25", "10/1/93 to 9/30/94", 1, 3.0, 5.0), (1994, "26", "10/1/94 to 9/30/95", 1, 2.0, 4.0),
    (1995, "27", "10/1/95 to 9/30/96", 1, 2.0, 4.0), (1996, "28", "10/1/96 to 9/30/97", 1, 5.0, 7.0),
    (1997, "29", "10/1/97 to 9/30/98", 1, 2.0, 4.0), (1998, "30", "10/1/98 to 9/30/99", 1, 2.0, 4.0),
    (1999, "31", "10/1/99 to 9/30/00", 1, 2.0, 4.0), (2000, "32", "10/1/00 to 9/30/01", 1, 4.0, 6.0),
    (2001, "33", "10/1/01 to 9/30/02", 1, 4.0, 6.0), (2002, "34", "10/1/02 to 9/30/03", 1, 2.0, 4.0),
    (2003, "35", "10/1/03 to 9/30/04", 1, 4.5, 7.5), (2004, "36", "10/1/04 to 9/30/05", 1, 3.5, 6.5),
    (2005, "37", "10/1/05 to 9/30/06", 1, 2.75, 5.5), (2006, "38", "10/1/06 to 9/30/07", 1, 4.25, 7.25),
    (2007, "39", "10/1/07 to 9/30/08", 1, 3.0, 5.75), (2008, "40", "10/1/08 to 9/30/09", 1, 4.5, 8.5),
    (2009, "41", "10/1/09 to 9/30/10", 1, 3.0, 6.0), (2010, "42", "10/1/10 to 9/30/11", 1, 2.25, 4.5),
    (2011, "43", "10/1/11 to 9/30/12", 1, 3.75, 7.25), (2012, "44", "10/1/12 to 9/30/13", 1, 2.0, 4.0),
    (2013, "45", "10/1/13 to 9/30/14", 1, 4.0, 7.75), (2014, "46", "10/1/14 to 9/30/15", 1, 1.0, 2.75),
    (2015, "47", "10/1/15 to 9/30/16", 1, 0.0, 2.0), (2016, "48", "10/1/16 to 9/30/17", 1, 0.0, 2.0),
    (2017, "49", "10/1/17 to 9/30/18", 1, 1.25, 2.0), (2018, "50", "10/1/18 to 9/30/19", 1, 1.5, 2.5),
    (2019, "51", "10/1/19 to 9/30/20", 1, 1.5, 2.5), (2020, "52", "10/1/20 to 9/30/21", 1, 0.0, 0.0),
    (2021, "53", "10/1/21 to 9/30/22", 1, 1.5, 2.5), (2022, "54", "10/1/22 to 9/30/23", 1, 3.25, 5.0),
    (2023, "55", "10/1/23 to 9/30/24", 1, 3.0, 2.75), (2024, "56", "10/1/24 to 9/30/25", 1, 2.75, 5.25),
    (2025, "57", "10/1/25 to 9/30/26", 1, 3.0, 4.5), (2026, "58", "10/1/26 to 9/30/27", 1, 0.0, 0.0),
]
# Two-year leases frozen for their first year only: the second year rose (2020: 1% in the second year).
TWO_YEAR_FIRST_YEAR_ONLY = {2020}
# The second-year rate of the two orders that split a two-year lease (the table above holds the first-year rate):
# 2020, 0% then 1%; 2023, 2.75% then 3.2%. Each is asserted against its row of the board's chart.
TWO_YEAR_SECOND_YEAR = {2020: 1.0, 2023: 3.2}

res: dict = {}


def pct_of(r: dict, nd=1) -> dict:
    """A survey share as percentages: estimate, 90% margin, records."""
    return {"pct": K.r(100 * r["est"], nd), "moe": K.r(100 * r["moe90"], nd), "n": r.get("n")}


def money_of(r: dict) -> dict:
    return {"usd": K.r(r["est"], 0), "moe": K.r(r["moe90"], 0), "n": r.get("n")}


def count_of(r: dict) -> dict:
    return {"count": K.r(r["est"], 0), "moe": K.r(r["moe90"], 0), "n": r.get("n")}


# ====================================================================================
# Saved sources: every quoted figure is found in its file before it is used
# ====================================================================================

BOOK = QuoteBook(REPORTS, context=100)    # every quoted figure is found in its saved file first
quote = BOOK.quote


def quoted_figures() -> dict:
    """The orders and the published anchors, each verified in its saved file."""
    q = {}
    o58, o57 = "rgb_adopted_guidelines_2026_27.html", "rgb_order_57_2025_26.html"
    q["orders"] = [
        quote("rgb_adopted_guidelines_2026_27", o58, "For a one -year lease commencing on or after October 1, 2026 and on or before September 30, 2027 : 0%"),
        quote("rgb_adopted_guidelines_2026_27", o58, "For a two -year lease commencing on or after October 1, 2026 and on or before September 30, 2027 : 0%"),
        quote("rgb_order_57_2025_26", o57, "For a one -year lease commencing on or after October 1, 2025 and on or before September 30, 2026 : 3%"),
        quote("rgb_order_57_2025_26", o57, "For a two -year lease commencing on or after October 1, 2025 and on or before September 30, 2026 : 4.5%"),
        quote("rgb_adopted_guidelines_2026_27", o58, "RENT STABILIZED LOFT GUIDELINES For one -year increase periods commencing on or after October 1, 2026 and on or before September 30, 2027 : 0%"),
        quote("rgb_adopted_guidelines_2026_27", o58, "Residential Class A (apartment) hotels – 0% 2) Lodging houses – 0% 3) Rooming houses"),
        # Hotels are in a separate order adopted the same day; the guideline sits on top of other lawful increases.
        quote("rgb_adopted_guidelines_2026_27", o58, "Apartment and Loft Order 58 and Hotel Order 56"),
        quote("rgb_adopted_guidelines_2026_27", o58, "Together with such further adjustments as may be authorized by law, the annual adjustment for leases for apartments shall be"),
        # The board's preliminary ranges (7 May 2026), the benchmark set beside Order #57 repeated.
        quote("rgb_proposed_guidelines_2026_27", "rgb_proposed_guidelines_2026_27.html", "promulgated by the NYC Rent Guidelines Board on May 7, 2026"),
        quote("rgb_proposed_guidelines_2026_27", "rgb_proposed_guidelines_2026_27.html",
              "For a one -year lease commencing on or after October 1, 2026 and on or before September 30, 2027 : 0%-2%"),
        quote("rgb_proposed_guidelines_2026_27", "rgb_proposed_guidelines_2026_27.html",
              "For a two -year lease commencing on or after October 1, 2026 and on or before September 30, 2027 : 0%-4%"),
        # The owners' challenge, as reported: the fact of the suit and that the freeze took effect while it is pending.
        quote("ny1_rgb_lawsuit_2026_07_23", "ny1_rgb_lawsuit_2026_07_23.html", "a group of building owners sued the city's Rent Guidelines Board"),
        quote("ny1_rgb_lawsuit_2026_07_23", "ny1_rgb_lawsuit_2026_07_23.html", "The landlords in the lawsuit want a judge to void the freeze. They want the board to meet again"),
        quote("ny1_rent_freeze_in_effect_2026_10_01", "ny1_rent_freeze_in_effect_2026_10_01.html", "The rent freeze for rent-stabilized apartments takes effect Thursday"),
        quote("ny1_rent_freeze_in_effect_2026_10_01", "ny1_rent_freeze_in_effect_2026_10_01.html", "The order stemmed from the lawsuit seeking to overturn the board"),
        quote("rgb_apartment_orders_chart", "rgb_apartment_orders_chart.pdf",
              "A separate vacancy allowance is not permitted under the Housing Stability and Tenant Protection Act of 2019", 1),
        quote("rgb_apartment_orders_chart", "rgb_apartment_orders_chart.pdf",
              "lease adjustments for one- and two-year renewal leases also apply to vacant apartments that become occupied during the term of the Order", 1),
    ]
    sif, ia = "2023-nychvs-selected-initial-findings.pdf", "nyc_rgb_income_affordability_study_2024.pdf"
    q["sif"] = [
        quote("nychvs_2023_initial_findings", sif, "The net rental vacancy rate in 2023 for all housing accommodations in New York City was 1.41 percent"),
        quote("nychvs_2023_initial_findings", sif, "In 2023, there were 2,324,000 renter-occupied households"),
        quote("nychvs_2023_initial_findings", sif, "Rent Stabilized 26% ±2% 34% ±2% 25% ±1% 15% ±1% 960,700"),
        quote("nychvs_2023_initial_findings", sif, "Market Rental 16% ±1% 19% ±1% 26% ±2% 40% ±2% 1,119,000"),
        quote("nychvs_2023_initial_findings", sif, "Public Housing 77% ±3% 16% ±2% 5% ±1% ** ** 167,700"),
        quote("nychvs_2023_initial_findings", sif, "In 2023, there were 996,600 rent stabilized units"),
        quote("nychvs_2023_initial_findings", sif, "Rent Stabilized $1,400 $1,542 $1,500"),
        quote("nychvs_2023_initial_findings", sif, "Market Rental $1,825 $2,010 $2,000"),
        quote("nychvs_2023_initial_findings", sif, "market renters had a median household income of $90,800"),
        quote("nychvs_2023_initial_findings", sif, "Public housing residents had the lowest household incomes with a median of $20,600"),
    ]
    q["ia"] = [
        quote("nyc_rgb_income_affordability_2024", ia, "median household income was $60,000, median contract rent was $1,500, and median gross rent was $1,570"),
        quote("nyc_rgb_income_affordability_2024", ia, "and 30% made $100,000 or more (as compared to 47% of market rentals)"),
        quote("nyc_rgb_income_affordability_2024", ia, "had a median gross rent-to-income ratio of 28.8%"),
        quote("nyc_rgb_income_affordability_2024", ia, "45.5%.41 This includes 18.3% considered moderately rent burdened and 27.2% considered severely rent burdened"),
    ]
    # The next survey: run with the University of Michigan, results not yet published (the page's last check, 2 October 2026).
    q["surveys"] = [
        quote("hpd_nychvs_2026_release", "hpd_nychvs_2026_release.html", "Institute for Social Research (ISR) to conduct the 2026 NYCHVS"),
        # The long view's two breaks: the 2017 recode of stabilization and the 2021 redesign.
        quote("nychvs_2017_initial_findings", "2017-hvs-initial-findings.pdf", "additional information on 62,000 units that registered with DHCR as permanently exempt"),
        quote("nychvs_2017_initial_findings", "2017-hvs-initial-findings.pdf",
              "any comparisons of the rent stabilized stock or tenant population between 2017 and earlier HVSs should be made with caution"),
        quote("nychvs_2021_initial_findings", "2021-nychvs-selected-initial-findings.pdf",
              "Questions related to income were substantially restructured in their phrasing and administration"),
    ]
    return {"quotes": q,
            "anchors": {"vacancy_pct": 1.41, "renter_households": 2_324_000, "stabilized_households": 960_700,
                        "market_households": 1_119_000, "public_housing_households": 167_700, "stabilized_units": 996_600,
                        "stabilized_median_contract_rent": 1500, "market_median_contract_rent": 2000,
                        "stabilized_median_income": 60_000, "market_median_income": 90_800, "public_housing_median_income": 20_600,
                        "stabilized_100k_pct": 30, "market_100k_pct": 47, "stabilized_median_gross_rent": 1570,
                        "stabilized_no_assistance_gross_rti_pct": 28.8, "stabilized_no_assistance_burdened_pct": 45.5,
                        "stabilized_no_assistance_severe_pct": 27.2}}


def orders_history() -> dict:
    """The renewal guidelines since 1968, each value found in its row of the board's chart."""
    pages = BOOK.pages("rgb_apartment_orders_chart.pdf")
    text = " ".join(pages)
    starts = [(m.start(), m.group(0)) for m in re.finditer(r"\d{1,2}/\d{1,2}/\d{2} to \d{1,2}/\d{1,2}/\d{2}", text)]
    rows = []
    for year, order, dates, nth, one, two in ORDERS:
        hits = [i for i, (_, s) in enumerate(starts) if s == dates]
        assert len(hits) >= nth, f"Order {order}: {dates!r} not found {nth} time(s) in the chart"
        i = hits[nth - 1]
        seg = text[starts[i][0]: starts[i + 1][0] if i + 1 < len(starts) else len(text)]
        vals = [float(x) for x in re.findall(r"(\d+(?:\.\d+)?)\s?%", seg)]
        assert one in vals and two in vals, f"Order {order} ({dates}): {one}% / {two}% not both in its row: {vals}"
        # Columns: every row prints its one-year rate before its two-year rate, except 2021's, whose
        # one-year cell ("0% for six months, then 1.5%") wraps around the two-year 2.5%.
        if one != two and year != 2021:
            assert vals.index(one) < vals.index(two), f"Order {order} ({dates}): one-year and two-year rates out of column order: {vals}"
        if year in TWO_YEAR_SECOND_YEAR:
            assert TWO_YEAR_SECOND_YEAR[year] in vals, f"Order {order}: second-year rate {TWO_YEAR_SECOND_YEAR[year]}% not in its row: {vals}"
        rows.append({"year": year, "order": order, "dates": dates, "one_year": one, "two_year": two})
    years = [r["year"] for r in rows]
    assert years == sorted(years) and len(rows) == 58, "one row per guideline year: 1968 (a two-year order), then 1970 to 2026"
    one_freezes = [r["year"] for r in rows if r["one_year"] == 0]
    two_full = [r["year"] for r in rows if r["two_year"] == 0 and r["year"] not in TWO_YEAR_FIRST_YEAR_ONLY]
    assert one_freezes == [2015, 2016, 2020, 2026] and two_full == [2026], (one_freezes, two_full)
    since = {y: next(r for r in rows if r["year"] == y) for y in (2023, 2024, 2025)}
    factor = float(np.prod([1 + since[y]["one_year"] / 100 for y in (2023, 2024, 2025)]))
    pre = [r["one_year"] for r in rows if r["year"] < 2015]
    recent = [r for r in rows if r["year"] >= 2014]
    top = max(recent, key=lambda r: r["one_year"])
    return {"table": rows, "one_year_freezes": one_freezes, "two_year_full_freezes": two_full,
            "one_year_max_since_2014": top["one_year"], "one_year_max_since_2014_year": top["year"],
            "one_year_6plus_years_before_1990": sum(r["one_year"] >= 6 for r in rows if r["year"] < 1990),
            "orders_before_1990": sum(r["year"] < 1990 for r in rows),
            "rolled_forward_factor": K.r(factor, 4), "rolled_forward_pct": K.r(100 * (factor - 1), 1),
            "orders_since_survey": [since[y]["order"] for y in (2023, 2024, 2025)],
            "one_year_max": max(r["one_year"] for r in rows), "one_year_max_year": max(rows, key=lambda r: r["one_year"])["year"],
            "one_year_median_before_2015": K.r(float(np.median(pre)), 2),
            "first_year": rows[0]["year"], "orders": len(rows),
            "two_year_second_year": {str(y): v for y, v in TWO_YEAR_SECOND_YEAR.items()}}


# ====================================================================================
# The households: who lives in stabilized apartments (2023 NYCHVS)
# ====================================================================================

def lease_phase(h: pd.DataFrame) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Order #57's rate for each household, and the share of its renewal falling in the freeze year.

    LEASE_LENGTH (codebook p. 115): 1 under a year, 2 one year, 3 more than one but less than
    two, 4 two years, 5 more than two, -2 no active lease, -1 not reported. The orders deem a
    tenancy of up to one year a one-year lease and one of over a year up to and including two
    a two-year lease, and set no rate for a longer one. The main estimate gives the rest (no
    active lease, none reported, longer than two years: `other`, about 15%) the one-year rate
    and a renewal inside the window; `lease_cases` imputes them or leaves them out instead.
    LEASE_START (month the current lease began; code 1 = January-February 2021, k >= 2 =
    March 2021 + (k - 2) months) fixes a two-year lease's phase; earlier or unreported starts
    count one half."""
    ll = h["LEASE_LENGTH"].to_numpy()
    two = np.isin(ll, [3, 4])
    rate = np.where(two, ORDER57["two_year"], ORDER57["one_year"])
    ls = h["LEASE_START"].to_numpy()
    month_index = np.where(ls == 1, 1, 2 + (ls - 2))
    year = 2021 + (month_index // 12)
    month = (month_index % 12) + 1
    order_year = np.where(month >= 10, year, year - 1)
    known = ls >= 1
    phase = np.where(known, ((FREEZE_ORDER_YEAR - order_year) % 2 == 0).astype(float), 0.5)   # a two-year lease's chance of renewing in the window
    in_window = np.where(two, phase, 1.0)
    other = np.isin(ll, [-2, -1, 5])
    return rate, in_window, two, other, phase


def household_estimates(orders: dict) -> tuple[dict, pd.DataFrame]:
    h, rw = H.households()
    T = H.Survey(h["FW"], rw)
    a, arw = H.allunits()
    A = H.Survey(a["FW"], arw)
    grp = h["group"]
    rs, mkt, ph = (grp == "Rent stabilized").to_numpy(), (grp == "Market rental").to_numpy(), (grp == "Public housing").to_numpy()
    renter = h["renter"].to_numpy()
    allhh = np.ones(len(h), bool)
    inc = h["inc"]
    rent = h["RENT_AMOUNT"].where(h["RENT_AMOUNT"] >= 0).astype(float)
    grent = h["GRENT"].where(h["GRENT"] > 0).astype(float)
    assist = (h["RENTASSIST"] == 1).to_numpy()
    voucher = (h["RENTASSIST_VOUCHER"] == 1).to_numpy()
    cat = h["RENTBURDEN_CAT"]
    b = pd.Series(recode.bands(inc, [hi for *_, hi in INC_BANDS[:-1]], [k for k, *_ in INC_BANDS], na=""), index=inc.index)
    out: dict = {}

    def rti(r_monthly, mask):
        ratio = (12 * r_monthly / inc).where(inc > 0).where(r_monthly.notna())
        return T.median(mask, ratio.to_numpy(dtype=float))

    # --- anchors: the survey reproduces HPD's and the board's published figures
    occ_r = ((a["OCC"] == 1) & a["CSR"].isin(H.RENTER_CSR)).to_numpy()
    vac = (a["OCC"] == 2).to_numpy()
    rows = [
        ("Net rental vacancy rate, 2023 (HPD)", 1.41, "pct", A.share(occ_r | vac, vac)),
        ("Renter households (HPD)", 2_324_000, "count", T.total(renter)),
        ("Rent-stabilized households (HPD)", 960_700, "count", T.total(rs)),
        ("Market-rental households (HPD)", 1_119_000, "count", T.total(mkt)),
        ("Public-housing households (HPD)", 167_700, "count", T.total(ph)),
        ("Median contract rent, stabilized (HPD)", 1500, "usd", T.median(rs, rent)),
        ("Median contract rent, market (HPD)", 2000, "usd", T.median(mkt, rent)),
        ("Median gross rent, stabilized (RGB)", 1570, "usd", T.median(rs, grent)),
        ("Median 2022 income, stabilized (HPD, RGB)", 60_000, "usd", T.median(rs, inc)),
        ("Median 2022 income, market (HPD)", 90_800, "usd", T.median(mkt, inc)),
        ("Stabilized households with income of $100,000 or more (RGB)", 30, "pct", T.share(rs, (inc >= 100_000).to_numpy())),
        ("Market households with income of $100,000 or more (RGB)", 47, "pct", T.share(mkt, (inc >= 100_000).to_numpy())),
        ("Median gross rent-to-income, stabilized without a voucher (RGB)", 28.8, "pct", rti(grent, rs & ~voucher)),
    ]
    anchors = []
    for label, pub, kind, r in rows:
        est = 100 * r["est"] if kind == "pct" else r["est"]
        moe = 100 * r["moe90"] if kind == "pct" else r["moe90"]
        # Published shares are rounded to whole points or one decimal; counts to the hundred.
        tol = max(moe, 0.5 if kind == "pct" and float(pub).is_integer() else 0.05 if kind == "pct" else 50 if kind == "count" else 5)
        anchors.append({"statistic": label, "published": pub, "reproduced": K.r(est, 2 if kind == "pct" else 0), "moe90": K.r(moe, 2 if kind == "pct" else 0),
                        "n": r.get("n"), "within_margin": bool(abs(est - pub) <= tol), "kind": kind})
    nmt = rs & cat.isin([1, 2, 3]).to_numpy()
    out["anchors"] = anchors
    out["anchor_burden"] = {"burdened_30plus": pct_of(T.share(nmt, cat.isin([1, 2]).to_numpy())),
                            "severe_50plus": pct_of(T.share(nmt, (cat == 1).to_numpy()))}

    # --- who lives there
    prof = {}
    for g, m in (("stabilized", rs), ("market", mkt), ("public", ph), ("renters", renter)):
        d = {"households": count_of(T.total(m))}
        for k, *_ in INC_BANDS:
            d[f"inc_{k}"] = pct_of(T.share(m, (b == k).to_numpy()), 0)
        d["inc_100k_plus"] = pct_of(T.share(m, (inc >= 100_000).to_numpy()), 1)
        d["inc_under_50k"] = pct_of(T.share(m, (inc < 50_000).to_numpy()), 1)
        d["count_100k_plus"] = count_of(T.total(m & (inc >= 100_000).to_numpy()))
        d["count_200k_plus"] = count_of(T.total(m & (inc >= 200_000).to_numpy()))
        d["median_income"] = money_of(T.median(m, inc))
        d["median_contract_rent"] = money_of(T.median(m, rent))
        nm = m & cat.isin([1, 2, 3]).to_numpy()
        d["burden_30plus"] = pct_of(T.share(nm, cat.isin([1, 2]).to_numpy()), 0)
        d["burden_50plus"] = pct_of(T.share(nm, (cat == 1).to_numpy()), 0)
        d["householder_65plus"] = pct_of(T.share(m, (h["resp_age"] >= 65).to_numpy()), 0)
        d["lives_alone"] = pct_of(T.share(m, (h["HHSIZE"] == 1).to_numpy()), 0)
        yrs = (2023 - h["HHFIRSTMOVEIN"]).where(h["HHFIRSTMOVEIN"] > 0).astype(float)
        d["median_years_in_unit"] = {"years": K.r(T.median(m, yrs)["est"], 0), "n": int(m.sum())}
        d["in_unit_10plus"] = pct_of(T.share(m, (yrs >= 10).to_numpy()), 0)
        d["in_unit_20plus"] = pct_of(T.share(m, (yrs >= 20).to_numpy()), 0)
        d["rental_assistance"] = pct_of(T.share(m, assist), 0)
        prof[g] = d
    out["profile"] = prof
    hud = h["HUDILFY23"].astype(float)
    out["hud"] = {"low_income": pct_of(T.share(rs, ((hud >= 0) & (hud <= 80)).to_numpy() | (inc <= 0).to_numpy()), 0),
                  "above_120": pct_of(T.share(rs, (hud > 120).to_numpy()), 0),
                  "above_165": pct_of(T.share(rs, (hud > 165).to_numpy()), 0),
                  "low_income_among_100k_plus": pct_of(T.share(rs & (inc >= 100_000).to_numpy(), ((hud >= 0) & (hud <= 80)).to_numpy()), 0),
                  "under_100k_among_low_income": pct_of(T.share(rs & (((hud >= 0) & (hud <= 80)).to_numpy() | (inc <= 0).to_numpy()),
                                                                (inc < 100_000).to_numpy()), 0)}
    # The dollar line behind "low income": HUDILFY23 is income over HUD's FY2023 limit for the household's size, in whole
    # percent, so income / (HUDILFY23 / 100) recovers the limit; the median over households of one size removes the rounding.
    lim = {}
    for size in (1, 2, 3, 4):
        m = ((h["HHSIZE"] == size) & (hud >= 20) & (hud <= 240) & (inc > 0)).to_numpy()
        full = float(np.median((inc / (hud / 100)).to_numpy()[m]))
        lim[str(size)] = {"limit_100": K.r(full, -2), "low_income_80": K.r(0.8 * full, -2), "n": int(m.sum())}
    out["hud"]["implied_limits"] = lim
    out["hud"]["no_or_negative_income"] = pct_of(T.share(rs, (inc <= 0).to_numpy()), 0)
    out["hud"]["low_income_with_positive_income"] = pct_of(T.share(rs, ((hud >= 0) & (hud <= 80) & (inc > 0)).to_numpy()), 0)
    # HUD's ratio applies only to households with positive income; the share among those alone:
    out["hud"]["low_income_among_positive_income"] = pct_of(T.share(rs & (inc > 0).to_numpy(), ((hud >= 0) & (hud <= 80)).to_numpy()), 0)
    out["rent_by_band"] = {k: money_of(T.median(rs & (b == k).to_numpy(), rent)) for k, *_ in INC_BANDS}
    out["burden_30plus_by_band"] = {k: pct_of(T.share(nmt & (b == k).to_numpy(), cat.isin([1, 2]).to_numpy()), 0) for k, *_ in INC_BANDS}
    # The survey counts rent paid on no income as severe burden; the bottom band's share without those households:
    out["burden_30plus_lt25k_positive_income"] = pct_of(T.share(nmt & (b == "lt25k").to_numpy() & (inc > 0).to_numpy(), cat.isin([1, 2]).to_numpy()), 0)
    out["diffs_stabilized_minus_market"] = {
        k: {key: K.r(100 * v, 1) if key != "excludes_zero" else v for key, v in T.share_diff(rs, x, mkt, x).items()}
        for k, x in (("inc_100k_plus", (inc >= 100_000).to_numpy()), ("householder_65plus", (h["resp_age"] >= 65).to_numpy()),
                     ("lives_alone", (h["HHSIZE"] == 1).to_numpy()), ("in_unit_10plus", ((2023 - h["HHFIRSTMOVEIN"]) >= 10).to_numpy() & (h["HHFIRSTMOVEIN"] > 0).to_numpy()))}
    nm_mkt = mkt & cat.isin([1, 2, 3]).to_numpy()
    for k, x in (("burden_30plus", cat.isin([1, 2]).to_numpy()), ("burden_50plus", (cat == 1).to_numpy())):
        out["diffs_stabilized_minus_market"][k] = {key: K.r(100 * v, 1) if key != "excludes_zero" else v
                                                   for key, v in T.share_diff(nmt, x, nm_mkt, x).items()}
    out["borough_100k_plus"] = {bn: pct_of(T.share(rs & (h["boro_name"] == bn).to_numpy(), (inc >= 100_000).to_numpy()), 0)
                                for bn in H.BORO.values()}
    hi = (inc >= 100_000).to_numpy()
    out["borough_diffs_100k_plus"] = {
        f"{x}_minus_{y}": {key: K.r(100 * v, 1) if key != "excludes_zero" else v
                           for key, v in T.share_diff(rs & (h["boro_name"] == x).to_numpy(), hi, rs & (h["boro_name"] == y).to_numpy(), hi).items()}
        for x, y in (("Manhattan", "Queens"), ("Manhattan", "Brooklyn"), ("Queens", "Brooklyn"), ("Brooklyn", "Bronx"), ("Queens", "Bronx"), ("Manhattan", "Bronx"))}

    # --- the freeze's value (computed)
    rate57, in_window, two, other, phase = lease_phase(h)
    rs_rent = rs & rent.notna().to_numpy()
    save = 12 * rent.fillna(0).to_numpy() * rate57          # annual saving at renewal, 2023 rents
    f = orders["rolled_forward_factor"]
    q = T.quantiles(allhh, inc, [0.2, 0.4, 0.6, 0.8])
    qi = np.searchsorted(np.array(q), inc.to_numpy(), side="right")
    tot = T.total(rs_rent, save)
    first = T.total(rs_rent, save * in_window)
    med = T.median(rs_rent, save)
    fz = {"counterfactual": "Order #57 repeated: 3% for a one-year lease, 4.5% for a two-year lease",
          "households_with_rent": count_of(T.total(rs_rent)),
          "lease_two_year": pct_of(T.share(rs, two), 0), "renews_in_freeze_year": pct_of(T.ratio(rs, in_window, np.ones(len(h))), 0),
          "median_saving": money_of(med), "mean_saving": money_of(T.mean(rs_rent, save)),
          "median_saving_rolled": {"usd": K.r(med["est"] * f, 0)},
          "total_all_renewed": {"usd": K.r(tot["est"], -6), "moe": K.r(tot["moe90"], -6)},
          "total_all_renewed_rolled": {"usd": K.r(tot["est"] * f, -6)},
          "total_freeze_year": {"usd": K.r(first["est"], -6), "moe": K.r(first["moe90"], -6)},
          "total_freeze_year_rolled": {"usd": K.r(first["est"] * f, -6)},
          "quintile_cuts": [K.r(x, 0) for x in q]}
    shares = {}
    for k, lab, lo, hi in INC_BANDS:
        m = (b == k).to_numpy()
        shares[k] = {"households": pct_of(T.share(rs_rent, m), 1), "dollars": pct_of(T.ratio(rs_rent, save * m, save), 1),
                     "median_saving": money_of(T.median(rs_rent & m, save)),
                     "saving_pct_income": pct_of(T.median(rs_rent & m & (inc > 0).to_numpy(), save / inc.where(inc > 0).to_numpy()), 1)}
    fifths = {}
    for k in range(5):
        m = qi == k
        fifths[f"q{k + 1}"] = {"households": pct_of(T.share(rs_rent, m), 1), "dollars": pct_of(T.ratio(rs_rent, save * m, save), 1)}
    fz["by_band"], fz["by_fifth"] = shares, fifths
    for key, m in (("100k_plus", (inc >= 100_000).to_numpy()), ("200k_plus", (inc >= 200_000).to_numpy()), ("under_50k", (inc < 50_000).to_numpy()),
                   ("under_100k", (inc < 100_000).to_numpy())):
        fz[f"households_{key}"] = pct_of(T.share(rs_rent, m), 1)
        fz[f"dollars_{key}"] = pct_of(T.ratio(rs_rent, save * m, save), 1)
    na = rs_rent & ~assist
    fz["no_assistance"] = {"dollars_100k_plus": pct_of(T.ratio(na, save * (inc >= 100_000).to_numpy(), save), 0),
                           "households_100k_plus": pct_of(T.share(na, (inc >= 100_000).to_numpy()), 0),
                           "pct_income_by_band": {k: pct_of(T.median(na & (b == k).to_numpy() & (inc > 0).to_numpy(), save / inc.where(inc > 0).to_numpy()), 1)
                                                  for k, *_ in INC_BANDS}}
    fz["assisted_by_band"] = {k: pct_of(T.share(rs_rent & (b == k).to_numpy(), assist), 0) for k, *_ in INC_BANDS}
    # The order's own reach: a one-year renewal in its window saves a year, a two-year renewal two years; a two-year
    # lease that began in the year before the freeze renews after the window closes and gets nothing from this order.
    fz["not_renewing_in_freeze_year"] = pct_of(T.ratio(rs, 1 - in_window, np.ones(len(h))), 0)
    span = save * in_window * np.where(two, 2, 1)
    order_tot = T.total(rs_rent, span)
    fz["total_order_leases"] = {"usd": K.r(order_tot["est"], -6), "moe": K.r(order_tot["moe90"], -6)}
    fz["total_order_leases_rolled"] = {"usd": K.r(order_tot["est"] * f, -6)}
    fz["order_weighted_dollars_100k_plus"] = pct_of(T.ratio(rs_rent, span * (inc >= 100_000).to_numpy(), span), 0)
    r0 = rent.fillna(0).to_numpy()
    fz["uniform_rate_dollars_100k_plus"] = pct_of(T.ratio(rs_rent, r0 * (inc >= 100_000).to_numpy(), r0), 0)
    fz["lease_two_year_100k_plus"] = pct_of(T.share(rs & (inc >= 100_000).to_numpy(), two), 0)
    fz["lease_two_year_under_100k"] = pct_of(T.share(rs & (inc < 100_000).to_numpy(), two), 0)
    fz["median_saving_monthly"] = {"usd": K.r(med["est"] / 12, 0)}
    fz["median_saving_monthly_rolled"] = {"usd": K.r(med["est"] * f / 12, 0)}
    mr = prof["stabilized"]["median_contract_rent"]["usd"]
    fz["example_at_median_rent"] = {"rent": mr, "one_year_monthly": K.r(mr * ORDER57["one_year"], 2), "two_year_monthly": K.r(mr * ORDER57["two_year"], 2)}
    # The households this order reaches: each weighted by the share of its renewal that falls in the window.
    W = H.Survey(h["FW"] * in_window, rw * in_window[:, None])
    reach = rs_rent & (in_window > 0)
    mw = W.median(reach, save)
    six = (inc >= 100_000).to_numpy()               # not `hi`: the band loop above rebinds it to a band edge
    fz["in_window"] = {"median_saving": money_of(mw), "median_saving_rolled": {"usd": K.r(mw["est"] * f, 0)},
                       "median_saving_monthly": {"usd": K.r(mw["est"] / 12, 0)},
                       "dollars_100k_plus": pct_of(T.ratio(rs_rent, save * in_window * six, save * in_window), 0),
                       "households_100k_plus": pct_of(T.ratio(rs, in_window * six, in_window), 0)}
    fz["assisted_share_of_dollars"] = pct_of(T.ratio(rs_rent, save * assist, save), 0)
    fz["no_or_negative_income_lt25k"] = pct_of(T.share(rs_rent & (b == "lt25k").to_numpy(), (inc <= 0).to_numpy()), 0)
    fz["no_lease_or_unreported"] = pct_of(T.share(rs, np.isin(h["LEASE_LENGTH"].to_numpy(), [-2, -1])), 0)
    fz["unknown_lease"] = pct_of(T.share(rs, other), 0)
    fz["unknown_lease_100k_plus"] = pct_of(T.share(rs & other, six), 0)
    fz["known_lease_100k_plus"] = pct_of(T.share(rs & ~other, six), 0)
    # Who gets this order: the distribution among the households modelled to renew in its window, each weighted
    # by the share of its renewal that falls there, and a year of each one's saving (outside review, 2026-10-02:
    # the title asks about this order, not about every household's next renewal; the all-household figures above
    # answer that wider question and stay in the FAQ).
    wsave = save * in_window
    ren = {"by_band": {}, "by_fifth": {}}
    for k, *_ in INC_BANDS:
        m = (b == k).to_numpy()
        ren["by_band"][k] = {"households": pct_of(T.ratio(rs_rent, in_window * m, in_window), 1),
                             "dollars": pct_of(T.ratio(rs_rent, wsave * m, wsave), 1),
                             "median_saving": money_of(W.median(reach & m, save)),
                             "saving_pct_income": pct_of(W.median(reach & m & (inc > 0).to_numpy(), save / inc.where(inc > 0).to_numpy()), 1)}
    for k in range(5):
        m = qi == k
        ren["by_fifth"][f"q{k + 1}"] = {"households": pct_of(T.ratio(rs_rent, in_window * m, in_window), 1),
                                        "dollars": pct_of(T.ratio(rs_rent, wsave * m, wsave), 1)}
    for key, m in (("100k_plus", six), ("200k_plus", (inc >= 200_000).to_numpy()), ("under_50k", (inc < 50_000).to_numpy()),
                   ("under_100k", (inc < 100_000).to_numpy())):
        ren[f"households_{key}"] = pct_of(T.ratio(rs_rent, in_window * m, in_window), 1)
        ren[f"dollars_{key}"] = pct_of(T.ratio(rs_rent, wsave * m, wsave), 1)
    lt25 = (b == "lt25k").to_numpy()
    ren["no_assistance_dollars_100k_plus"] = pct_of(T.ratio(na, wsave * six, wsave), 0)
    ren["no_assistance_households_100k_plus"] = pct_of(T.ratio(na, in_window * six, in_window), 0)
    ren["no_assistance_pct_income_lt25k"] = pct_of(W.median(na & (in_window > 0) & lt25 & (inc > 0).to_numpy(), save / inc.where(inc > 0).to_numpy()), 1)
    ren["uniform_rate_dollars_100k_plus"] = pct_of(T.ratio(rs_rent, r0 * in_window * six, r0 * in_window), 0)
    ren["assisted_lt25k"] = pct_of(T.ratio(rs_rent & lt25, in_window * assist, in_window), 0)
    ren["no_or_negative_income_lt25k"] = pct_of(T.ratio(rs_rent & lt25, in_window * (inc <= 0).to_numpy(), in_window), 0)
    ren["assisted_share_of_dollars"] = pct_of(T.ratio(rs_rent, wsave * assist, wsave), 0)
    ren["rental_assistance"] = pct_of(T.ratio(rs_rent, in_window * assist, in_window), 0)
    # A group's share of the savings minus its share of the households, computed on every replicate: both shares
    # come from the same sample, so the margin of the gap is not the two margins combined (outside review, 2026-10-02).
    wts = np.column_stack([T.w, T.rw])
    sv_w = (wsave * rs_rent)[:, None] * wts
    hh_w = (in_window * rs_rent)[:, None] * wts

    def gap(m: np.ndarray) -> dict:
        d = sv_w[m].sum(0) / sv_w.sum(0) - hh_w[m].sum(0) / hh_w.sum(0)
        moe = T.moe(float(d[0]), d[1:])
        return {"pts": K.r(100 * float(d[0]), 1), "moe": K.r(100 * moe, 1), "excludes_zero": bool(abs(float(d[0])) > moe)}

    ren["gaps"] = {key: gap(m) for key, m in (("100k_plus", six), ("200k_plus", (inc >= 200_000).to_numpy()),
                                              ("under_50k", (inc < 50_000).to_numpy()), ("under_100k", (inc < 100_000).to_numpy()))}
    for k, *_ in INC_BANDS:
        ren["by_band"][k]["gap"] = gap((b == k).to_numpy())
    for k in range(5):
        ren["by_fifth"][f"q{k + 1}"]["gap"] = gap(qi == k)
    fz["renewing"] = ren
    # The ~15% without a standard lease, three ways. A, the main estimate: the one-year rate and a renewal
    # in the window. B: lease length imputed from the one-/two-year mix of known leases in the household's
    # income band (expected values). C: left out.
    r_one, r_two, r0 = ORDER57["one_year"], ORDER57["two_year"], rent.fillna(0).to_numpy()
    p_two = np.zeros(len(h))
    for k, *_ in INC_BANDS:
        m = (b == k).to_numpy()
        p_two[m] = T.share(rs & m & ~other, two)["est"]
    sv_a = 12 * r0 * rate57
    sv_b = np.where(other, 12 * r0 * ((1 - p_two) * r_one + p_two * r_two), sv_a)
    win_b = np.where(other, 1 - 0.5 * p_two, in_window)
    year_a, span_a = sv_a * in_window, sv_a * in_window * np.where(two, 2, 1)
    year_b = np.where(other, 12 * r0 * ((1 - p_two) * r_one + 0.5 * p_two * r_two), year_a)
    span_b = np.where(other, sv_b, span_a)       # an imputed two-year renewal in the window keeps the freeze both years
    cases = {"A": (sv_a, in_window, year_a, span_a, np.ones(len(h), bool)),
             "B": (sv_b, win_b, year_b, span_b, np.ones(len(h), bool)),
             "C": (sv_a, in_window, year_a, span_a, ~other)}
    fz["lease_cases"] = {}
    for name, (sv, win, yr, sp, keep) in cases.items():
        base = rs_rent & keep
        Wc = H.Survey(h["FW"] * win, rw * win[:, None])
        fz["lease_cases"][name] = {"renews_in_freeze_year": pct_of(T.ratio(rs & keep, win, np.ones(len(h))), 0),
                                   "dollars_100k_plus": pct_of(T.ratio(base, sv * six, sv), 1),
                                   "renewing_dollars_100k_plus": pct_of(T.ratio(base, yr * six, yr), 1),
                                   "renewing_households_100k_plus": pct_of(T.ratio(base, win * six, win), 1),
                                   "median_saving_renewing": money_of(Wc.median(base & (win > 0), sv)),
                                   "year_renewing": {"usd": K.r(T.total(base, yr)["est"], -6)},
                                   "total_order_leases": {"usd": K.r(T.total(base, sp)["est"], -6)}}
    # Other benchmarks: the board's own preliminary ranges of 7 May 2026 (0-2% one-year, 0-4% two-year).
    fz["alternatives"] = {}
    for key, (a1, a2) in {"prelim_top": (0.02, 0.04), "prelim_mid": (0.01, 0.02)}.items():
        sv = 12 * r0 * np.where(two, a2, a1)
        t_alt = T.total(rs_rent, sv * in_window * np.where(two, 2, 1))
        fz["alternatives"][key] = {"one_year": K.r(100 * a1, 1), "two_year": K.r(100 * a2, 1),
                                   "total_order_leases": {"usd": K.r(t_alt["est"], -6), "moe": K.r(t_alt["moe90"], -6)},
                                   "year_renewing": {"usd": K.r(T.total(rs_rent, sv * in_window)["est"], -6)},
                                   "median_saving_renewing": money_of(W.median(reach, sv)),
                                   "dollars_100k_plus": pct_of(T.ratio(rs_rent, sv * six, sv), 0),
                                   "renewing_dollars_100k_plus": pct_of(T.ratio(rs_rent, sv * in_window * six, sv * in_window), 0)}
    g0 = grent.fillna(0).to_numpy()
    fz["gross_rent_base_dollars_100k_plus"] = pct_of(T.ratio(rs & grent.notna().to_numpy(), 12 * g0 * rate57 * (inc >= 100_000).to_numpy(), 12 * g0 * rate57), 0)
    fz["sweep_total_all_renewed"] = {f"{100 * r:g}": K.r(T.total(rs_rent, 12 * rent.fillna(0).to_numpy() * r)["est"], -6) for r in SWEEP}
    out["freeze"] = fz

    # --- the companion data: renter households, the variables used
    data = pd.DataFrame({
        "control": h["CONTROL"], "weight_FW": h["FW"], "group": grp, "borough": h["boro_name"], "income_2022": inc,
        "income_band": b.map({k: lab for k, lab, *_ in INC_BANDS}), "income_pct_of_hud_limit_fy2023": hud,
        "contract_rent_2023": rent, "gross_rent_2023": grent, "lease_length_code": h["LEASE_LENGTH"], "lease_start_code": h["LEASE_START"],
        "rental_assistance": assist.astype(int), "section8_voucher": voucher.astype(int), "rent_burden_category": cat,
        "householder_age": h["resp_age"], "household_size": h["HHSIZE"], "first_move_in_year": h["HHFIRSTMOVEIN"],
        "order57_rate": np.where(rs, rate57, np.nan), "annual_saving_at_2023_rent": np.where(rs_rent, np.round(save, 2), np.nan),
        "renews_in_freeze_year": np.where(rs, in_window, np.nan),
    })[renter].reset_index(drop=True)
    return out, data


# ====================================================================================
# The long view: every survey from 1993 to 2023, in 2022 dollars
# ====================================================================================

LONG_GROUPS = {"stabilized": "Rent stabilized", "market": "Private unregulated"}


def long_view() -> tuple[dict, pd.DataFrame]:
    """Thirty years of the survey (card nychvs_income): the share of stabilized and market-rate
    renter households with incomes of $100,000 or more in 2022 dollars, their medians, and
    six-figure households' share of the stabilized rent, which is how a freeze of one rate on
    every lease would split its dollars. Each survey is a separate snapshot; 1993-2014, 2017 and
    2021-2023 are separate stretches (the 2017 recode of stabilization, the 2021 redesign), and
    a change is named only within one."""
    years = [y for y in HI.YEARS if HI.available(y)]
    assert years[0] == 1993 and years[-1] == 2023 and len(years) == 11, years
    checks = [c for fn in (HI.check_published_tables, HI.check_published_reports, HI.check_gvf_against_replicates,
                           HI.check_unreported_1991) for c in fn()]
    assert all(c.passed for c in checks), [c.detail for c in checks if not c.passed]
    six, med, tilt, rows = {}, {}, {}, []
    for y in years:
        t = HI.rent_tilt(y)
        tilt[str(y)] = {"households_pct": K.r(t["households_pct"]), "households_moe": K.r(t["households_moe"]), "rent_pct": K.r(t["rent_pct"]),
                        "rent_moe": K.r(t["rent_moe"]), "rent_share": K.r(t["rent_pct"] / 100, 3),   # as a fraction, for the findings' words
                        "tilt": K.r(t["tilt"]), "tilt_moe": K.r(t["tilt_moe"]), "n": t["n"]}
        for k, g in LONG_GROUPS.items():
            s, m = HI.six_figure(y, g), HI.median_income(y, g)
            six.setdefault(k, {})[str(y)] = {"pct": K.r(s["pct"]), "moe": K.r(s["moe"]), "margin": s["margin"], "n": s["n"],
                                             "households": K.r(s["households"], -2)}
            med.setdefault(k, {})[str(y)] = {"usd": K.r(m["median"], 0), "usd_2022": K.r(m["median_2022"], -2),
                                             "moe_2022": K.r(m["moe_2022"], -2)}
            rows.append({"survey_year": y, "income_year": HI.income_year(y), "stretch": HI.stretch(y), "group": g,
                         "records": s["n"], "households": round(s["households"]), "six_figure_pct": K.r(s["pct"], 2),
                         "six_figure_moe": K.r(s["moe"], 2), "margin_method": s["margin"], "threshold_nominal": round(s["threshold"]),
                         "median_income_nominal": K.r(m["median"], 0), "median_income_2022": K.r(m["median_2022"], 0),
                         "six_figure_share_of_rent_pct": K.r(t["rent_pct"], 2) if k == "stabilized" else None})

    def change(k: str, a: int, b: int) -> dict:
        """b minus a for one group: two independent surveys within one stretch."""
        assert HI.stretch(a) == HI.stretch(b), (a, b)
        x, z = six[k][str(a)], six[k][str(b)]
        d, moe = z["pct"] - x["pct"], float(np.hypot(x["moe"], z["moe"]))
        return {"from": a, "to": b, "pts": K.r(d), "moe": K.r(moe), "excludes_zero": bool(abs(d) > moe)}

    def gap(y: int) -> dict:
        """Market minus stabilized in one survey (two disjoint groups; margins combined as independent)."""
        a, b = six["stabilized"][str(y)], six["market"][str(y)]
        return {"pts": K.r(b["pct"] - a["pct"]), "moe": K.r(float(np.hypot(a["moe"], b["moe"])))}

    def gap_change(a: int, b: int) -> dict:
        assert HI.stretch(a) == HI.stretch(b), (a, b)
        ga, gb = gap(a), gap(b)
        d, moe = gb["pts"] - ga["pts"], float(np.hypot(ga["moe"], gb["moe"]))
        return {"from": a, "to": b, "pts": K.r(d), "moe": K.r(moe), "excludes_zero": bool(abs(d) > moe)}

    def med_change(k: str, a: int, b: int) -> dict:
        """Median income in 2022 dollars, b minus a, where both files carry replicate weights."""
        assert HI.stretch(a) == HI.stretch(b), (a, b)
        x, z = med[k][str(a)], med[k][str(b)]
        d, moe = z["usd_2022"] - x["usd_2022"], float(np.hypot(x["moe_2022"], z["moe_2022"]))
        return {"from": a, "to": b, "usd": d, "moe": K.r(moe, -2), "excludes_zero": bool(abs(d) > moe)}

    tilts = [v["tilt"] for v in tilt.values()]
    cpi = BLS.cpi().set_index("year")["index_value"]
    from nycdata import hvs_history as HH          # the file names and the positions the two readers share
    layout = {str(y): {"file": HH.LAYOUT[y][0], "control_status": list(HH.LAYOUT[y][2]), "weight": HH.LAYOUT[y][3], "tenure": HH.LAYOUT[y][4],
                       "first_replicate_weight": HH.LAYOUT[y][5], "income": list(HI.INCOME[y]), "contract_rent": [HI.RENT[y], 5]}
              for y in years if y in HH.LAYOUT}
    out = {
        "years": years, "dollars": HI.DOLLARS,
        "stretches": {s: [y for y in years if HI.stretch(y) == s] for s in dict.fromkeys(HI.stretch(y) for y in years)},
        "stretch": {str(y): HI.stretch(y) for y in years},
        "income_year": {str(y): HI.income_year(y) for y in years},
        "threshold_nominal": {str(y): round(HI.threshold(y)) for y in years},
        "cpi_ny_all_items": {str(HI.income_year(y)): float(cpi[HI.income_year(y)]) for y in years},
        "six_figure": six, "median_income": med, "rent_tilt": tilt,
        "changes": {"stabilized_1993_2014": change("stabilized", 1993, 2014), "market_1993_2014": change("market", 1993, 2014),
                    "stabilized_2021_2023": change("stabilized", 2021, 2023), "market_2021_2023": change("market", 2021, 2023)},
        "gaps": {str(y): gap(y) for y in years},
        "gap_changes": {"1993_2014": gap_change(1993, 2014), "2021_2023": gap_change(2021, 2023)},
        "median_changes": {"stabilized_2021_2023": med_change("stabilized", 2021, 2023), "market_2021_2023": med_change("market", 2021, 2023)},
        "stabilized_1993_2014_range": [min(six["stabilized"][str(y)]["pct"] for y in years if HI.stretch(y) == "1993-2014"),
                                       max(six["stabilized"][str(y)]["pct"] for y in years if HI.stretch(y) == "1993-2014")],
        "tilt_range": [min(tilts), max(tilts)], "tilt_all_positive": bool(min(tilts) > 0),
        "market_households_growth_1993_2014_pct": K.r(100 * (six["market"]["2014"]["households"] / six["market"]["1993"]["households"] - 1)),
        "tilt_max_replicate_moe": max(v["tilt_moe"] for v in tilt.values() if v["tilt_moe"] is not None),
        "formula_moe_range": [min(v["moe"] for d in six.values() for v in d.values() if v["margin"] == "gvf"),
                              max(v["moe"] for d in six.values() for v in d.values() if v["margin"] == "gvf")],
        "unreported_1991_pct": K.r(next(c.value for c in checks if c.name == "unreported_1991")),
        "layout": layout,
        "checks": [{"name": c.name, "passed": c.passed, "detail": c.detail, "source": c.source} for c in checks],
    }
    return out, pd.DataFrame(rows)


# ====================================================================================
# Bar-chart cubes (the Chart component) and the manifest
# ====================================================================================

def chart_cubes(hv: dict) -> dict:
    P_ = hv["profile"]
    bands_series = [[k, lab] for k, lab, *_ in INC_BANDS]
    charts = {}
    charts["income_mix"] = {
        "title": "Stabilized tenants run from the poorest to six-figure earners",
        "subtitle": "Households by 2022 income, rent-stabilized and market-rate rentals, 2023",
        "note": "2023 New York City Housing and Vacancy Survey, households by their 2022 income. Ticks are 90% margins of error from the survey's replicate weights.",
        "controls": [], "metrics": [{"key": "pair", "label": "Share of households", "format": "pct", "axis": "Share of households (%)"}],
        "series": bands_series, "paired": [["rs", "Rent stabilized"], ["mkt", "Market rental"]],
        "cells": {"": {k: {"n": P_["stabilized"][f"inc_{k}"]["n"], "rs": P_["stabilized"][f"inc_{k}"]["pct"], "rs_moe": P_["stabilized"][f"inc_{k}"]["moe"],
                           "mkt": P_["market"][f"inc_{k}"]["pct"], "mkt_moe": P_["market"][f"inc_{k}"]["moe"]} for k, *_ in INC_BANDS}},
        "provenance": "2023 NYCHVS public use files"}
    sb = hv["freeze"]["renewing"]["by_band"]
    charts["saving"] = {
        "title": "Relative to income, the saving shrinks as income rises",
        "subtitle": "The freeze's yearly saving for rent-stabilized households renewing in its window: 2023 rents against 2022 incomes",
        "note": "Computed: 12 months of the 2023 contract rent times Order #57's rate (3% one-year, 4.5% two-year), the benchmark the freeze is set against, "
                "for the households modelled to renew between October 2026 and September 2027 (by their 2023 lease dates). "
                "The share of income is each household's saving over its 2022 income, the median within each band; households with no or negative income are left out of it. "
                "Ticks are 90% sampling margins; they leave out the uncertainty of the model itself (who renews, and what rent and lease a household has in 2026).",
        "controls": [], "metrics": [{"key": "pct_income", "label": "Saving as a share of income", "format": "pct", "axis": "Median saving as a share of 2022 income (%)"},
                                    {"key": "usd", "label": "Saving in dollars", "format": "dollar", "axis": "Median saving a year at 2023 rents ($)"}],
        "series": bands_series,
        "cells": {"": {k: {"n": sb[k]["median_saving"]["n"], "pct_income": sb[k]["saving_pct_income"]["pct"], "pct_income_moe": sb[k]["saving_pct_income"]["moe"],
                           "usd": sb[k]["median_saving"]["usd"], "usd_moe": sb[k]["median_saving"]["moe"]} for k, *_ in INC_BANDS}},
        "provenance": "2023 NYCHVS public use files; RGB Orders #57 and #58"}
    bo = hv["borough_100k_plus"]
    charts["boroughs"] = {
        "title": "Six-figure stabilized tenants are rarest in the Bronx",
        "subtitle": "Rent-stabilized households with 2022 income of $100,000 or more, by borough, 2023",
        "note": "Staten Island has too few stabilized households in the survey to show. Manhattan and Queens differ by less than the margin of error. Ticks are 90% margins of error.",
        "controls": [], "metrics": [{"key": "share", "label": "Share of stabilized households", "format": "pct", "axis": "Income of $100,000 or more (%)"}],
        "series": [[k, k] for k in ("Bronx", "Brooklyn", "Queens", "Manhattan", "Staten Island")],
        "cells": {"": {k: ({"n": v["n"], "thin": True} if v["n"] < THIN else {"n": v["n"], "share": v["pct"], "share_moe": v["moe"]}) for k, v in bo.items()}},
        "provenance": "2023 NYCHVS public use files"}
    rows = [("householder_65plus", "Householder 65 or older"), ("lives_alone", "Lives alone"), ("in_unit_10plus", "In the apartment 10 years or more"),
            ("in_unit_20plus", "In the apartment 20 years or more"), ("burden_30plus", "Rent over 30% of income"), ("burden_50plus", "Rent over 50% of income")]
    charts["tenants"] = {
        "title": "Older, more often alone, and longer in place",
        "subtitle": "Rent-stabilized and market-rate renter households, 2023",
        "note": "Rent burden is contract rent against 2022 income among households outside means-tested housing (no Section 8 voucher or public housing), the board's convention. Ticks are 90% margins of error.",
        "controls": [], "metrics": [{"key": "pair", "label": "Share of households", "format": "pct", "axis": "Share of households (%)"}],
        "series": [[k, lab] for k, lab in rows], "paired": [["rs", "Rent stabilized"], ["mkt", "Market rental"]],
        "cells": {"": {k: {"n": P_["stabilized"][k]["n"], "rs": P_["stabilized"][k]["pct"], "rs_moe": P_["stabilized"][k]["moe"],
                           "mkt": P_["market"][k]["pct"], "mkt_moe": P_["market"][k]["moe"]} for k, _ in rows}},
        "provenance": "2023 NYCHVS public use files"}
    return {"charts": charts, "min_records": 1, "source": "2023 New York City Housing and Vacancy Survey public use files; NYC Rent Guidelines Board orders."}


def main(argv=None) -> int:
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.parse_args(argv)
    res["published"] = quoted_figures()
    res["orders"] = orders_history()
    hv, data = household_estimates(res["orders"])
    res["hvs"] = hv
    res["reproduction"] = [{k: a[k] for k in ("statistic", "published", "reproduced", "within_margin", "kind")} for a in hv["anchors"]]
    res["reproduction_source"] = "HPD's 2023 Selected Initial Findings and the Rent Guidelines Board's 2024 Income and Affordability Study"
    res["reproduction_note"] = ("Tolerance: the survey estimate's 90% margin, or the published figure's own rounding where that is wider "
                                "(half a point for whole-number shares). Each published figure is checked against the saved report before it is compared.")
    res["as_of"] = "2023 survey (incomes 2022); orders to #58 (leases from 1 October 2026); long view, surveys of 1993 to 2023"
    assert all(a["within_margin"] for a in hv["anchors"]), [a for a in hv["anchors"] if not a["within_margin"]]
    res["long_view"], long_rows = long_view()
    charts = chart_cubes(hv)
    PATHS.json("results.json", res, nd=K.size_decimals)
    PATHS.json("charts.json", charts, compact=True, nd=K.size_decimals)
    from article_kit import fig_from_cube, grouped_bars  # static twins of the bar charts, for readers without scripts
    plain = json.loads(json.dumps(K.clean(charts, nd=K.size_decimals)))["charts"]
    figures = ((2, "income_mix"), (5, "saving"), (6, "boroughs"), (7, "tenants"))     # Figures 1, 3 and 4 are template charts (viz.json)
    for stale in [*OUT.glob("fig*.svg"), *(OUT / "pdf").glob("fig*.pdf")]:
        if stale.stem not in {f"fig{n}_{cid}" for n, cid in figures}:
            stale.unlink()                                 # a figure renumbered: drop the old file (the build drops the served copy)
    for n, cid in figures:
        if cid == "saving":                               # shares under 10%: one decimal, as the interactive chart prints them
            mm, row = plain[cid]["metrics"][0], plain[cid]["cells"][""]
            rows = [(s[1], [(row[s[0]][mm["key"]], row[s[0]].get(mm["key"] + "_moe") or 0.0)]) for s in plain[cid]["series"]]
            grouped_bars(rows, OUT / f"fig{n}_{cid}.svg", plain[cid]["title"], plain[cid]["subtitle"], [mm["label"]], fmt="{:.1f}%",
                     height=1.6 + 0.3 * len(rows), xlabel=mm.get("axis"))
        else:
            fig_from_cube(plain[cid], OUT / f"fig{n}_{cid}.svg")
        PATHS.serve(OUT / f"fig{n}_{cid}.svg")
    PATHS.csv("data.csv", data)
    PATHS.csv("long_view.csv", long_rows)
    PATHS.serve(Path(__file__), "analysis.py")
    # The survey loaders, which analysis.py imports, published beside it so a reader can inspect every step
    # (outside review, 2026-10-02): the 2023 loader, the reader of the 1993-2021 files for the long view (its file
    # names and positions are also in results.json, long_view.layout), and the replicate-weight variance code both use.
    served = {"hvs.py": H.__file__, "hvs_income.py": HI.__file__, "replicates.py": replicates.__file__}
    for name, module in served.items():
        PATHS.serve(Path(module), name)                     # a loader no longer served is removed by the build
    lv = res["long_view"]
    print(f"long view: stabilized six-figure {lv['six_figure']['stabilized']['1993']['pct']}% (1993) -> "
          f"{lv['six_figure']['stabilized']['2017']['pct']}% (2017); market {lv['six_figure']['market']['1993']['pct']}% -> "
          f"{lv['six_figure']['market']['2017']['pct']}%; rent tilt {lv['tilt_range'][0]} to {lv['tilt_range'][1]} points")
    fz = hv["freeze"]
    print(f"median saving {fz['median_saving']['usd']} (rolled {fz['median_saving_rolled']['usd']}); $100k+ households "
          f"{fz['households_100k_plus']['pct']}% get {fz['dollars_100k_plus']['pct']}% of the dollars; total {fz['total_all_renewed']['usd']:,.0f}")
    print(f"outputs -> {PATHS.out} and {PATHS.site}")
    return 0


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