"""The six-figure renter — every number and figure in the article, with the working. py src/article_six_figure_renter.py # default seed, writes outputs py src/article_six_figure_renter.py --seed 7 # different bootstrap draws This file is the article's notebook. The article itself is short and keeps to what was found; the how and the why live here, section by section, next to the code that computes them. Read the docstrings top to bottom and you have the method. ====================================================================================== 1. The question, and the two datasets that answer it ====================================================================================== How many New York City renter households earning six figures pay 30% or more of their income in rent, who are they, and what does "burdened" mean at that income? "Six figures" means household income of $100,000 or more (2024 dollars). The first two versions of this article drew the line at $150,000; a reader pointed out that six figures means $100,000, and the analysis was redone at that line. The $150,000 and $200,000 cuts are still reported once, as context. Two datasets are used, and the article is explicit about which one each number comes from. * **The census microdata directly.** The 2020–2024 five-year American Community Survey Public Use Microdata Sample (PUMS) for the study area, with its 80 replicate weights. Every count and share in the article that is labelled "census" is a weighted tabulation of these records with a proper margin of error. Anyone with the PUMS files can reproduce them. * **The SparkyData synthetic population** (release 0.2, 500,000 people in 195,799 households, seed 20240101). Whole households are sampled from the same PUMS and re-weighted to published margins, so income, rent, size, borough, move-in year, building age and every person's earnings arrive *together* from a real disclosure-protected record. The article uses it to show that a synthetic population reproduces the census's own answer within its margin, and as the joined, queryable database the profile is built from. Nothing in the article depends on a modelled layer (see section 8 for why the wealth layer was dropped). ====================================================================================== 2. Families and singles first; roommates separately ====================================================================================== Household income is the sum of every member's income, related or not, so three graduates earning $60,000 each are a $180,000 "six-figure household". That is a different thing from a person or a family earning six figures, and the article keeps them apart: * **Families and singles** ("families" in the code): the census's family households (married couples and other families) plus people living alone. In seven of ten of these households above $100,000, one person alone earns $100,000 or more. This is the article's main subject. * **Roommate households**: the census's nonfamily households where the householder does not live alone (HHT codes 5 and 7). Unmarried partners fall in here too, because the census cannot tell a partner from a roommate at this level, so the category is an upper bound on roommates in the everyday sense. The article gives them their own section. ====================================================================================== 3. Direct census estimates: successive difference replication ====================================================================================== PUMS ships 80 replicate weights per household. For any statistic θ computed with the full weight, the same statistic is recomputed with each replicate weight r, and the variance is Var(θ) = (4 / 80) · Σ_r (θ_r − θ)² The Census Bureau's published margins are 90% intervals, so MOE = 1.645 · SE, and that is the convention this script uses on the census side. Counts (weighted sums) and shares (weighted means) get the same treatment. Medians are reported without a margin; they are weighted medians of the records in the group. ====================================================================================== 4. The synthetic population: what it adds, and how its uncertainty is measured ====================================================================================== The synthetic sample is unweighted, so its uncertainty comes from a household-level bootstrap: 2,000 resamples with replacement of the households in the group, the statistic recomputed each time, and the 5th and 95th percentiles reported as a **90% interval** — the same level as the census margins. The bootstrap is seeded. The synthetic population is a fifth the size of the PUMS pool for the study area. It does not add information about the joint distribution that the PUMS lacks — it is sampled from the PUMS. What it adds is a single coherent database in which households, the people in them, their earnings, their housing and (in other layers) their spending and wealth are already joined. For an article that uses only inherited variables, the honest description is: the census gives the count, the synthetic population reproduces it, and the profile could be computed from either. Counts from the synthetic sample are scaled to real households by the ratio of census-weighted NYC renter households with positive income to synthetic ones (about 41 real households per synthetic one). ====================================================================================== 5. The reproduction check: what the bands are and what "held out" means ====================================================================================== The renter cost-burden targets were computed from PUMS after the build and were never calibration inputs. Reproducing them is *not* an out-of-sample test — the synthetic households are sampled from the same PUMS, so the joint is inherited by construction. It is a check that calibration (re-weighting to region, tenure, size, income quintile, householder age and the person margins) did not distort the rent-by-income joint. That is worth checking and it passes. The band on each row comes from the release's verification framework: band = sqrt( MOE_pub² + B_samp² ), B_samp = 1.96 · sqrt( p(1−p) / (n / deff) ) where MOE_pub is the published 90% margin, n is the size of the statistic's own denominator in the synthetic sample, and deff ≈ 1.96 is the design effect the calibration reported (it arises from weight variation). The two terms are at different confidence levels (90% and 95%), which changes bands by about 2%. The median-rent row uses the framework's quantile band. ====================================================================================== 6. Definitions ====================================================================================== * **Gross rent** (GRNTP): contract rent plus utilities and fuel if paid separately. * **Rent burden** (GRPIP): gross rent as a percentage of household income, the census's own variable, top-coded at 101. Burdened = 30 or more; severely burdened = 50 or more. Households with zero or negative income are excluded. * **Move-in period** (MV): when the householder moved into the unit, seven bands from "12 months or less" to "30 years or more". * **People and bedrooms**: NP (persons in the household) and BDSP (bedrooms). * **The stabilised stock, by proxy.** PUMS does not record rent regulation. Rent stabilisation in New York covers buildings of six or more units built before 1974 (plus newer buildings with tax abatements). The closest PUMS cut is buildings with five or more units (BLD ≥ 6, whose categories are 5–9, 10–19, 20–49, 50+) built before 1970 (YRBLT ≤ 1960, whose categories run by decade). The article calls this "the stock rent stabilisation covers" and says it is a proxy. Half of all NYC renter households live in it. * **Owner cost burden** (OCPIP): selected monthly owner costs (mortgage payment including principal, property taxes, insurance, utilities, condo fees) as a percentage of household income, same 30% threshold, owners with a mortgage. * **Income bands** are in constant 2024 dollars (ADJINC applied). ====================================================================================== 7. The sanity check on long-tenure rents ====================================================================================== A six-figure household thirty years into a lease pays a median gross rent of about $1,500. A reader asked whether that is plausible and why the owner has not evicted them. It is plausible, and the owner cannot: * The 2023 NYC Housing and Vacancy Survey puts the citywide median rent at $1,641, the median rent-stabilised rent at $1,500 and rent-controlled at $988 (2023 dollars). A long-tenured household in a stabilised unit paying about $1,500 is the typical stabilised tenant, not an anomaly. * Rent-stabilised tenants have a statutory right to renew their lease. Before June 2019 a unit could be deregulated if the tenants' income exceeded $200,000 for two consecutive years *and* the rent was above a threshold; the Housing Stability and Tenant Protection Act of 2019 repealed that, so income is no longer a ground for deregulation or non-renewal at all. A $150,000 household was never over the old income line either. * The Furman Center reports an average tenure of eight years in stabilised units against three in unregulated ones, and that about 65% of the stabilised stock is pre-1974 buildings. * In the PUMS itself, the share of six-figure family and single households living in pre-1970 buildings of five or more units rises from about 40% among those who moved in within a year to about 80% among those in place thirty years or more, and the householder's median age rises from 33 to 68. Long tenure, old multi-unit buildings and low rents go together, which is what stabilisation predicts. The general tenure discount (landlords raising sitting tenants' rents more slowly than the market) is folded in and cannot be separated here. ====================================================================================== 8. Why the wealth layer is not in the article ====================================================================================== The first version reported liquid assets and net worth from the release's modelled balance-sheet layer. The layer assigns each household to a Survey of Consumer Finances donor cell by the household's income *rank within the metro*, mapped onto *national* income percentile bands, so a $190,000 New York household draws from the national p60–80 cell (median income about $116,000). Most of the group also drew from fallback cells that drop education or tenure, and rent is not a conditioning variable. The numbers are still computed below under `not_in_article` so the reasoning can be checked. ====================================================================================== 9. What changed across versions ====================================================================================== * v1 (2026-09-02): 100,000-person release; $150,000 line; wealth section. * v2 (2026-09-03): 500,000-person release; census side with replicate weights; retiree and single-earner findings withdrawn (the second came from PUMS WIF, workers in *family*, undefined for nonfamily households); wealth withdrawn. * v3 (2026-09-03): line moved to $100,000; families and singles separated from roommate households; move-in analysis extended to people, bedrooms, householder age and building stock; long-tenure rents sanity-checked against the HVS and the 2019 rent law. Outputs (all regenerable; nothing is edited by hand): output/articles/six_figure_renter/results.json every number quoted in the prose output/articles/six_figure_renter/pums_direct.csv the census-side estimates with margins output/articles/six_figure_renter/reproduction.csv the held-out reproduction table output/articles/six_figure_renter/fig*.svg figures output/articles/six_figure_renter/data.csv companion dataset (re-randomised ids) site/public/articles/six-figure-renter/ the same, copied for the website """ from __future__ import annotations import argparse import json import shutil import sqlite3 from pathlib import Path import numpy as np import pandas as pd from _paths import P, INTERIM, TARGETS, DB_PATH, OUTPUT import constraints_nyc as C import decode as D from sparkylib import db as sdb from sparkylib.series import SeriesRegistry BASE_YEAR = 2024 SLUG = "six-figure-renter" OUT = OUTPUT / "articles" / "six_figure_renter" SITE_OUT = P.root / "site" / "public" / "articles" / SLUG BANDS = D.RENT_BURDEN_INCOME_BANDS BAND_LABELS = { "lt35k": "under $35k", "35k_75k": "$35k–75k", "75k_100k": "$75k–100k", "100k_150k": "$100k–150k", "150k_200k": "$150k–200k", "200k_plus": "$200k+", } HIGH_INCOME = 100_000 # "six figures" BURDEN = 30.0 # gross rent (or owner costs) as % of household income SEVERE = 50.0 Z90 = 1.645 # ACS margins are 90% intervals; both sides use 90% CI = 0.90 REPS = 2000 MV_BANDS = [(1, "12 months or less"), (2, "13 to 23 months"), (3, "2 to 4 years"), (4, "5 to 9 years"), (5, "10 to 19 years"), (6, "20 to 29 years"), (7, "30 years or more")] MV_SHORT = {1: "<1 yr", 2: "1–2 yrs", 3: "2–4 yrs", 4: "5–9 yrs", 5: "10–19 yrs", 6: "20–29 yrs", 7: "30+ yrs"} OWNER_BANDS = [("100k_150k", "$100k–150k", 100_000, 150_000), ("150k_200k", "$150k–200k", 150_000, 200_000), ("200k_plus", "$200k+", 200_000, float("inf"))] OLD_YEARS = {"1939 or earlier", "1940 to 1949", "1950 to 1959", "1960 to 1969"} BIG_BUILDINGS = {"5-9 Apartments", "10-19 Apartments", "20-49 Apartments", "50 or more apartments"} GROUPS = (("families", "Families and singles"), ("roommates", "Roommate households")) INK = "#6f6e69" GRID = "#d5d4cd" BLUE = "#2a78d6" ORANGE = "#eb6834" FIG_W = 6.4 # ------------------------------------------------------------------ census side def load_pums() -> tuple[pd.DataFrame, np.ndarray]: """Study-area PUMS housing records with replicate weights, geography, and per-household person summaries (highest individual income, adults, earners).""" h = pd.read_csv(INTERIM / "metro_housing.csv.gz", dtype={"STATE": str, "PUMA": str}, usecols=["SERIALNO", "PUMA", "STATE", "WGTP", "NP", "TEN", "HINCP", "GRNTP", "GRPIP", "OCPIP", "MV", "HHT", "HHLDRAGEP", "R18", "BDSP", "BLD", "YRBLT"]) rw = np.load(INTERIM / "metro_housing_repwts.npy") assert len(rw) == len(h), "replicate weights must align row-for-row with the CSV" geo = pd.read_csv(P.reference / "geography_pumas.csv", dtype={"state_fips": str, "puma": str}) geo["puma_geoid"] = geo["state_fips"].str.zfill(2) + geo["puma"].str.zfill(5) h["puma_geoid"] = h["STATE"].str.zfill(2) + h["PUMA"].str.zfill(5) h = h.merge(geo[["puma_geoid", "in_nyc", "borough"]], on="puma_geoid", how="left") pp = pd.read_csv(INTERIM / "metro_person.csv.gz", usecols=["SERIALNO", "PINCP", "AGEP"]) pp["adult"] = (pp["AGEP"] >= 18).astype(int) pp["earner"] = (pp["PINCP"].fillna(0) >= 1000).astype(int) agg = pp.groupby("SERIALNO").agg(max_pinc=("PINCP", "max"), n_adults=("adult", "sum"), n_earners=("earner", "sum")) h = h.merge(agg, left_on="SERIALNO", right_index=True, how="left") h["roommates"] = h["HHT"].isin([5, 7]) h["alone"] = h["HHT"].isin([4, 6]) h["old_multi"] = (h["YRBLT"] <= 1960) & (h["BLD"] >= 6) # pre-1970, 5+ units h["residual"] = h["HINCP"] - 12 * h["GRNTP"].fillna(0) h["ratio"] = 12 * h["GRNTP"].fillna(0) / h["HINCP"].where(h["HINCP"] > 0) * 100 return h, rw class Census: """Weighted estimates with successive-difference-replication margins (90%).""" def __init__(self, h: pd.DataFrame, rw: np.ndarray): self.h, self.rw, self.w = h, rw, h["WGTP"].to_numpy(dtype=float) def share(self, mask, num) -> tuple[float, float]: mask = np.asarray(mask); num = np.asarray(num) w, rw = self.w[mask], self.rw[mask] x = num[mask].astype(float) full = (w * x).sum() / w.sum() reps = (rw * x[:, None]).sum(0) / rw.sum(0) return float(full), float(Z90 * np.sqrt(4 / 80 * ((reps - full) ** 2).sum())) def count(self, mask) -> tuple[float, float]: mask = np.asarray(mask) full = self.w[mask].sum() reps = self.rw[mask].sum(0) return float(full), float(Z90 * np.sqrt(4 / 80 * ((reps - full) ** 2).sum())) def median(self, mask, col: str) -> float: mask = np.asarray(mask) v = self.h.loc[mask, col].to_numpy(dtype=float) w = self.w[mask] keep = ~np.isnan(v) v, w = v[keep], w[keep] o = np.argsort(v) cw = np.cumsum(w[o]) return float(v[o][np.searchsorted(cw, cw[-1] / 2)]) def mean(self, mask, col: str) -> float: mask = np.asarray(mask) return float(np.average(self.h.loc[mask, col].to_numpy(dtype=float), weights=self.w[mask])) TRAITS = [ ("manhattan", "Lives in Manhattan"), ("recent_mover", "Moved in within 2 years"), ("long_tenure", "In the same home 10+ years"), ("lives_alone", "Lives alone"), ("married_couple", "Married couple"), ("has_children", "Children under 18"), ("has_100k_earner", "Someone in it earns $100k+"), ("two_plus_earners", "Two or more earners"), ("head_under_35", "Householder under 35"), ("head_65_plus", "Householder 65 or older"), ("rent_3000_plus", "Gross rent $3,000+"), ("old_multi", "Pre-1970 building, 5+ units"), ] def census_estimates(cs: Census) -> tuple[dict, pd.DataFrame]: h = cs.h rows: list[dict] = [] def put(key, label, value, moe, n, kind): rows.append({"key": key, "label": label, "value": value, "moe90": moe, "n_records": int(n), "kind": kind}) nyc = ((h["in_nyc"] == 1) & (h["TEN"] == 3) & (h["NP"] > 0) & (h["HINCP"] > 0) & h["GRPIP"].notna()).to_numpy() burdened = (h["GRPIP"] >= BURDEN).to_numpy() severe = (h["GRPIP"] >= SEVERE).to_numpy() inc = h["HINCP"].to_numpy(dtype=float) mv = h["MV"].to_numpy() rm = h["roommates"].to_numpy() hi = nyc & (inc >= HIGH_INCOME) groups = {"families": hi & ~rm, "roommates": hi & rm} res: dict = {} n, m = cs.count(nyc); res["nyc_renter_households"] = {"value": n, "moe90": m} put("nyc_renters", "NYC renter households with positive income", n, m, nyc.sum(), "count") n, m = cs.count(nyc & burdened); res["nyc_burdened_households"] = {"value": n, "moe90": m} put("nyc_burdened", "NYC rent-burdened renter households", n, m, (nyc & burdened).sum(), "count") # The count at three lines, all household types (context). res["lines"] = {} for line, key in ((100_000, "100k"), (150_000, "150k"), (200_000, "200k")): mk = nyc & (inc >= line) p, pm = cs.share(mk, burdened); c, cm = cs.count(mk); b, bm = cs.count(mk & burdened); s, sm = cs.share(mk, severe) res["lines"][key] = {"households": c, "households_moe90": cm, "burdened_share": p, "burdened_share_moe90": pm, "burdened_households": b, "burdened_households_moe90": bm, "severe_share": s, "severe_share_moe90": sm, "n_records": int(mk.sum()), "n_burdened_records": int((mk & burdened).sum())} put(f"share_{key}", f"Burdened share, NYC renters ${line:,}+", p, pm, mk.sum(), "share") put(f"count_{key}", f"Burdened NYC renter households ${line:,}+", b, bm, (mk & burdened).sum(), "count") # The gradient for figure 1. grad = [] for lo, hi_, lab in BANDS: mk = nyc & (inc >= lo) & (inc < hi_) p, pm = cs.share(mk, burdened); c, cm = cs.count(mk); b, bm = cs.count(mk & burdened) grad.append({"band": lab, "label": BAND_LABELS[lab], "households": c, "households_moe90": cm, "burdened_share": p, "burdened_share_moe90": pm, "burdened_households": b, "burdened_households_moe90": bm, "n_records": int(mk.sum())}) put(f"grad_{lab}", f"Burdened share, income {BAND_LABELS[lab]}", p, pm, mk.sum(), "share") res["gradient"] = grad total_b = sum(g["burdened_households"] for g in grad) res["burdened_share_by_band"] = {g["band"]: g["burdened_households"] / total_b for g in grad} # Trait columns (census side). trait_num = { "manhattan": (h["borough"] == "Manhattan").to_numpy(), "recent_mover": np.isin(mv, [1, 2]), "long_tenure": mv >= 5, "lives_alone": h["alone"].to_numpy(), "married_couple": (h["HHT"] == 1).to_numpy(), "has_children": (h["R18"] == 1).to_numpy(), "has_100k_earner": (h["max_pinc"] >= HIGH_INCOME).to_numpy(), "two_plus_earners": (h["n_earners"] >= 2).to_numpy(), "head_under_35": (h["HHLDRAGEP"] < 35).to_numpy(), "head_65_plus": (h["HHLDRAGEP"] >= 65).to_numpy(), "rent_3000_plus": (h["GRNTP"] >= 3000).to_numpy(), "old_multi": h["old_multi"].to_numpy(), } # The two groups. res["groups"] = {} hb_total = cs.count(hi & burdened)[0] for key, label in GROUPS: g = groups[key] c, cm = cs.count(g); b, bm = cs.count(g & burdened); p, pm = cs.share(g, burdened); s, sm = cs.share(g, severe) d = {"label": label, "households": c, "households_moe90": cm, "share_of_100k": c / res["lines"]["100k"]["households"], "burdened_share": p, "burdened_share_moe90": pm, "burdened_households": b, "burdened_households_moe90": bm, "share_of_burdened_100k": b / hb_total, "severe_share": s, "severe_share_moe90": sm, "n_records": int(g.sum()), "n_burdened_records": int((g & burdened).sum()), "medians": {}, "profile": []} for grp, mk in (("burdened", g & burdened), ("unburdened", g & ~burdened)): d["medians"][grp] = {"household_income": cs.median(mk, "HINCP"), "gross_rent": cs.median(mk, "GRNTP"), "burden_pct": cs.median(mk, "GRPIP"), "max_person_income": cs.median(mk, "max_pinc"), "household_size": cs.median(mk, "NP"), "n_records": int(mk.sum())} for tkey, tlabel in TRAITS: pb, mb = cs.share(g & burdened, trait_num[tkey]); pu, mu = cs.share(g & ~burdened, trait_num[tkey]) d["profile"].append({"trait": tkey, "label": tlabel, "burdened": pb, "burdened_moe90": mb, "unburdened": pu, "unburdened_moe90": mu, "distinguishable": bool(abs(pb - pu) > mb + mu)}) put(f"{key}_share", f"Burdened share, {label} $100k+", p, pm, g.sum(), "share") put(f"{key}_count", f"Burdened households, {label} $100k+", b, bm, (g & burdened).sum(), "count") res["groups"][key] = d res["profile"] = res["groups"]["families"]["profile"] # Move-in staircase for families and singles, with who lives there. fam = groups["families"] res["move_in"] = [] for code, label in MV_BANDS: mk = fam & (mv == code) p, pm = cs.share(mk, burdened); s, _ = cs.share(fam, mv == code); om, omm = cs.share(mk, h["old_multi"].to_numpy()) res["move_in"].append({ "mv": code, "label": label, "share_of_group": s, "burdened_share": p, "burdened_share_moe90": pm, "median_gross_rent": cs.median(mk, "GRNTP"), "median_income": cs.median(mk, "HINCP"), "median_people": cs.median(mk, "NP"), "mean_people": cs.mean(mk, "NP"), "share_alone": cs.share(mk, h["alone"].to_numpy())[0], "share_with_children": cs.share(mk, (h["R18"] == 1).to_numpy())[0], "median_bedrooms": cs.median(mk, "BDSP"), "share_0_1_bedroom": cs.share(mk, (h["BDSP"] <= 1).to_numpy())[0], "share_3plus_bedroom": cs.share(mk, (h["BDSP"] >= 3).to_numpy())[0], "median_head_age": cs.median(mk, "HHLDRAGEP"), "share_head_65_plus": cs.share(mk, (h["HHLDRAGEP"] >= 65).to_numpy())[0], "old_multi_share": om, "old_multi_moe90": omm, "n_records": int(mk.sum())}) put(f"mv_{code}", f"Burdened share, families/singles $100k+, moved in {label}", p, pm, mk.sum(), "share") man = (h["borough"] == "Manhattan").to_numpy() res["move_in_manhattan"] = [{"mv": code, "label": label, "burdened_share": cs.share(fam & man & (mv == code), burdened)[0], "burdened_share_moe90": cs.share(fam & man & (mv == code), burdened)[1], "median_gross_rent": cs.median(fam & man & (mv == code), "GRNTP"), "n_records": int((fam & man & (mv == code)).sum())} for code, label in MV_BANDS] room = groups["roommates"] res["move_in_roommates"] = [{"mv": code, "label": label, "share_of_group": cs.share(room, mv == code)[0], "burdened_share": cs.share(room & (mv == code), burdened)[0], "burdened_share_moe90": cs.share(room & (mv == code), burdened)[1], "median_gross_rent": cs.median(room & (mv == code), "GRNTP"), "median_people": cs.median(room & (mv == code), "NP"), "n_records": int((room & (mv == code)).sum())} for code, label in MV_BANDS] # Families vs roommates, settled vs recent. rec = np.isin(mv, [1, 2]) res["recent_x_group"] = {} for key, _ in GROUPS: for name, mk in (("settled", groups[key] & ~rec), ("recent", groups[key] & rec)): p, pm = cs.share(mk, burdened) res["recent_x_group"][f"{key}_{name}"] = {"burdened_share": p, "burdened_share_moe90": pm, "n_records": int(mk.sum())} put(f"{key}_{name}", f"Burdened share, {key} {name}", p, pm, mk.sum(), "share") # Owners, same yardstick. owner_b = (h["OCPIP"] >= BURDEN).to_numpy() own_all = ((h["in_nyc"] == 1) & (h["TEN"] == 1) & (h["NP"] > 0) & h["OCPIP"].notna()).to_numpy() res["owners"] = [] for key, label, lo, hi_ in OWNER_BANDS: mo = own_all & (inc >= lo) & (inc < hi_); mr = nyc & (inc >= lo) & (inc < hi_) po, pmo = cs.share(mo, owner_b); pr, pmr = cs.share(mr, burdened) res["owners"].append({"band": key, "label": label, "owner_burdened_share": po, "owner_moe90": pmo, "owner_households": cs.count(mo)[0], "owner_burdened_households": cs.count(mo & owner_b)[0], "renter_burdened_share": pr, "renter_moe90": pmr, "renter_households": cs.count(mr)[0], "n_owner_records": int(mo.sum()), "n_renter_records": int(mr.sum())}) put(f"owner_{key}", f"Cost-burdened share, owners with a mortgage {label}", po, pmo, mo.sum(), "share") mo = own_all & (inc >= HIGH_INCOME) po, pmo = cs.share(mo, owner_b) res["owners_100k_plus"] = {"burdened_share": po, "moe90": pmo, "households": cs.count(mo)[0], "burdened_households": cs.count(mo & owner_b)[0], "n_records": int(mo.sum())} put("owner_100k", "Cost-burdened share, owners with a mortgage $100k+", po, pmo, mo.sum(), "share") # What burdened means: residual income for burdened families/singles. fb = fam & burdened med_resid = cs.median(fb, "residual") p, pm = cs.share(nyc, h["residual"].to_numpy() < med_resid) res["residual_income"] = {"burdened_families_median_after_rent": med_resid, "share_of_nyc_renters_with_less": p, "moe90": pm} put("resid_share", "Share of NYC renters with less left after rent than the burdened families/singles median", p, pm, nyc.sum(), "share") # Near the line, families/singles. ratio = h["ratio"].to_numpy() p, pm = cs.share(fam, (ratio >= 25) & (ratio < 30)) res["near_the_line"] = {"share_25_to_30": p, "moe90": pm} for pct in (10, 20): p2, pm2 = cs.share(fam, ratio * (1 + pct / 100) >= BURDEN) res["near_the_line"][f"burdened_if_rent_up_{pct}pct"] = p2 res["near_the_line"][f"burdened_if_rent_up_{pct}pct_moe90"] = pm2 # Sanity check on the stabilised stock proxy. p, pm = cs.share(nyc, h["old_multi"].to_numpy()) long_all = nyc & (mv == 7) res["stabilised_stock_check"] = { "share_of_all_nyc_renters_in_old_multi": p, "moe90": pm, "all_income_30plus_years_median_rent": cs.median(long_all, "GRNTP"), "all_income_30plus_years_old_multi_share": cs.share(long_all, h["old_multi"].to_numpy())[0], "hvs_2023": {"citywide_median_rent": 1641, "stabilised_median_rent": 1500, "market_median_rent": 2000, "controlled_median_rent": 988, "note": "2023 NYCHVS Selected Initial Findings, 2023 dollars"}} return res, pd.DataFrame(rows) # --------------------------------------------------------------- synthetic side def load_synthetic(conn: sqlite3.Connection) -> pd.DataFrame: hh = sdb.read_frame(conn, """ SELECT h.household_id, h.puma_geoid, h.borough, h.tenure, h.rent_regulation, h.household_size, h.n_children, h.household_type, h.household_income, h.monthly_rent, h.rent_burden_pct, h.owner_cost_pct, h.year_moved_in, h.bedrooms, h.vehicles, h.year_built, h.building_type, b.liquid_assets, b.net_worth, b.scf_reference_group, p.age AS head_age, p.education_band AS head_education, pu.puma_name FROM household h LEFT JOIN balance_sheet b ON b.household_id = h.household_id LEFT JOIN person p ON p.household_id = h.household_id AND p.relationship_to_head = 'Reference person' LEFT JOIN puma pu ON pu.puma_geoid = h.puma_geoid WHERE h.in_nyc = 1 """).drop_duplicates("household_id").reset_index(drop=True) ppl = sdb.read_frame(conn, """ SELECT p.household_id, MAX(COALESCE(p.personal_income, 0)) AS max_person_income, SUM(CASE WHEN p.age >= 18 THEN 1 ELSE 0 END) AS n_adults, SUM(CASE WHEN COALESCE(p.personal_income, 0) >= 1000 THEN 1 ELSE 0 END) AS n_earners FROM person p JOIN household h ON h.household_id = p.household_id WHERE h.in_nyc = 1 GROUP BY p.household_id """) hh = hh.merge(ppl, on="household_id", how="left") for col in ("household_size", "n_children", "vehicles", "bedrooms", "n_adults", "n_earners"): hh[col] = pd.to_numeric(hh[col], errors="coerce") ht = hh["household_type"].fillna("") hh["roommates"] = ht.str.contains("Not living alone") hh["lives_alone"] = ht.str.contains("Living alone") hh["married_couple"] = ht.str.startswith("Married") hh["old_multi"] = hh["year_built"].isin(OLD_YEARS) & hh["building_type"].isin(BIG_BUILDINGS) hh["mv"] = hh["year_moved_in"].map({lab: code for code, lab in MV_BANDS}) return hh def load_household_view(conn: sqlite3.Connection) -> pd.DataFrame: import verify as V return V.load_household_view(conn) def bootstrap(values, fn, rng, reps: int = REPS, ci: float = CI) -> tuple[float, float, float]: values = np.asarray(values, dtype=float) values = values[~np.isnan(values)] if len(values) == 0: return float("nan"), float("nan"), float("nan") point = float(fn(values)) idx = rng.integers(0, len(values), size=(reps, len(values))) stats = np.array([fn(values[i]) for i in idx]) lo, hi = np.quantile(stats, [(1 - ci) / 2, 1 - (1 - ci) / 2]) return point, float(lo), float(hi) def share_ci(mask, rng) -> dict: p, lo, hi = bootstrap(np.asarray(mask, dtype=float), np.mean, rng) return {"share": p, "ci_lo": lo, "ci_hi": hi} def median_ci(values, rng) -> dict: p, lo, hi = bootstrap(values, np.median, rng) return {"median": p, "ci_lo": lo, "ci_hi": hi} def reproduction_table(conn: sqlite3.Connection, registry: SeriesRegistry) -> pd.DataFrame: view = load_household_view(conn) deff = float(conn.execute("SELECT design_effect FROM build_run ORDER BY run_id DESC LIMIT 1").fetchone()[0]) cs = C.build_household_set(registry, BASE_YEAR, n=len(view), design_effect=deff) wanted = [c for c in cs.all if c.id.startswith(("hh.renter_burden", "hh.renter_severe", "hh.renter_share", "hh.median_gross_rent")) and ".by_area." not in c.id] rows = [] for c in wanted: got = c.statistic(view, None); band = c.effective_band(view); v = registry[c.id].at(BASE_YEAR) rows.append({"target_id": c.id, "description": c.description, "published": float(c.target), "published_moe90": float(v.moe) if v.moe is not None else None, "synthetic": float(got), "diff": float(got - c.target), "band": float(band), "diff_over_band": float(abs(got - c.target) / band) if band else None, "within_band": bool(abs(got - c.target) <= band), "provenance": c.provenance, "source_id": c.source_id}) return pd.DataFrame(rows) def synthetic_estimates(hh: pd.DataFrame, census: dict, rng: np.random.Generator) -> tuple[dict, pd.DataFrame]: r = hh[(hh["tenure"] == "Rented") & hh["rent_burden_pct"].notna() & (hh["household_income"] > 0)].copy() r["burdened"] = r["rent_burden_pct"] >= BURDEN r["severe"] = r["rent_burden_pct"] >= SEVERE r["residual"] = r["household_income"] - 12 * r["monthly_rent"] r["ratio"] = 12 * r["monthly_rent"] / r["household_income"] * 100 inc = r["household_income"].to_numpy(dtype=float) scale = census["nyc_renter_households"]["value"] / len(r) res: dict = {"n_nyc_renters": int(len(r)), "real_households_per_synthetic": scale, "bootstrap": {"reps": REPS, "interval": CI}} res["lines"] = {} for line, key in ((100_000, "100k"), (150_000, "150k"), (200_000, "200k")): sub = r[inc >= line]; s = share_ci(sub["burdened"], rng) res["lines"][key] = {"n_synthetic": int(len(sub)), "n_burdened_synthetic": int(sub["burdened"].sum()), "households_scaled": len(sub) * scale, "burdened_households_scaled": sub["burdened"].sum() * scale, "burdened_share": s["share"], "ci_lo": s["ci_lo"], "ci_hi": s["ci_hi"], "severe_share": float(sub["severe"].mean())} res["gradient"] = [] for lo, hi_, lab in BANDS: sub = r[(inc >= lo) & (inc < hi_)]; s = share_ci(sub["burdened"], rng) res["gradient"].append({"band": lab, "label": BAND_LABELS[lab], "n_synthetic": int(len(sub)), "burdened_share": s["share"], "ci_lo": s["ci_lo"], "ci_hi": s["ci_hi"]}) hi = r[inc >= HIGH_INCOME].copy() hi["recent_mover"] = hi["mv"].isin([1, 2]) hi["long_tenure"] = hi["mv"] >= 5 hi["manhattan"] = hi["borough"] == "Manhattan" hi["has_100k_earner"] = hi["max_person_income"] >= HIGH_INCOME hi["two_plus_earners"] = hi["n_earners"] >= 2 hi["has_children"] = hi["n_children"].fillna(0) > 0 hi["head_under_35"] = hi["head_age"] < 35 hi["head_65_plus"] = hi["head_age"] >= 65 hi["rent_3000_plus"] = hi["monthly_rent"] >= 3000 groups = {"families": ~hi["roommates"].to_numpy(), "roommates": hi["roommates"].to_numpy()} b = hi["burdened"].to_numpy() res["groups"] = {} for key, label in GROUPS: g = groups[key]; sub = hi[g]; s = share_ci(sub["burdened"], rng) d = {"label": label, "n_synthetic": int(g.sum()), "n_burdened_synthetic": int((g & b).sum()), "households_scaled": g.sum() * scale, "burdened_households_scaled": (g & b).sum() * scale, "share_of_100k": float(g.mean()), "burdened_share": s["share"], "ci_lo": s["ci_lo"], "ci_hi": s["ci_hi"], "share_of_burdened_100k": float((g & b).sum() / b.sum()), "severe_share": float(sub["severe"].mean()), "medians": {}, "profile": []} for grp, mk in (("burdened", g & b), ("unburdened", g & ~b)): s2 = hi[mk] d["medians"][grp] = {"n_synthetic": int(mk.sum()), "household_income": median_ci(s2["household_income"], rng), "gross_rent": median_ci(s2["monthly_rent"], rng), "burden_pct": median_ci(s2["rent_burden_pct"], rng), "max_person_income": median_ci(s2["max_person_income"], rng), "household_size": median_ci(s2["household_size"], rng)} for tkey, tlabel in TRAITS: sb = share_ci(hi.loc[g & b, tkey], rng); su = share_ci(hi.loc[g & ~b, tkey], rng) d["profile"].append({"trait": tkey, "label": tlabel, "burdened": sb["share"], "burdened_lo": sb["ci_lo"], "burdened_hi": sb["ci_hi"], "unburdened": su["share"], "unburdened_lo": su["ci_lo"], "unburdened_hi": su["ci_hi"]}) res["groups"][key] = d res["profile"] = res["groups"]["families"]["profile"] fam = groups["families"] res["move_in"] = [] for code, label in MV_BANDS: mk = fam & (hi["mv"] == code).to_numpy(); sub = hi[mk]; s = share_ci(sub["burdened"], rng) res["move_in"].append({"mv": code, "label": label, "n_synthetic": int(mk.sum()), "share_of_group": float(mk.sum() / fam.sum()), "burdened_share": s["share"], "ci_lo": s["ci_lo"], "ci_hi": s["ci_hi"], "median_gross_rent": float(sub["monthly_rent"].median()), "median_income": float(sub["household_income"].median()), "median_people": float(sub["household_size"].median()), "mean_people": float(sub["household_size"].mean()), "share_alone": float(sub["lives_alone"].mean()), "share_with_children": float(sub["has_children"].mean()), "median_bedrooms": float(sub["bedrooms"].median()), "share_0_1_bedroom": float((sub["bedrooms"] <= 1).mean()), "median_head_age": float(sub["head_age"].median()), "share_head_65_plus": float(sub["head_65_plus"].mean()), "old_multi_share": float(sub["old_multi"].mean())}) room = groups["roommates"] res["move_in_roommates"] = [] for code, label in MV_BANDS: mk = room & (hi["mv"] == code).to_numpy(); sub = hi[mk]; s = share_ci(sub["burdened"], rng) res["move_in_roommates"].append({"mv": code, "label": label, "n_synthetic": int(mk.sum()), "burdened_share": s["share"], "ci_lo": s["ci_lo"], "ci_hi": s["ci_hi"], "median_gross_rent": float(sub["monthly_rent"].median()), "median_people": float(sub["household_size"].median())}) rec = hi["recent_mover"].to_numpy() res["recent_x_group"] = {} for key, _ in GROUPS: for name, mk in (("settled", groups[key] & ~rec), ("recent", groups[key] & rec)): res["recent_x_group"][f"{key}_{name}"] = {**share_ci(hi.loc[mk, "burdened"], rng), "n_synthetic": int(mk.sum())} res["burdened_families_top_areas"] = hi.loc[fam & b, "puma_name"].value_counts().head(10).to_dict() res["burdened_families_by_borough"] = hi.loc[fam & b, "borough"].value_counts(normalize=True).round(4).to_dict() own = hh[(hh["tenure"] == "Owned with a mortgage") & hh["owner_cost_pct"].notna() & (hh["household_income"] > 0)] res["owners"] = [] for key, label, lo, hi_ in OWNER_BANDS: so = own[(own["household_income"] >= lo) & (own["household_income"] < hi_)] sr = r[(inc >= lo) & (inc < hi_)] o = share_ci(so["owner_cost_pct"] >= BURDEN, rng); rr = share_ci(sr["burdened"], rng) res["owners"].append({"band": key, "label": label, "n_owner_synthetic": int(len(so)), "owner_burdened_share": o["share"], "owner_lo": o["ci_lo"], "owner_hi": o["ci_hi"], "renter_burdened_share": rr["share"], "renter_lo": rr["ci_lo"], "renter_hi": rr["ci_hi"]}) so = own[own["household_income"] >= HIGH_INCOME] res["owners_100k_plus"] = {**share_ci(so["owner_cost_pct"] >= BURDEN, rng), "n_synthetic": int(len(so))} fb = fam & b med_resid = float(hi.loc[fb, "residual"].median()) res["residual_income"] = {"burdened_families_median_after_rent": med_resid, "share_of_nyc_renters_with_less": float((r["residual"] < med_resid).mean())} ratio = hi.loc[fam, "ratio"] res["near_the_line"] = {"share_25_to_30": float(((ratio >= 25) & (ratio < 30)).mean()), "burdened_if_rent_up_10pct": float((ratio * 1.1 >= BURDEN).mean()), "burdened_if_rent_up_20pct": float((ratio * 1.2 >= BURDEN).mean())} res["stabilised_stock_check"] = {"share_of_all_nyc_renters_in_old_multi": float(r["old_multi"].mean()), "all_income_30plus_years_median_rent": float(r.loc[r["mv"] == 7, "monthly_rent"].median()), "all_income_30plus_years_old_multi_share": float(r.loc[r["mv"] == 7, "old_multi"].mean())} hi["months_of_rent"] = hi["liquid_assets"].clip(lower=0) / hi["monthly_rent"] nia = {"why": "Not reported in the article: the wealth layer maps metro income rank onto national SCF income bands, " "most of this group drew from fallback cells, and months-of-rent is a rent denominator effect. " "See the module docstring, section 8."} for grp, mk in (("burdened", fb), ("unburdened", fam & ~b)): sub = hi[mk] nia[grp] = {"liquid_assets_median": median_ci(sub["liquid_assets"], rng), "net_worth_median": median_ci(sub["net_worth"], rng), "months_of_rent_median": median_ci(sub["months_of_rent"], rng), "scf_cells_used": sub["scf_reference_group"].value_counts().head(6).to_dict()} res["not_in_article"] = nia return res, hi # ------------------------------------------------------------------------ figures def _mpl(): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"svg.fonttype": "none", "font.family": "sans-serif", "font.sans-serif": ["Inter", "Segoe UI", "Helvetica Neue", "Arial", "DejaVu Sans"], "font.size": 11}) return plt def _save(fig, path: Path): """SVG for the site plus a PDF twin; a .pdf path writes the PDF only.""" if path.suffix == ".pdf": fig.savefig(path, format="pdf", transparent=True) return fig.savefig(path, format="svg", transparent=True) pdf_dir = path.parent / "pdf" pdf_dir.mkdir(exist_ok=True) fig.savefig(pdf_dir / path.with_suffix(".pdf").name, format="pdf", transparent=True) def _style(ax, title: str, subtitle: str | None = None, xlabel: str | None = None): ax.set_facecolor("none") for s in ("top", "right", "left"): ax.spines[s].set_visible(False) ax.spines["bottom"].set_color(GRID) ax.tick_params(colors=INK, labelsize=10.5, length=0) ax.yaxis.grid(False) ax.xaxis.grid(True, color=GRID, linewidth=1) ax.set_axisbelow(True) ax.set_title(title, loc="left", fontsize=13, color=INK, fontweight="bold", pad=18) if subtitle: ax.text(0, 1.03, subtitle, transform=ax.transAxes, fontsize=10, color=INK, va="bottom") if xlabel: ax.set_xlabel(xlabel, color=INK, fontsize=10.5) def fig_where(census: dict, path: Path): plt = _mpl() g = census["gradient"] fig, ax = plt.subplots(figsize=(FIG_W, 3.9), dpi=100) fig.patch.set_alpha(0) y = np.arange(len(g))[::-1] total = np.array([r["households"] for r in g]) / 1000 burd = np.array([r["burdened_households"] for r in g]) / 1000 ax.barh(y, total, height=0.55, color=GRID, label="All renter households") ax.barh(y, burd, height=0.55, color=BLUE, label="Rent-burdened (30%+ of income)") ax.set_yticks(y); ax.set_yticklabels([r["label"] for r in g]) for yi, bv, tv in zip(y, burd, total): ax.text(tv + 8, yi, f"{bv:,.0f}k of {tv:,.0f}k", va="center", fontsize=10, color=INK) ax.set_xlim(0, total.max() * 1.42); ax.set_xticks([0, 200, 400, 600]) _style(ax, "Where the burdened renters are", "NYC renter households by income, thousands (census, 2020–24)") ax.legend(loc="lower right", frameon=False, fontsize=10, labelcolor=INK) fig.tight_layout(); _save(fig, path); plt.close(fig) def fig_move_in(census: dict, synth: dict, path: Path, figsize: tuple[float, float] = (FIG_W, 3.9)): """The staircase: burdened share of six-figure families and singles by move-in period, census with margin and synthetic point, median rent written at the right.""" plt = _mpl() c = census["move_in"]; s = {r["mv"]: r for r in synth["move_in"]} fig, ax = plt.subplots(figsize=figsize, dpi=100) fig.patch.set_alpha(0) y = np.arange(len(c))[::-1] p = np.array([r["burdened_share"] for r in c]) * 100 m = np.array([r["burdened_share_moe90"] for r in c]) * 100 sp = np.array([s[r["mv"]]["burdened_share"] for r in c]) * 100 rent = np.array([r["median_gross_rent"] for r in c]) xmax = 30 ax.hlines(y, p - m, p + m, color=ORANGE, linewidth=2) ax.scatter(p, y, s=64, color=ORANGE, zorder=3, edgecolor="white", linewidth=2, label="Census (90% margin)") ax.scatter(sp, y, s=40, color=BLUE, zorder=4, edgecolor="white", linewidth=1.5, label="Synthetic population") for yi, rv in zip(y, rent): ax.text(xmax - 0.4, yi, f"${rv:,.0f}", va="center", ha="right", fontsize=10, color=INK) ax.text(xmax - 0.4, y[0] + 0.85, "median rent", va="bottom", ha="right", fontsize=9, color=INK, style="italic") ax.set_yticks(y); ax.set_yticklabels([MV_SHORT[r["mv"]] for r in c]) ax.set_xlim(0, xmax); ax.set_xticks([0, 10, 20]); ax.set_xticklabels(["0%", "10%", "20%"]) _style(ax, "The burden is a move-in date", "Six-figure families and singles renting in NYC, by how long ago they moved in", xlabel="Rent-burdened (%)") h_in = figsize[1] fig.legend(*ax.get_legend_handles_labels(), loc="lower center", ncol=2, frameon=False, fontsize=10, labelcolor=INK, bbox_to_anchor=(0.5, 0.0)) fig.subplots_adjust(left=0.13 * FIG_W / figsize[0], right=0.97, top=1 - 0.6 / h_in, bottom=0.9 / h_in) _save(fig, path); plt.close(fig) def fig_groups(census: dict, synth: dict, path: Path): """Four bars: families/singles and roommate households, settled and recent.""" plt = _mpl() order = [("families_settled", "Families and singles, settled"), ("families_recent", "Families and singles, recent"), ("roommates_settled", "Roommates, settled"), ("roommates_recent", "Roommates, recent")] c = census["recent_x_group"]; s = synth["recent_x_group"] fig, ax = plt.subplots(figsize=(FIG_W, 3.2), dpi=100) fig.patch.set_alpha(0) y = np.arange(len(order))[::-1] p = np.array([c[k]["burdened_share"] for k, _ in order]) * 100 m = np.array([c[k]["burdened_share_moe90"] for k, _ in order]) * 100 sp = np.array([s[k]["share"] for k, _ in order]) * 100 ax.barh(y, p, height=0.55, color=ORANGE, label="Census") ax.errorbar(p, y, xerr=m, fmt="none", ecolor="white", elinewidth=1.5, capsize=0) ax.scatter(sp, y, s=40, color=BLUE, zorder=4, edgecolor="white", linewidth=1.5, label="Synthetic population") for yi, pv in zip(y, p): ax.text(pv + 0.5, yi, f"{pv:.0f}%", va="center", fontsize=10, color=INK) ax.set_yticks(y); ax.set_yticklabels([lab for _, lab in order]) ax.set_xlim(0, 32); ax.set_xticks([0, 10, 20, 30]); ax.set_xticklabels(["0%", "10%", "20%", "30%"]) _style(ax, "Who is over the line", "Six-figure NYC renters who are rent-burdened (recent = moved in within 2 years)", xlabel="Rent-burdened (%)") ax.legend(loc="upper right", frameon=False, fontsize=10, labelcolor=INK) fig.tight_layout(); _save(fig, path); plt.close(fig) def fig_profile(census: dict, path: Path): plt = _mpl() rows = census["profile"] order = np.argsort([-(r["burdened"] - r["unburdened"]) for r in rows]); rows = [rows[i] for i in order] fig, ax = plt.subplots(figsize=(FIG_W, 5.6), dpi=100) fig.patch.set_alpha(0) y = np.arange(len(rows))[::-1] b = np.array([r["burdened"] for r in rows]) * 100; bm = np.array([r["burdened_moe90"] for r in rows]) * 100 u = np.array([r["unburdened"] for r in rows]) * 100; um = np.array([r["unburdened_moe90"] for r in rows]) * 100 ax.hlines(y, u, b, color=GRID, linewidth=2, zorder=1) ax.hlines(y, b - bm, b + bm, color=BLUE, linewidth=4, alpha=0.35, zorder=2) ax.hlines(y, u - um, u + um, color=ORANGE, linewidth=4, alpha=0.35, zorder=2) ax.scatter(u, y, s=64, color=ORANGE, zorder=3, edgecolor="white", linewidth=2, label="Not burdened") ax.scatter(b, y, s=64, color=BLUE, zorder=3, edgecolor="white", linewidth=2, label="Rent-burdened") ax.set_yticks(y); ax.set_yticklabels([r["label"] + ("" if r["distinguishable"] else " *") for r in rows]) ax.set_xlim(0, 100); ax.set_xticks([0, 25, 50, 75, 100]); ax.set_xticklabels(["0%", "25%", "50%", "75%", "100%"]) _style(ax, "Who the burdened six-figure family or single is", "Share of $100k+ NYC renter households with each trait (census)") ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.1), ncol=2, frameon=False, fontsize=10, labelcolor=INK) fig.text(0, 0.0, "* margins overlap: the data cannot tell these two groups apart on this trait.", fontsize=9.5, color=INK, va="bottom") fig.tight_layout(rect=(0, 0.04, 1, 1)); _save(fig, path); plt.close(fig) def fig_owners(census: dict, path: Path, figsize: tuple[float, float] = (FIG_W, 3.4)): plt = _mpl() rows = census["owners"] fig, ax = plt.subplots(figsize=figsize, dpi=100) fig.patch.set_alpha(0) y = np.arange(len(rows))[::-1] o = np.array([r["owner_burdened_share"] for r in rows]) * 100; om = np.array([r["owner_moe90"] for r in rows]) * 100 r_ = np.array([r["renter_burdened_share"] for r in rows]) * 100; rm = np.array([r["renter_moe90"] for r in rows]) * 100 h = 0.32 ax.barh(y + h / 2 + 0.02, o, height=h, color=ORANGE, label="Owners with a mortgage") ax.barh(y - h / 2 - 0.02, r_, height=h, color=BLUE, label="Renters") ax.errorbar(o, y + h / 2 + 0.02, xerr=om, fmt="none", ecolor="white", elinewidth=1.5, capsize=0) ax.errorbar(r_, y - h / 2 - 0.02, xerr=rm, fmt="none", ecolor="white", elinewidth=1.5, capsize=0) for yi, ov, rv in zip(y, o, r_): ax.text(ov + 1.2, yi + h / 2 + 0.02, f"{ov:.0f}%", va="center", fontsize=10, color=INK) ax.text(rv + 1.2, yi - h / 2 - 0.02, f"{rv:.0f}%", va="center", fontsize=10, color=INK) ax.set_yticks(y); ax.set_yticklabels([r["label"] for r in rows]) ax.set_xlim(0, 66); ax.set_xticks([0, 20, 40, 60]); ax.set_xticklabels(["0%", "20%", "40%", "60%"]) _style(ax, "Same income, same yardstick: owners are more burdened", "Share paying 30%+ of income in housing costs, NYC households (census)") ax.legend(loc="lower right", frameon=False, fontsize=10, labelcolor=INK) fig.tight_layout(); _save(fig, path); plt.close(fig) # ------------------------------------------------------------------- companion data def companion_dataset(hi: pd.DataFrame, rng: np.random.Generator) -> pd.DataFrame: cols = ["borough", "household_income", "monthly_rent", "rent_burden_pct", "burdened", "roommates", "household_type", "household_size", "n_adults", "n_earners", "max_person_income", "n_children", "head_age", "head_education", "year_moved_in", "bedrooms", "year_built", "building_type", "old_multi", "rent_regulation", "vehicles"] df = hi[cols].copy().rename(columns={"roommates": "roommate_household", "old_multi": "pre1970_building_5plus_units"}) df.insert(0, "row_id", [f"sfr-{i:05d}" for i in rng.permutation(len(df)) + 1]) df = df.sort_values("row_id").reset_index(drop=True) df["is_synthetic"] = 1 return df # ------------------------------------------------------------------------- main def main(argv: list[str] | None = None) -> int: ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--seed", type=int, default=20240101) ap.add_argument("--db", type=Path, default=DB_PATH) args = ap.parse_args(argv) rng = np.random.default_rng(args.seed) OUT.mkdir(parents=True, exist_ok=True) print("census side: PUMS with replicate weights ...") h, rw = load_pums() census, pums_rows = census_estimates(Census(h, rw)) conn = sqlite3.connect(args.db) registry = SeriesRegistry.load(TARGETS) run = conn.execute("SELECT run_id, seed, n_people, n_households, tool_version, started_utc " "FROM build_run ORDER BY run_id DESC LIMIT 1").fetchone() print("reproduction table ...") repro = reproduction_table(conn, registry) print(repro[["target_id", "published", "synthetic", "band", "within_band"]].to_string(index=False)) print("synthetic side: release households ...") hh = load_synthetic(conn) synth, hi = synthetic_estimates(hh, census, rng) res = {"article": SLUG, "definitions": {"high_income_line": HIGH_INCOME, "burden_threshold_pct": BURDEN, "severe_threshold_pct": SEVERE, "dollar_year": BASE_YEAR, "groups": dict(GROUPS), "stabilised_stock_proxy": "buildings with 5+ units (BLD >= 6) built before 1970 (YRBLT <= 1960)", "census_margin": "90% (successive difference replication, 80 replicate weights)", "synthetic_interval": f"{int(CI * 100)}% percentile bootstrap, {REPS} household resamples"}, "census": census, "synthetic": synth, "reproduction": repro.to_dict(orient="records"), "build": {"run_id": run[0], "seed": run[1], "n_people": run[2], "n_households": run[3], "tool_version": run[4], "started_utc": run[5]}, "bootstrap_seed": args.seed} print("figures ...") fig_where(census, OUT / "fig1_where_the_burdened_are.svg") fig_move_in(census, synth, OUT / "fig2_move_in_date.svg") fig_groups(census, synth, OUT / "fig3_families_and_roommates.svg") fig_owners(census, OUT / "fig4_owners_vs_renters.svg") fig_profile(census, OUT / "figS1_profile.svg") companion_dataset(hi, rng).to_csv(OUT / "data.csv", index=False) pums_rows.to_csv(OUT / "pums_direct.csv", index=False) repro.to_csv(OUT / "reproduction.csv", index=False) (OUT / "results.json").write_text(json.dumps(res, indent=1, default=float), encoding="utf-8") for old in ("fig1_burden_by_income.svg", "fig2_where_the_burdened_are.svg", "fig4_months_of_rent.svg", "fig3_profile.svg", "fig3_recent_and_roommates.svg"): for folder in (OUT, SITE_OUT, OUT / "pdf"): (folder / old).unlink(missing_ok=True) (folder / old).with_suffix(".pdf").unlink(missing_ok=True) SITE_OUT.mkdir(parents=True, exist_ok=True) for f in OUT.glob("*"): if f.suffix in (".svg", ".csv", ".json"): shutil.copy2(f, SITE_OUT / f.name) shutil.copy2(Path(__file__), SITE_OUT / "analysis.py") cf = census["groups"]["families"]; sf = synth["groups"]["families"] print(f"\ncensus: six-figure families/singles burdened {cf['burdened_share']:.1%} ± {cf['burdened_share_moe90']:.1%}" f" -> {cf['burdened_households']:,.0f} ± {cf['burdened_households_moe90']:,.0f} households") print(f"synthetic: {sf['burdened_share']:.1%} [{sf['ci_lo']:.1%}, {sf['ci_hi']:.1%}] from {sf['n_burdened_synthetic']} of {sf['n_synthetic']}") print(f"outputs -> {OUT.relative_to(P.root)} and {SITE_OUT.relative_to(P.root)}") return 0 if __name__ == "__main__": raise SystemExit(main())