"""Not just the cold: thirty winters of heat in New York: the analysis behind a SparkyData article. py src/build_article.py not-just-the-cold Was last winter's record in heat complaints just the cold, and how does it look against thirty years of New York winters? Three records answer it, each set beside the others rather than joined (none observes the same household or building as another): - HPD's complaint file (card hpd_heat_complaints): every no-heat complaint to the City since February 2003, by month, borough and building, cross-checked against 311's own record (card nyc311_heat). Public housing is all but absent from it. - NOAA's Central Park record (card noaa_central_park): heating degree days, so each heat season's complaints can be set against how cold it was. - The Housing and Vacancy Survey (cards nychvs_history, nychvs_2023): whether a renter household's heat broke down for six hours or more, every survey from 1991 to 2023, by kind of rental, public housing included. The question changed in 2021, so the series is two stretches, compared within each. The complaints are administrative counts (no sampling margin); the unit is HPD's problem record with a code beginning NO HEAT, which tracks 311's requests within about 2% a season. "The cold alone" is a computed benchmark: the season's degree days times the complaints per degree day of the sixteen winters 2003-04 to 2018-19, before the pandemic. Its noise band comes from how far those sixteen winters strayed from it; two regression models, the equally cold winters, a per-renter-household version and a count without duplicates bound it. Survey shares carry 90% margins: replicate weights from 2011, the Census Bureau's generalized variance formula before (nycdata.hvs_history). Quoted report figures go through A.quotes(), so each one is found in its saved source. Writes results.json, viz.json, charts.json, data.csv and build.json into output/articles/not_just_the_cold/ and site/public/articles/not-just-the-cold/, and copies this script there as analysis.py. Deterministic: no sampling, and SEED is kept for any resampling. """ from __future__ import annotations import numpy as np import pandas as pd from scipy import stats as st from _paths import RAW from article_kit import Article, lines, r from nycdata import city, weather from nycdata import hvs_history as hh from sparkylib import models from sparkylib.stats import moe_difference SLUG = "not-just-the-cold" SEED = 20261002 FIRST, LATEST = 2003, 2025 # complete heat seasons 2003-04 to 2025-26 BASE = list(range(2003, 2019)) # the sixteen winters 2003-04 to 2018-19, before the pandemic BOOT = 20_000 # bootstrap draws for the excess KNOTS = (25, 35, 45, 55) # daily mean temperature (F) where the response to cold may bend NIGHT_KNOTS = (10, 20) # daily minimum temperature (F): the bitter nights Z90 = 1.645 PEER_BAND = 0.05 # "about as cold": within 5% in heating degree days WEATHER_FROM = 1996 # thirty winters, 1996-97 to 2025-26 HIGH = 20 # a high-complaint building: 20 or more in a season LOOKBACK = 10 # seasons searched for a repeat REPEAT_FROM = 2013 # the first season with ten earlier seasons in the record TOP_DROP = 100 # concentration checked without the 100 busiest buildings SEASON_AXIS = "Heat season (October to May)" # the time axis of every chart ticked by season BOROUGHS = city.BOROUGHS MONTHS = [(10, "Oct", "October"), (11, "Nov", "November"), (12, "Dec", "December"), (1, "Jan", "January"), (2, "Feb", "February"), (3, "Mar", "March"), (4, "Apr", "April"), (5, "May", "May")] NYCHA_ADDR = RAW / "housing" / "nycha_residential_addresses.csv" LOTS = RAW / "housing" / "hpd_heat_complaints_by_lot_2025_26.csv" BY_BUILDING = RAW / "housing" / "hpd_heat_complaints_by_building_2025_26.csv" REGISTRY = RAW / "housing" / "hpd_buildings_registry_units.csv" PROV = "NYC HPD, Housing Maintenance Code Complaints and Problems (NYC Open Data)" PROV_HVS = "U.S. Census Bureau and NYC HPD, New York City Housing and Vacancy Survey, 1991-2023" def label(season: int) -> str: return f"{season}-{str(season + 1)[2:]}" def k(x: float) -> str: """A count in thousands, as the prose prints it: 304,578 -> '305,000'.""" return f"{int(round(x, -3)):,}" def pct(x: float) -> str: """A whole percentage, halves rounded up (48.5 -> '49%'), as a reader would round it.""" return f"{int(np.floor(x + 0.5))}%" def num(x: float) -> str: return f"{int(round(x)):,}" # ------------------------------------------------------------------ complaints and the cold def renter_households_by_season(seasons_: list[int]) -> pd.Series: """Renter households outside public housing in each heat season (public-housing tenants report to NYCHA, not to HPD): the survey's weighted counts (every cycle, 1991-2023) interpolated to the January of the season, held at the 2023 count after.""" years = [y for y in hh.YEARS if hh.available(y)] pts = pd.Series({y: hh.renter_households(y, public_housing=False) for y in years}) x = [s + 1 for s in seasons_] return pd.Series(np.interp(x, pts.index.to_numpy(float), pts.to_numpy()), index=seasons_) def seasons() -> pd.DataFrame: """One row per complete heat season: no-heat complaints (and those not flagged as duplicates), every heat and hot-water complaint, heating degree days, renter households, and the benchmarks for the cold.""" p = city.hpd_heat() c = p.dropna(subset=["season"]).groupby("season")["problems"].sum() a = city.hpd_heat(all_categories=True) ca = a.dropna(subset=["season"]).groupby("season")["problems"].sum() det = city.hpd_heat_detail() nd = det[~det["duplicate"]].groupby("season")["problems"].sum() w = weather.heat_seasons() s = pd.DataFrame({"complaints": c, "all_heat_hot_water": ca}) s.index = s.index.astype(int) s = s.loc[FIRST:LATEST].join(w[["hdd", "coldest_f", "nights_10f"]]) s["nondup"] = nd.reindex(s.index) if s[["hdd", "nondup"]].isna().any().any(): raise ValueError("a heat season has no degree days or no duplicate split: refresh the sources") s["renters"] = renter_households_by_season(list(s.index)) base = s.loc[BASE] rate = base["complaints"].sum() / base["hdd"].sum() s["expected"] = rate * s["hdd"] s["ratio_rate"] = s["complaints"] / s["expected"] # noise: how far the base winters strayed from the benchmark, on a log scale sd = float(np.log(base["complaints"] / (rate * base["hdd"])).std(ddof=1)) # per renter household rate_h = base["complaints"].sum() / (base["hdd"] * base["renters"]).sum() s["ratio_household"] = s["complaints"] / (rate_h * s["hdd"] * s["renters"]) # without the complaints flagged as duplicates rate_nd = base["nondup"].sum() / base["hdd"].sum() s["ratio_nondup"] = s["nondup"] / (rate_nd * s["hdd"]) # model A: log complaints on log degree days, fitted on the base seasons X = np.c_[np.ones(len(base)), np.log(base["hdd"])] coef = np.linalg.lstsq(X, np.log(base["complaints"].astype(float)), rcond=None)[0] s["ratio_loglog"] = s["complaints"] / np.exp(coef[0] + coef[1] * np.log(s["hdd"])) # model B: monthly, log complaints on log degree days with a term for each month s["ratio_monthly"] = monthly_model(p) # matched winters: the base seasons within PEER_BAND of the season's degree days peers, lo, hi, n = [], [], [], [] for y in s.index: m = base[(abs(base["hdd"] / s.loc[y, "hdd"] - 1) <= PEER_BAND) & (base.index != y)] peers.append(m["complaints"].mean() if len(m) else np.nan) lo.append(m["complaints"].min() if len(m) else np.nan) hi.append(m["complaints"].max() if len(m) else np.nan) n.append(len(m)) s["peer_mean"], s["peer_min"], s["peer_max"], s["peer_n"] = peers, lo, hi, n s["ratio_peers"] = s["complaints"] / s["peer_mean"] s.attrs.update(rate=rate, sd=sd, elasticity=coef[1], rate_household=rate_h, rate_nondup=rate_nd) return s def monthly_model(p: pd.DataFrame) -> pd.Series: """Season totals over a monthly model fitted on the base seasons: log complaints = month effect + b * log degree days. Returns actual / predicted for each season.""" m = p.groupby("month", as_index=False)["problems"].sum().rename(columns={"problems": "complaints"}) w = weather.monthly() m = m.merge(w[["month", "hdd"]], on="month", how="inner") m["season"] = city.heat_season(m["month"]) m = m.dropna(subset=["season"]) m["season"] = m["season"].astype(int) m = m[(m["season"] >= FIRST) & (m["season"] <= LATEST)] mos = [x[0] for x in MONTHS] def design(x): return np.c_[np.log(x["hdd"].clip(lower=1)), np.column_stack([(x["month"].dt.month == q).astype(float) for q in mos])] fit = m[m["season"].isin(BASE)] coef = np.linalg.lstsq(design(fit), np.log(fit["complaints"].astype(float)), rcond=None)[0] m["pred"] = np.exp(design(m) @ coef) g = m.groupby("season")[["complaints", "pred"]].sum() return g["complaints"] / g["pred"] def complaint_months() -> pd.DataFrame: """No-heat complaints by month of each complete season (October to May).""" p = city.hpd_heat() m = p.dropna(subset=["season"]).groupby("month", as_index=False)["problems"].sum().rename(columns={"problems": "complaints"}) m["season"] = city.heat_season(m["month"]).astype(int) return m[(m["season"] >= FIRST) & (m["season"] <= LATEST)] # ------------------------------------------------------------------ buildings def concentration() -> pd.DataFrame: """Per season: buildings with any no-heat complaint, the busiest tenth's share (with and without the TOP_DROP busiest buildings), how many buildings made up half of all complaints, and buildings with HIGH or more.""" h = city.hpd_heat_buildings() h = h[h["season"].between(FIRST, LATEST)] rows = [] for s, g in h.groupby("season"): n = np.sort(np.repeat(g["n"].to_numpy(), g["buildings"].to_numpy()))[::-1] # one entry per building, largest first def top10(v): c = np.cumsum(v) return 100 * c[int(np.ceil(0.1 * len(v))) - 1] / c[-1] tot, B = n.sum(), len(n) cum = np.cumsum(n) rows.append({"season": int(s), "buildings": B, "complaints": int(tot), "per_building": tot / B, "top10_share": top10(n), "top10_share_wo_top": top10(n[TOP_DROP:]), "half_buildings": int(np.searchsorted(cum, tot / 2) + 1), "high_buildings": int((n >= HIGH).sum()), "high_share": 100 * n[n >= HIGH].sum() / tot, "busiest": int(n[0])}) out = pd.DataFrame(rows).set_index("season") out["half_share_of_buildings"] = 100 * out["half_buildings"] / out["buildings"] return out def repeats(cz: pd.DataFrame) -> pd.DataFrame: """For each season from REPEAT_FROM: of the buildings with HIGH or more complaints, how many had HIGH or more in at least one of the LOOKBACK seasons before; and the share of ALL the season's no-heat complaints sent by such repeat buildings.""" t = city.hpd_heat_buildings_ten_plus() t = t[t["n"] >= HIGH] rows = [] for s in range(REPEAT_FROM, LATEST + 1): now = t[t["season"] == s] before = set(t.loc[t["season"].between(s - LOOKBACK, s - 1), "building_id"]) last = set(t.loc[t["season"] == s - 1, "building_id"]) rep = now[now["building_id"].isin(before)] rows.append({"season": s, "buildings": len(now), "repeat": len(rep), "repeat_pct": 100 * len(rep) / len(now), "repeat_last_pct": 100 * now["building_id"].isin(last).sum() / len(now), "repeat_share_of_all": 100 * rep["n"].sum() / cz.loc[s, "complaints"]}) return pd.DataFrame(rows).set_index("season") def building_sizes() -> dict: """Last season's complaints per building against the apartments HPD's registry gives each building: are the high-complaint buildings just the big ones?""" c = pd.read_csv(BY_BUILDING, dtype={"building_id": str}) g = pd.read_csv(REGISTRY, dtype={"buildingid": str}).rename(columns={"buildingid": "building_id"}) m = c.merge(g[["building_id", "legalclassa", "managementprogram"]], on="building_id", how="left") matched = m["legalclassa"].notna() out = {"buildings": int(len(m)), "matched_pct": 100 * matched.mean(), "complaints_matched_pct": 100 * m.loc[matched, "n"].sum() / m["n"].sum(), "nycha_buildings": int((m["managementprogram"] == "NYCHA").sum()), "nycha_complaints": int(m.loc[m["managementprogram"] == "NYCHA", "n"].sum())} m = m[m["legalclassa"] > 0] hi_, lo_ = m[m["n"] >= HIGH], m[m["n"] < HIGH] out.update({"high_median_apartments": float(hi_["legalclassa"].median()), "other_median_apartments": float(lo_["legalclassa"].median()), "high_per_apartment": hi_["n"].sum() / hi_["legalclassa"].sum(), "other_per_apartment": lo_["n"].sum() / lo_["legalclassa"].sum()}) out["per_apartment_ratio"] = out["high_per_apartment"] / out["other_per_apartment"] return out def borough_complaints() -> pd.DataFrame: p = city.hpd_heat() b = p.dropna(subset=["season"]).groupby(["season", "borough"])["problems"].sum().unstack() b.index = b.index.astype(int) return b.loc[FIRST:LATEST, BOROUGHS] def public_housing_complaints() -> dict: """No-heat complaints in 2025-26 on lots in NYCHA's published residential address list.""" a = pd.read_csv(NYCHA_ADDR, dtype=str) lots = set(a["borough_block_lot"].dropna().str.replace(r"\D", "", regex=True)) c = pd.read_csv(LOTS, dtype={"bbl": str}) c["bbl"] = c["bbl"].fillna("").str.replace(r"\.0$", "", regex=True) on = c[c["bbl"].isin(lots)] return {"complaints": int(on["complaints"].sum()), "lots": int(len(on)), "all": int(c["complaints"].sum()), "nycha_lots": len(lots), "developments": int(a["development"].nunique())} # ------------------------------------------------------------------ the survey def survey() -> tuple[pd.DataFrame, dict]: """Heating breakdowns of six hours or more among renter households, by group and survey, and the comparisons the prose states (each with its 90% range, within one question).""" s = hh.breakdown_series() v = {(int(x.year), x.group): x for x in s.itertuples()} def change(g, a, b): d = v[(b, g)].pct - v[(a, g)].pct m = moe_difference(v[(a, g)].moe, v[(b, g)].moe) return {"from": r(v[(a, g)].pct), "to": r(v[(b, g)].pct), "pts": r(d), "lo": r(d - m), "hi": r(d + m)} ch = { "renters_1991_2017": change("All renters", 1991, 2017), "renters_2002_2005": change("All renters", 2002, 2005), "renters_2008_2011": change("All renters", 2008, 2011), "stabilized_1991_2017": change("Rent stabilized", 1991, 2017), "private_1991_2017": change("Private unregulated", 1991, 2017), "public_2002_2017": change("Public housing", 2002, 2017), "public_2008_2017": change("Public housing", 2008, 2017), "public_1991_2002": change("Public housing", 1991, 2002), "public_1991_2017": change("Public housing", 1991, 2017), "renters_2021_2023": change("All renters", 2021, 2023), "private_2021_2023": change("Private unregulated", 2021, 2023), "stabilized_2021_2023": change("Rent stabilized", 2021, 2023), "public_2021_2023": change("Public housing", 2021, 2023), } def gap(y): """Public housing minus private unregulated within one survey: from the replicate estimates where the file has them (their covariance kept), else the variance formula treating the two as independent, an approximation.""" e_pub, e_priv = hh.breakdown_est(y, "Public housing"), hh.breakdown_est(y, "Private unregulated") if e_pub is not None: g = e_pub - e_priv return g.value, g.moe return (v[(y, "Public housing")].pct - v[(y, "Private unregulated")].pct, moe_difference(v[(y, "Public housing")].moe, v[(y, "Private unregulated")].moe)) gaps = {str(y): {"pts": r(gap(y)[0]), "lo": r(gap(y)[0] - gap(y)[1]), "hi": r(gap(y)[0] + gap(y)[1]), "margin": "replicate" if hh.breakdown_est(y, "Public housing") is not None else "independence approximation"} for y in (1991, 2017, 2021, 2023)} q = hh.breakdown_est(2023, "Public housing") / hh.breakdown_est(2023, "Private unregulated") ch["public_to_private_ratio_2023"] = {"ratio": r(q.value, 2), "lo": r(q.lo, 2), "hi": r(q.hi, 2)} ch["public_minus_private"] = gaps for a_, b_ in ((1991, 2017), (2021, 2023)): (d1, m1), (d2, m2) = gap(a_), gap(b_) dd, mm = d2 - d1, moe_difference(m1, m2) ch[f"gap_change_{a_}_{b_}"] = {"pts": r(dd), "lo": r(dd - mm), "hi": r(dd + mm), "ratio": r(d2 / d1, 2)} return s, ch def daily_model() -> dict: """The weather allowed for day by day: a Poisson model of each day's no-heat complaints on that day's and the two days before's mean temperature and on the night's low and the night before's (linear splines, so complaints may climb more steeply as it gets colder), with terms for the month and the day of the week. Fitted on every heat-season day of the base winters, it predicts each later season's total; each base winter is predicted from a model fitted without it, and the scatter of those predictions gives the reference band.""" c = city.hpd_heat_daily().set_index("day")["problems"] w = weather.daily().set_index("date") days = pd.date_range(f"{FIRST}-09-29", f"{LATEST + 1}-05-31", freq="D") # two days early, for the lags d = pd.DataFrame(index=days) d["y"] = c.reindex(days).fillna(0) d["tmean"] = (w["tmax_f"] + w["tmin_f"]).reindex(days) / 2 d["tmin"] = w["tmin_f"].reindex(days) d["season"] = weather.season_of(pd.Series(days, index=days)) d["tmean1"], d["tmean2"], d["tmin1"] = d["tmean"].shift(1), d["tmean"].shift(2), d["tmin"].shift(1) d = d[d.index >= f"{FIRST}-10-01"].dropna(subset=["season"]) if d[["tmean", "tmean1", "tmean2", "tmin", "tmin1"]].isna().any().any(): raise ValueError("a heat-season day has no temperature in NOAA's daily file") d["season"] = d["season"].astype(int) monthly = city.hpd_heat().dropna(subset=["season"]).groupby("season")["problems"].sum() by_day = d.groupby("season")["y"].sum() if not (by_day == monthly.reindex(by_day.index)).all(): raise ValueError("the daily and monthly complaint files disagree on a season's total") def design(x): cols = [np.ones(len(x))] cols += [(x.index.month == m).astype(float) for m in (10, 11, 12, 2, 3, 4, 5)] cols += [(x.index.dayofweek == k_).astype(float) for k_ in range(1, 7)] for col, knots in (("tmean", KNOTS), ("tmean1", KNOTS), ("tmean2", (35, 50)), ("tmin", NIGHT_KNOTS), ("tmin1", NIGHT_KNOTS)): cols += list(models.hinge_basis(x[col].to_numpy(), knots).T) return np.column_stack(cols) def fit_predict(train, target): tr = d[d["season"].isin(train)] beta = models.poisson_fit(design(tr), tr["y"].to_numpy(float)) te = d[d["season"] == target] return float(te["y"].sum()), float(np.exp(design(te) @ beta).sum()) got = {s_: fit_predict([b for b in BASE if b != s_] if s_ in BASE else BASE, s_) for s_ in range(FIRST, LATEST + 1)} sd = float(np.std(np.log([got[s_][0] / got[s_][1] for s_ in BASE]), ddof=1)) act, pred = got[LATEST] return {"actual": act, "predicted": pred, "excess_pct": 100 * (act / pred - 1), "lo_pct": 100 * (act / (pred * np.exp(Z90 * sd)) - 1), "hi_pct": 100 * (act / (pred * np.exp(-Z90 * sd)) - 1), "sd_log_left_out": sd, "days": int(len(d)), "ratios": {label(s_): a_ / p_ for s_, (a_, p_) in got.items()}} def outcomes(s: pd.DataFrame) -> pd.DataFrame: """How each season's first complaints were closed (HPD's duplicates take the result of the first complaint about the same condition, so they are left out): how many ended in a heat violation, the share of those an inspector checked that did, and how many ended with heat restored or no access. The share of all no-heat complaints is the Comptroller's measure.""" o = city.hpd_heat_outcomes() o = o[~o["duplicate"] & o["season"].between(FIRST, LATEST)] p = o.pivot_table(index="season", columns="outcome", values="problems", aggfunc="sum").fillna(0) for col in ("violation", "violation on file", "no violation", "heat not required", "restored", "no access"): if col not in p: p[col] = 0.0 p.index = p.index.astype(int) n = p.sum(axis=1) out = pd.DataFrame({"first_complaints": n, "violations": p["violation"], "on_file": p["violation on file"], "inspected": p["violation"] + p["violation on file"] + p["no violation"], "heat_not_required": p["heat not required"], "restored_pct": 100 * p["restored"] / n, "no_access_pct": 100 * p["no access"] / n}) out["violation_rate_inspected"] = 100 * (out["violations"] + out["on_file"]) / out["inspected"] out["share_of_all_complaints"] = 100 * out["violations"] / s["complaints"].reindex(out.index) return out def benchmark_checks(s: pd.DataFrame, rng) -> dict: """How firm the excess is: the same benchmark with LaGuardia's and JFK's weather, a regression prediction interval, a bootstrap of the base winters, and any drift within them.""" base = s.loc[BASE] act = float(s.loc[LATEST, "complaints"]) out = {"stations": {}} for stn in ("laguardia", "jfk"): w = weather.heat_seasons(station=stn)["hdd"].reindex(s.index) rate = base["complaints"].sum() / w.loc[BASE].sum() ratio = s["complaints"] / (rate * w) out["stations"][stn] = {"hdd": w.loc[LATEST], "excess_pct": 100 * (ratio.loc[LATEST] - 1), "above_base_range": [label(y) for y in ratio.index if ratio[y] > ratio.loc[BASE].max()]} # log complaints on log degree days over the base winters: a 90% prediction interval for last winter y, x = np.log(base["complaints"].to_numpy(float)), np.log(base["hdd"].to_numpy(float)) X = np.c_[np.ones_like(x), x] beta = np.linalg.lstsq(X, y, rcond=None)[0] res = y - X @ beta x0 = np.array([1.0, np.log(s.loc[LATEST, "hdd"])]) se = float(np.sqrt(res @ res / (len(y) - 2) * (1 + x0 @ np.linalg.inv(X.T @ X) @ x0))) t = float(st.t.ppf(0.95, len(y) - 2)) pred = float(x0 @ beta) out["loglog_interval"] = {"excess_pct": 100 * (act / np.exp(pred) - 1), "lo_pct": 100 * (act / np.exp(pred + t * se) - 1), "hi_pct": 100 * (act / np.exp(pred - t * se) - 1)} # bootstrap: resample the base winters, re-estimate the rate, add one resampled winter's scatter cb, wb = base["complaints"].to_numpy(float), base["hdd"].to_numpy(float) ex = np.empty(BOOT) for i in range(BOOT): j = rng.integers(0, len(cb), len(cb)) rate = cb[j].sum() / wb[j].sum() e = rng.choice(np.log(cb[j] / (rate * wb[j]))) ex[i] = 100 * (act / (rate * s.loc[LATEST, "hdd"] * np.exp(e)) - 1) lo, med, hi = np.percentile(ex, [5, 50, 95]) out["bootstrap"] = {"lo_pct": lo, "median_pct": med, "hi_pct": hi, "draws": BOOT} # a drift within the base winters would move the benchmark itself tr = st.linregress(np.array(BASE, float), np.log(cb / (cb.sum() / wb.sum() * wb))) m = float(st.t.ppf(0.95, len(BASE) - 2)) * tr.stderr out["base_trend"] = {"pct_a_year": 100 * tr.slope, "lo": 100 * (tr.slope - m), "hi": 100 * (tr.slope + m)} return out def repeats_by_rate(cz: pd.DataFrame) -> dict: """The repeat-building check with high-complaint buildings defined by complaints per apartment rather than by a count: half a complaint or more per apartment in a building of 20 or more apartments, or one or more per apartment in a building of 10 or more. Either way a building needs ten complaints, so every one of them is in the ten-plus file. Apartment counts are today's (HPD's register as retrieved), applied to every season: a current-size check. registry_match_pct gives, season by season, the share of ten-plus buildings the register gives an apartment count.""" t = city.hpd_heat_buildings_ten_plus() g = pd.read_csv(REGISTRY, dtype={"buildingid": str}).rename(columns={"buildingid": "building_id"}) t = t.merge(g[["building_id", "legalclassa"]], on="building_id", how="left") seen = t[t["season"].between(REPEAT_FROM - LOOKBACK, LATEST)] cov = seen.groupby("season")["legalclassa"].apply(lambda x: 100 * float((x > 0).mean())) out = {"registry_match_pct": {label(int(s_)): r(v) for s_, v in cov.items()}} for name, (per, size) in {"half_per_apartment_20plus": (0.5, 20), "one_per_apartment_10plus": (1.0, 10)}.items(): h = t[(t["legalclassa"] >= size) & (t["n"] >= per * t["legalclassa"])] rows = {} for s_ in range(REPEAT_FROM, LATEST + 1): now = h[h["season"] == s_] before = set(h.loc[h["season"].between(s_ - LOOKBACK, s_ - 1), "building_id"]) again = now[now["building_id"].isin(before)] rows[label(s_)] = {"buildings": int(len(now)), "repeat_pct": r(100 * len(again) / len(now)), "repeat_share_of_all": r(100 * again["n"].sum() / cz.loc[s_, "complaints"])} out[name] = rows return out def borough_rates(bo: pd.DataFrame) -> pd.DataFrame: """No-heat complaints per 1,000 renter households outside public housing, by borough and season: the survey's counts (every cycle 1991-2023) interpolated to the January of the season and held at the 2023 count after.""" years = [y for y in hh.YEARS if hh.available(y)] rb = pd.DataFrame({y: hh.renter_households_by_borough(y, public_housing=False) for y in years}).T x = [y + 1 for y in bo.index] den = pd.DataFrame({b: np.interp(x, rb.index.to_numpy(float), rb[b].to_numpy()) for b in BOROUGHS}, index=bo.index) out = 1000 * bo[list(BOROUGHS)] / den out["New York City"] = 1000 * bo[list(BOROUGHS)].sum(axis=1) / den.sum(axis=1) return out # ------------------------------------------------------------------ quotes def quotes(A: Article) -> dict: """Every figure the prose quotes from a report, found in its saved copy.""" Q = A.quotes() rr = "rental-ripoffs-hearing-report-072026.pdf" nh = "nycha_heat_season_2025_26_release.html" nc = "nycha_extreme_cold_release_2026-01-23.html" no = "nycha_heat_season_2026_27_release.html" ct = "nyc_comptroller_turn_up_the_heat_2025.html" hp = "hpd_heat_and_hot_water_page.html" got = { "city_cold": Q.quote("nyc_rental_ripoff_report_2026", rr, "due at least in part to extreme cold weather"), "city_small_group": Q.quote("nyc_rental_ripoff_report_2026", rr, "small group of properties"), "city_owners": Q.quote("nyc_rental_ripoff_report_2026", rr, "a persistent set of noncompliant building owners"), "city_300k": Q.quote("nyc_rental_ripoff_report_2026", rr, "received over 300,000 heat or hot water complaints"), "city_duplicates": Q.quote("nyc_rental_ripoff_report_2026", rr, "will no longer be considered duplicates of each other"), "city_same_apartment": Q.quote("nyc_rental_ripoff_report_2026", rr, "Multiple complaints from the same apartment will still be considered duplicates"), "city_non_anonymous": Q.quote("nyc_rental_ripoff_report_2026", rr, "attempt an inspection for all apartments that file non-anonymous complaints"), "city_primary": Q.quote("nyc_rental_ripoff_report_2026", rr, "all duplicate complaints are closed with the same result"), "city_underreported": Q.quote("nyc_rental_ripoff_report_2026", rr, "housing quality issues are likely underreported"), "nycha_outages": Q.quote("nycha_heat_season_2025_26_release", nh, "2,206 out of 2,240"), "nycha_hours": Q.quote("nycha_heat_season_2025_26_release", nh, "7.8-hour average restoration time"), "nycha_fewer": Q.quote("nycha_heat_season_2026_27_release", no, "six percent total decrease in heat or hot water outages from the 2021-2022 season"), "nycha_contact": Q.quote("nycha_extreme_cold_release_2026_01", nc, "Customer Contact Center at 718-707-7771"), "nycha_mynycha": Q.quote("nycha_extreme_cold_release_2026_01", nc, "MyNYCHA"), "nycha_boilers": Q.quote("nycha_heat_season_2026_27_release", no, "approximately 29 years old"), "nycha_life": Q.quote("nycha_heat_season_2026_27_release", no, "approximately 25 to 30 years"), "hpd_season_count": Q.quote("hpd_heat_and_hot_water_page", hp, "344,437"), "comptroller_1283": Q.quote("nyc_comptroller_turn_up_the_heat_2025", ct, "there are 1,283 buildings"), "comptroller_901": Q.quote("nyc_comptroller_turn_up_the_heat_2025", ct, "approximately 70%, or 901 of those 1,283 buildings"), "comptroller_3pct": Q.quote("nyc_comptroller_turn_up_the_heat_2025", ct, "an average of about 3% of complaints resulted in a violation"), "comptroller_4_5pct": Q.quote("nyc_comptroller_turn_up_the_heat_2025", ct, "increased modestly to 4.5%"), "ll86_night": Q.quote("nyc_local_law_86_2017", "nyc_local_law_86_2017.pdf", "sixty-two"), "ll86_effective": Q.quote("nyc_local_law_86_2017", "nyc_local_law_86_2017.pdf", "take effect on October 1, 2017"), "law_day": Q.quote("hpd_heat_and_hot_water_page", hp, "68 degrees"), "law_night": Q.quote("hpd_heat_and_hot_water_page", hp, "62 degrees"), "law_55": Q.quote("hpd_heat_and_hot_water_page", hp, "55 degrees"), "law_hot_water": Q.quote("hpd_heat_and_hot_water_page", hp, "120 degrees"), "int252_day": Q.quote("nyc_council_int_0252_2026", "nyc_council_int_0252_2026.html", "from 68 degrees to 70 degrees"), "int252_night": Q.quote("nyc_council_int_0252_2026", "nyc_council_int_0252_2026.html", "from 62 degrees to 66 degrees"), "codebook_2021": Q.quote("nychvs_2021_codebook", "2021-nychvs-puf-user-guide-codebook-v2.1.pdf", "October 2019 through May of 2020"), "question_changed": Q.quote("nychvs_2021_initial_findings", "2021-nychvs-selected-initial-findings.pdf", "heating breakdown measure has changed"), "hpd_2021_stabilized": Q.quote("nychvs_2021_initial_findings", "2021-nychvs-selected-initial-findings.pdf", "162,800 ± 16,410 17%"), "hpd_2021_private": Q.quote("nychvs_2021_initial_findings", "2021-nychvs-selected-initial-findings.pdf", "51,800 ± 10,910 5%"), "hpd_2023_stabilized": Q.quote("nychvs_2023_initial_findings", "2023-nychvs-selected-initial-findings.pdf", "188,800 ±12,360 20%"), "hpd_2023_private": Q.quote("nychvs_2023_initial_findings", "2023-nychvs-selected-initial-findings.pdf", "103,500 ±12,090 9%"), "hvs_2026": Q.quote("hpd_nychvs_2026_release", "hpd_nychvs_2026_release.html", "University of Michigan"), "press_80k": Q.quote("gothamist_heat_record_2026_02", "gothamist_heat_record_2026-02-05.html", "Nearly 80,000"), "press_37k": Q.quote("gothamist_heat_record_2026_02", "gothamist_heat_record_2026-02-05.html", "about 37,000"), } return got # ------------------------------------------------------------------ main def main(argv: list[str] | None = None) -> int: """Run the analysis and write every file (no options).""" A = Article(SLUG, __file__, seed=SEED) s = seasons() last, prev = s.loc[LATEST], s.loc[LATEST - 1] base = s.loc[BASE] w30 = weather.heat_seasons().loc[WEATHER_FROM:LATEST] # ---- the season and the cold colder_since = max(y for y in s.index if y < LATEST and s.loc[y, "hdd"] > last["hdd"]) peers = base[(abs(base["hdd"] / last["hdd"] - 1) <= PEER_BAND)] ratios = {"rate": last["ratio_rate"], "loglog": last["ratio_loglog"], "monthly": last["ratio_monthly"], "peers": last["ratio_peers"], "nondup": last["ratio_nondup"], "household": last["ratio_household"]} excess = {k_: 100 * (v - 1) for k_, v in ratios.items()} band = [100 * (np.exp(np.log(last["ratio_rate"]) + z * s.attrs["sd"]) - 1) for z in (-Z90, Z90)] nondup_record_since = min(y for y in s.index if y > 2019 and all( s.loc[q, "nondup"] > s.loc[FIRST:q - 1, "nondup"].max() for q in range(y, LATEST + 1))) A.results["season"] = { "label": label(LATEST), "complaints": int(last["complaints"]), "previous": int(prev["complaints"]), "all_heat_hot_water": int(last["all_heat_hot_water"]), "nondup": int(last["nondup"]), "record_before": int(s.loc[FIRST:LATEST - 1, "complaints"].max()), "record_before_season": label(int(s.loc[FIRST:LATEST - 1, "complaints"].idxmax())), "old_record": int(s.loc[FIRST:LATEST - 2, "complaints"].max()), "old_record_season": label(int(s.loc[FIRST:LATEST - 2, "complaints"].idxmax())), "ratios": {label(y): r(v, 2) for y, v in s["ratio_rate"].items()}, "ratios_household": {label(y): r(v, 2) for y, v in s["ratio_household"].items()}, "ratios_nondup": {label(y): r(v, 2) for y, v in s["ratio_nondup"].items()}, "base": f"{label(BASE[0])} to {label(BASE[-1])}", "base_n": len(BASE), "base_ratio_range": [r(base["ratio_rate"].min(), 2), r(base["ratio_rate"].max(), 2)], "base_ratio_range_nondup": [r(base["ratio_nondup"].min(), 2), r(base["ratio_nondup"].max(), 2)], "base_ratio_range_household": [r(base["ratio_household"].min(), 2), r(base["ratio_household"].max(), 2)], "base_above_fifth": [label(y) for y in base.index if base.loc[y, "ratio_rate"] > 1.2], "above_base_range": [label(y) for y in s.index if s.loc[y, "ratio_rate"] > base["ratio_rate"].max()], "rise_on_previous_pct": r(100 * (last["complaints"] / prev["complaints"] - 1)), "hdd": r(last["hdd"], 0), "colder_since": label(colder_since), "base_rate": r(s.attrs["rate"]), "rate": r(last["complaints"] / last["hdd"]), "expected": int(round(last["expected"])), "noise_sd_log": r(s.attrs["sd"], 3), "peers": {"n": int(len(peers)), "min": int(peers["complaints"].min()), "max": int(peers["complaints"].max()), "seasons": [label(y) for y in peers.index]}, "excess": {"pct": r(excess["rate"]), "lo_pct": r(band[0]), "hi_pct": r(band[1]), "by_check": {k_: r(v) for k_, v in excess.items()}, "elasticity_loglog": r(s.attrs["elasticity"], 2), "note": "lo_pct and hi_pct: 90% band from the base winters' scatter around the benchmark"}, "nondup_record_since": label(nondup_record_since), "nights_10f": int(last["nights_10f"]), "coldest_f": int(last["coldest_f"]), "cold_snap_peers": {label(y): {"nights_10f": int(s.loc[y, "nights_10f"]), "complaints": int(s.loc[y, "complaints"])} for y in s.index if y < LATEST and s.loc[y, "nights_10f"] >= last["nights_10f"]}, } A.results["weather30"] = {"from": label(WEATHER_FROM), "colder_than": int((w30["hdd"] < last["hdd"]).sum()), "of": int(len(w30) - 1), "rank_coldest": int((w30["hdd"] > last["hdd"]).sum()) + 1, "min": r(w30["hdd"].min(), 0), "min_season": label(int(w30["hdd"].idxmin())), "max": r(w30["hdd"].max(), 0), "max_season": label(int(w30["hdd"].idxmax()))} A.claim("season.excess", "up") A.printed("season.complaints", k(last["complaints"])) A.printed("season.expected", k(last["expected"])) A.printed("season.record_before", k(A.results["season"]["record_before"])) A.printed("season.old_record", k(A.results["season"]["old_record"])) A.printed("season.peers.min", k(peers["complaints"].min())) A.printed("season.peers.max", k(peers["complaints"].max())) A.printed("season.excess.pct", pct(excess["rate"])) A.printed("season.excess.lo_pct", pct(band[0])) A.printed("season.excess.hi_pct", pct(band[1])) A.printed("season.excess.by_check.household", pct(excess["household"])) A.printed("season.excess.by_check.nondup", pct(excess["nondup"])) A.printed("season.excess.by_check.loglog", pct(excess["loglog"])) A.printed("season.excess.by_check.monthly", pct(excess["monthly"])) A.printed("season.excess.by_check.peers", pct(excess["peers"])) A.printed("season.base_rate", f"{s.attrs['rate']:.1f}") A.printed("season.hdd", num(last["hdd"])) A.printed("season.nondup", k(last["nondup"])) A.printed("season.all_heat_hot_water", num(last["all_heat_hot_water"])) chk = benchmark_checks(s, A.rng) A.results["season"]["checks"] = { "stations": {k_: {"hdd": r(v["hdd"], 0), "excess_pct": r(v["excess_pct"]), "above_base_range": v["above_base_range"]} for k_, v in chk["stations"].items()}, **{k_: {kk: (r(vv) if isinstance(vv, float) else vv) for kk, vv in chk[k_].items()} for k_ in ("loglog_interval", "bootstrap", "base_trend")}} for stn in ("laguardia", "jfk"): A.printed(f"season.checks.stations.{stn}.excess_pct", pct(chk["stations"][stn]["excess_pct"])) for key in ("loglog_interval", "bootstrap"): A.printed(f"season.checks.{key}.lo_pct", pct(chk[key]["lo_pct"])) A.printed(f"season.checks.{key}.hi_pct", pct(chk[key]["hi_pct"])) over = int(round(last["complaints"] - last["expected"])) assert over < 100_000, "the prose says 'nearly 100,000' more than the benchmark" A.results["season"]["excess_count"] = over A.printed("season.excess_count", f"{int(round(over, -4)):,}") snaps = A.results["season"]["cold_snap_peers"] A.results["season"]["cold_snap_max"] = max(v["complaints"] for v in snaps.values()) A.printed("season.cold_snap_max", k(A.results["season"]["cold_snap_max"])) A.results["season"]["peers"]["mean"] = int(round(peers["complaints"].mean())) A.printed("season.peers.mean", f"{int(round(peers['complaints'].mean(), -4)):,}") dm = daily_model() A.results["season"]["daily_model"] = {**{k_: (r(v, 3 if k_.startswith("sd") else 1) if isinstance(v, float) else v) for k_, v in dm.items() if k_ != "ratios"}, "ratios": {k_: r(v, 2) for k_, v in dm["ratios"].items()}, "note": "lo_pct and hi_pct: 90% reference band from the base winters, each predicted " "by a model fitted without it"} for key in ("excess_pct", "lo_pct", "hi_pct"): A.printed(f"season.daily_model.{key}", pct(dm[key])) A.printed("season.daily_model.predicted", k(dm["predicted"])) # ---- the record month, and the press and City figures it reproduces mo = complaint_months() jan = mo[(mo["month"].dt.year == LATEST + 1) & (mo["month"].dt.month == 1)]["complaints"].sum() every = city.hpd_heat().groupby("month", as_index=False)["problems"].sum() # every month, Feb 2003 on prior_max = every[every["month"] < f"{LATEST}-10-01"].sort_values("problems").iloc[-1].rename({"problems": "complaints"}) c311 = city.heat_311_citywide() jan311 = int(c311[c311["month"] == f"{LATEST + 1}-01-01"]["n"].sum()) allc = city.hpd_heat(all_categories=True) y2025 = int(allc[allc["month"].dt.year == 2025]["problems"].sum()) mo["mon"] = mo["month"].dt.month month_ranks = {} for q, short, full in MONTHS: # the same calendar month in every year of the file, the partial first season included x = every[every["month"].dt.month == q].set_index("month")["problems"] latest = x[x.index.year == (LATEST if q >= 10 else LATEST + 1)].iloc[0] month_ranks[full] = int((x > latest).sum()) + 1 A.results["record_month"] = {"month": f"January {LATEST + 1}", "complaints": int(jan), "previous_max": int(prior_max["complaints"]), "previous_max_month": prior_max["month"].strftime("%B %Y"), "calls_311_all": jan311, "year2025_all": y2025, "ranks": month_ranks} A.printed("record_month.calls_311_all", f"{jan311:,}") A.printed("record_month.complaints", num(jan)) A.printed("record_month.previous_max", num(prior_max["complaints"])) # ---- buildings cz = concentration() rp = repeats(cz) b0, b1 = cz.loc[FIRST], cz.loc[LATEST] sz = building_sizes() A.results["buildings"] = { "first": {"season": label(FIRST), "buildings": int(b0["buildings"]), "top10_share": r(b0["top10_share"]), "top10_share_wo_top": r(b0["top10_share_wo_top"]), "half_buildings": int(b0["half_buildings"]), "half_share": r(b0["half_share_of_buildings"]), "per_building": r(b0["per_building"])}, "latest": {"season": label(LATEST), "buildings": int(b1["buildings"]), "top10_share": r(b1["top10_share"]), "top10_share_wo_top": r(b1["top10_share_wo_top"]), "half_buildings": int(b1["half_buildings"]), "half_share": r(b1["half_share_of_buildings"]), "high_buildings": int(b1["high_buildings"]), "per_building": r(b1["per_building"]), "busiest": int(b1["busiest"])}, "repeats": {label(y): {k_: (r(v) if isinstance(v, float) else int(v)) for k_, v in row.items()} for y, row in rp.iterrows()}, "repeat_pct_before": {"min": r(rp.loc[:LATEST - 1, "repeat_pct"].min()), "max": r(rp.loc[:LATEST - 1, "repeat_pct"].max()), "from": label(REPEAT_FROM), "to": label(LATEST - 1)}, "new_high_buildings": int(rp.loc[LATEST, "buildings"] - rp.loc[LATEST, "repeat"]), "sizes": {k_: (r(v, 2) if isinstance(v, float) else v) for k_, v in sz.items()}, "high": HIGH, "lookback": LOOKBACK, "top_drop": TOP_DROP, } A.printed("buildings.first.top10_share", pct(b0["top10_share"])) A.printed("buildings.latest.top10_share", pct(b1["top10_share"])) A.printed("buildings.latest.high_buildings", f"{int(b1['high_buildings']):,}") A.printed("buildings.first.per_building", f"{b0['per_building']:.1f}") A.printed("buildings.latest.per_building", f"{b1['per_building']:.1f}") A.printed("buildings.first.top10_share_wo_top", pct(b0["top10_share_wo_top"])) A.printed("buildings.latest.top10_share_wo_top", pct(b1["top10_share_wo_top"])) A.printed(f"buildings.repeats.{label(REPEAT_FROM)}.repeat_share_of_all", pct(rp.loc[REPEAT_FROM, "repeat_share_of_all"])) A.printed(f"buildings.repeats.{label(LATEST)}.repeat_share_of_all", pct(rp.loc[LATEST, "repeat_share_of_all"])) A.printed(f"buildings.repeats.{label(LATEST)}.repeat_pct", pct(rp.loc[LATEST, "repeat_pct"])) A.printed("buildings.repeat_pct_before.min", pct(rp.loc[:LATEST - 1, "repeat_pct"].min())) A.printed("buildings.repeat_pct_before.max", pct(rp.loc[:LATEST - 1, "repeat_pct"].max())) A.printed("buildings.new_high_buildings", f"{int(rp.loc[LATEST, 'buildings'] - rp.loc[LATEST, 'repeat']):,}") A.printed("buildings.sizes.high_median_apartments", num(sz["high_median_apartments"])) A.printed("buildings.sizes.other_median_apartments", num(sz["other_median_apartments"])) A.printed("buildings.sizes.high_per_apartment", f"{sz['high_per_apartment']:.2f}") A.printed("buildings.sizes.other_per_apartment", f"{sz['other_per_apartment']:.2f}") by_rate = repeats_by_rate(cz) A.results["buildings"]["by_rate"] = by_rate for s_ in (REPEAT_FROM, LATEST): A.printed(f"buildings.by_rate.half_per_apartment_20plus.{label(s_)}.repeat_share_of_all", pct(by_rate["half_per_apartment_20plus"][label(s_)]["repeat_share_of_all"])) oc = outcomes(s) ob = oc.loc[BASE] A.results["outcomes"] = { "seasons": {label(y): {k_: (r(v) if k_.endswith(("pct", "inspected", "complaints")) and isinstance(v, float) and not k_ in ("first_complaints", "inspected") else int(round(v))) for k_, v in row.items()} for y, row in oc.iterrows()}, "latest": {"violations": int(oc.loc[LATEST, "violations"]), "violation_rate_inspected": r(oc.loc[LATEST, "violation_rate_inspected"]), "share_of_all_complaints": r(oc.loc[LATEST, "share_of_all_complaints"]), "no_access_pct": r(oc.loc[LATEST, "no_access_pct"])}, "base": {"mean_violations": r(ob["violations"].mean(), 0), "violation_rate_inspected": r(100 * ob["violations"].sum() / ob["inspected"].sum()), "share_of_all_complaints": r(100 * ob["violations"].sum() / s.loc[BASE, "complaints"].sum()), "no_access_pct": r((oc.loc[BASE, "no_access_pct"] * oc.loc[BASE, "first_complaints"]).sum() / oc.loc[BASE, "first_complaints"].sum())}, "peak_no_access_pct": r(oc["no_access_pct"].max()), "peak_no_access_season": label(int(oc["no_access_pct"].idxmax())), } A.printed("outcomes.latest.violations", f"{int(round(oc.loc[LATEST, 'violations'], -2)):,}") A.printed("outcomes.base.mean_violations", f"{int(round(ob['violations'].mean(), -2)):,}") A.printed("outcomes.latest.share_of_all_complaints", f"{oc.loc[LATEST, 'share_of_all_complaints']:.1f}%") A.printed("outcomes.base.share_of_all_complaints", f"{100 * ob['violations'].sum() / s.loc[BASE, 'complaints'].sum():.0f}%") det = city.hpd_heat_detail() def detail(season): d = det[det["season"] == season] tot = d["problems"].sum() bw = d["minor_category"].eq("ENTIRE BUILDING") return {"duplicate_pct": r(100 * d.loc[d["duplicate"], "problems"].sum() / tot), "building_wide_duplicate_pct": r(100 * d.loc[d["duplicate"] & bw, "problems"].sum() / tot), "apartment_duplicate_pct": r(100 * d.loc[d["duplicate"] & ~bw, "problems"].sum() / tot), "apartment_only_pct": r(100 * d.loc[d["minor_category"] == "APARTMENT ONLY", "problems"].sum() / tot)} A.results["duplicates"] = {"latest": detail(LATEST), "first_new_scheme": {"season": label(2014), **detail(2014)}, "first": detail(FIRST)} A.printed("duplicates.latest.duplicate_pct", pct(A.results["duplicates"]["latest"]["duplicate_pct"])) A.printed("duplicates.latest.building_wide_duplicate_pct", pct(A.results["duplicates"]["latest"]["building_wide_duplicate_pct"])) bo = borough_complaints() br = borough_rates(bo) A.results["boroughs"] = {b: {"first": int(bo.loc[FIRST, b]), "latest": int(bo.loc[LATEST, b]), "rate_first": r(br.loc[FIRST, b], 1), "rate_latest": r(br.loc[LATEST, b], 1), "rate_change_pct": r(100 * (br.loc[LATEST, b] / br.loc[FIRST, b] - 1))} for b in BOROUGHS} A.results["boroughs"]["New York City"] = {"rate_first": r(br.loc[FIRST, "New York City"], 1), "rate_latest": r(br.loc[LATEST, "New York City"], 1)} A.results["boroughs"]["unit"] = "no-heat complaints per 1,000 renter households outside public housing" for b in ("Bronx", "Brooklyn", "Manhattan", "New York City"): A.printed(f"boroughs.{b}.rate_first", num(br.loc[FIRST, b])) A.printed(f"boroughs.{b}.rate_latest", num(br.loc[LATEST, b])) # ---- public housing in the complaints ph = public_housing_complaints() ph["registry_buildings"] = sz["nycha_buildings"] ph["registry_complaints"] = sz["nycha_complaints"] A.results["public_housing_complaints"] = ph A.printed("public_housing_complaints.complaints", f"{ph['complaints']}") A.printed("public_housing_complaints.registry_complaints", f"{ph['registry_complaints']}") # ---- the survey sv, ch = survey() A.results["survey"] = {"series": {f"{int(x.year)}|{x.group}": {"pct": r(x.pct), "moe": r(x.moe), "n": int(x.n), "winter": x.winter, "margin": x.margin} for x in sv.itertuples()}, "changes": ch, "nonresponse": {str(y): {g: r(hh.nonresponse(y, g)) for g in hh.GROUPS} for y in hh.YEARS if hh.available(y)}} for key in ("renters_1991_2017", "renters_2002_2005", "renters_2008_2011", "stabilized_1991_2017", "private_1991_2017", "public_2002_2017", "public_2008_2017", "renters_2021_2023", "private_2021_2023", "gap_change_1991_2017"): A.claim(f"survey.changes.{key}", "up" if ch[key]["pts"] > 0 else "down") A.claim("survey.changes.public_minus_private.2023", "up") for yy in (2021, 2023): c = hh.called_311(yy) A.results["survey"][f"called_311_{yy}"] = {k_: (r(v) if isinstance(v, float) else v) for k_, v in c.items()} for g in ("Public housing", "Rent stabilized"): cg = hh.called_311(yy, g) A.results["survey"].setdefault("called_311_by_group", {})[f"{yy}|{g}"] = {"pct": r(cg["pct"]), "moe": r(cg["moe"])} pv = {(int(x.year), x.group): x.pct for x in sv.itertuples()} points = (((1991, "All renters"), "renters_1991"), ((2017, "All renters"), "renters_2017"), ((1991, "Public housing"), "public_1991"), ((2002, "Public housing"), "public_2002"), ((2017, "Public housing"), "public_2017"), ((2023, "Public housing"), "public_2023"), ((2023, "Private unregulated"), "private_2023"), ((2021, "All renters"), "renters_2021"), ((2023, "All renters"), "renters_2023"), ((2021, "Public housing"), "public_2021"), ((1991, "Rent stabilized"), "stabilized_1991"), ((2017, "Rent stabilized"), "stabilized_2017"), ((1991, "Private unregulated"), "private_1991"), ((2017, "Private unregulated"), "private_2017"), ((2021, "Rent stabilized"), "stabilized_2021"), ((2023, "Rent stabilized"), "stabilized_2023"), ((2021, "Private unregulated"), "private_2021")) for (yy, g), name in points: A.results["survey"].setdefault("points", {})[name] = r(pv[(yy, g)]) A.printed(f"survey.points.{name}", pct(pv[(yy, g)])) A.printed("survey.called_311_2023.pct", pct(A.results["survey"]["called_311_2023"]["pct"])) A.printed("survey.called_311_2021.pct", pct(A.results["survey"]["called_311_2021"]["pct"])) A.printed("survey.nonresponse.2008.All renters", pct(A.results["survey"]["nonresponse"]["2008"]["All renters"])) # ---- the published tables the survey files reproduce (Census Series IA Table 48, HPD's reports), # in the shape the site's Reproduce box reads (src/data/validation.ts) names = {"sif2021_all": ("All households with a heating breakdown, 2021 (HPD Table 11)", 22_390), "sif2021_stabilized": ("Rent-stabilized households with a heating breakdown, 2021 (HPD)", 16_410), "sif2021_private": ("Private unregulated households with a heating breakdown, 2021 (HPD)", 10_910), "sif2017_no_breakdown": ("Renter households with no heating breakdown, 2017 (HPD Table 18)", None)} rows = [] for x in hh.check_published_reports() + hh.check_published_tables()[::-1]: if x.passed is None: continue if x.name.startswith("table48_"): lab, moe = f"Renter households with a heating breakdown, {x.name[-4:]} (Census Table 48)", None else: lab, moe = names[x.name] rows.append({"statistic": lab, "published": x.target, "published_moe90": moe, "reproduced": r(x.value, 0), "kind": "count", "within_margin": bool(x.passed), "source": x.source}) A.results["reproduction"] = rows A.results["reproduction_source"] = "the Census Bureau's and HPD's published survey tables" A.results["reproduction_note"] = ("Census tables must match to the household (they do, 1999 to 2014); HPD's 2021 " "figures must fall within their published 90% margins. The first eight rows are " "shown; all of them, and the complaint-file checks, are in results.json.") A.results["as_of"] = "2 October 2026" # ---- the checks the page cites (dataset cards), run here and recorded; a failure stops the build checks = city.check_hpd_heat() + weather.check_degree_days() + hh.check_gvf_parameters() \ + hh.check_published_tables() + hh.check_published_reports() A.results["validation"] = [{"name": x.name, "passed": x.passed, "detail": x.detail} for x in checks] failed = [x.name for x in checks if x.passed is False] if failed: raise ValueError(f"validation failed: {failed}") pc = city.hpd_heat().dropna(subset=["season"]).groupby("season")[["complaints", "problems"]].sum() pc.index = pc.index.astype(int) gap_ = 100 * (1 - pc.loc[FIRST:LATEST, "complaints"] / pc.loc[FIRST:LATEST, "problems"]) A.results["season"]["distinct_complaints"] = {"latest_below_pct": r(gap_.loc[LATEST]), "max_below_pct": r(gap_.max()), "max_season": label(int(gap_.idxmax()))} # R12: the overnight standard changed inside the base winters (Local Law 86 of 2017) A.results["law"] = {"overnight_from": "55 degrees when it was below 40 outside", "overnight_to": "62 degrees", "effective": "2017-10-01", "ratio_2017_18": r(s.loc[2017, "ratio_rate"], 2), "ratio_2018_19": r(s.loc[2018, "ratio_rate"], 2)} # ---- quotes A.results["quotes"] = quotes(A) # ---- the current reading and the long view (owner, 2026-10-02) A.latest(f"{LATEST + 1}-05", f"HPD no-heat complaints, heat season {label(LATEST)} (October {LATEST} to May {LATEST + 1})") A.long_view(1991, 2023, "NYC Housing and Vacancy Survey: renter households whose heating broke down, every survey " "1991-2023 (the question changed in 2021, so compared within 1991-2017 and 2021-2023), " "beside thirty winters of Central Park degree days, 1996-97 to 2025-26") # ---- charts xs = [[str(y), label(y)] for y in s.index] ticks = [str(y) for y in s.index if y in (2003, 2008, 2013, 2018, LATEST)] A.chart("seasons_vs_cold", lines( title="Last winter's complaints ran far ahead of its cold", subtitle=f"No-heat complaints to the City each heat season, October to May, {label(FIRST)} to {label(LATEST)}, " f"against the number the season's cold would have brought at the rate of the {len(BASE)} winters " f"before the pandemic", x=xs, x_ticks=ticks, x_axis=SEASON_AXIS, series=[["complaints", "No-heat complaints", "Complaints"], ["expected", f"At the {label(BASE[0])} to {label(BASE[-1])} rate per degree day of cold", "The cold alone"]], cells={"": {"complaints": {str(y): int(v) for y, v in s["complaints"].items()}, "expected": {str(y): int(round(v)) for y, v in s["expected"].items()}}}, format="count", axis="Complaints in the season", note=(f"Complaints to HPD whose problem is no heat, with or without hot water, received October to May; public " f"housing is all but absent from them. The second line is each season's heating degree days at Central " f"Park times {s.attrs['rate']:.1f} complaints per degree day, the rate of the {len(BASE)} winters from " f"{label(BASE[0])} to {label(BASE[-1])} taken together. A computed benchmark, not a forecast."), provenance=PROV + "; NOAA NCEI, Central Park")) A.chart("concentration", lines( title="The busiest tenth of buildings sends a growing share of the complaints", subtitle=f"Share of each heat season's no-heat complaints from the tenth of buildings with the most, " f"{label(FIRST)} to {label(LATEST)}", x=xs, x_ticks=ticks, x_axis=SEASON_AXIS, series=[["top10", "From the busiest tenth of the buildings with any complaint", "Busiest tenth"]], cells={"": {"top10": {str(y): r(v) for y, v in cz["top10_share"].items()}}}, format="pct", axis="Share of the season's complaints", note=(f"Each season, the buildings with at least one no-heat complaint are ranked by their number of complaints " f"({int(b0['buildings']):,} buildings in {label(FIRST)}, {int(b1['buildings']):,} in {label(LATEST)}); " f"the line is the share sent by the top tenth. Buildings are HPD's building records, and none is named."), provenance=PROV)) rx = [[str(y), label(y)] for y in range(REPEAT_FROM, LATEST + 1)] A.chart("repeat_share", lines( title="Repeat buildings send a growing share of the complaints", subtitle=f"Share of each heat season's no-heat complaints from repeat buildings, high-complaint that winter " f"and in one of the ten before, {label(REPEAT_FROM)} to {label(LATEST)}", x=rx, x_ticks=[str(y) for y in (REPEAT_FROM, 2016, 2019, 2022, LATEST)], x_axis=SEASON_AXIS, series=[["count", "Twenty or more complaints in the season", "20+ complaints"], ["per_apt", "Half a complaint or more per apartment (buildings of 20 apartments or more)", "Per apartment"]], cells={"": {"count": {str(y): r(rp.loc[y, "repeat_share_of_all"]) for y in rp.index}, "per_apt": {str(y): by_rate["half_per_apartment_20plus"][label(y)]["repeat_share_of_all"] for y in rp.index}}}, format="pct", axis="Share of the season's no-heat complaints", note=("A repeat building is one over the mark that winter and in at least one of the ten winters before, so " "the series starts in 2013-14, the first winter with ten before it in the record. Apartments come from " "HPD's register of the buildings it oversees. Buildings are HPD's building records, and none is named."), provenance=PROV + "; NYC HPD, Buildings Subject to HPD Jurisdiction")) A.chart("violations", lines( title="More than twice as many no-heat complaints ended in a new violation last winter as before the pandemic", subtitle=f"First complaints of no heat (not duplicates) that HPD closed with a new heat violation, each heat " f"season, {label(FIRST)} to {label(LATEST)}", x=xs, x_ticks=ticks, x_axis=SEASON_AXIS, series=[["violations", "First complaints closed with a new heat violation", "New violations"]], cells={"": {"violations": {str(y): int(v) for y, v in oc["violations"].items()}}}, format="count", axis="First complaints closed with a new violation", note=("Complaints HPD flagged as duplicates take the first complaint's result and are left out. A violation " "is counted when HPD's closing note says a new one was issued; the few closed with one already on file " "are not. How often inspectors got in changed over these years, so the questions give the rate among " "the complaints inspected while heat was required."), provenance=PROV)) years = [y for y in hh.YEARS if hh.available(y)] xsv = [[str(y), str(y)] for y in years if y <= 2017] + [["break", "question changed"]] + \ [[str(y), str(y)] for y in years if y > 2017] cells = {} for g in hh.GROUPS: c = {} for y in years: x = sv[(sv["year"] == y) & (sv["group"] == g)].iloc[0] c[str(y)], c[f"{y}_moe"] = r(x.pct), r(x.moe) cells[g] = c A.chart("survey_groups", lines( title="Heat grew more reliable in private rentals; public housing turned the other way", subtitle="Renter households whose heating broke down for six hours or more that winter, by kind of rental, every " "survey from 1991 to 2023; the question changed in 2021, so the lines break after 2017", x=xsv, x_ticks=["1991", "1999", "2008", "2017", "2023"], series=[["Public housing", "Public housing", "Public housing"], ["Rent stabilized", "Rent stabilized", "Stabilized"], ["All renters", "All renter households", "All renters"], ["Private unregulated", "Private, unregulated (market-rate) rentals", "Unregulated"]], cells={"": cells}, format="pct", axis="Share of renter households", note=("Each survey through 2017 asks about the winter of its interviews; from 2021 it asks about one fixed, " "earlier heat season (October 2019 to May 2020 for the 2021 survey, October 2021 to May 2022 for 2023), " "worded differently, so the line breaks there: compare years within a stretch. Ticks are 90% margins of " "error. Households that did not answer (through 2017) or did not live there that winter (from 2021) are " "left out, so from 2021 the shares run a little above HPD's published ones, which count those households " "as having no breakdown. In 2017 about 62,000 apartments registered as exempt moved from rent stabilized " "to unregulated; HPD cautions against comparing those two groups across that change."), provenance=PROV_HVS, # every survey point has a margin: replicate weights from 2011, the Bureau's variance formula before moe_note="ticks are 90% margins of error")) # FAQ charts order = [m_[0] for m_ in MONTHS] series = [] for y, g in mo.groupby("season"): vals = [int(g.loc[g["mon"] == q, "complaints"].sum()) for q in order] series.append({"key": str(y), "label": label(y), "values": vals, "total": int(sum(vals))}) A.chart("season_months", { "kind": "cycle", "title": f"January {LATEST + 1} was the busiest month in the record", "subtitle": f"No-heat complaints each month of every heat season, {label(FIRST)} to {label(LATEST)}, " f"with {label(LATEST)} drawn over the range of the earlier seasons", "x": [m_[1] for m_ in MONTHS], "x_full": [m_[2] for m_ in MONTHS], "series": series, "highlight": str(LATEST), "band": True, "format": "count", "axis": "No-heat complaints in the month", "note": "Complaints to HPD whose problem is no heat, by the month they were received.", "provenance": PROV}) A.chart("weather30", lines( title=f"Last winter was the coldest since {label(colder_since)}, and colder than most of the last thirty", subtitle=f"Heating degree days at Central Park, October to May, {label(WEATHER_FROM)} to {label(LATEST)}", x=[[str(y), label(y)] for y in w30.index], x_ticks=[str(y) for y in w30.index if y in (1996, 2005, 2015, LATEST)], x_axis=SEASON_AXIS, series=[["hdd", "Heating degree days, October to May", "Degree days"]], cells={"": {"hdd": {str(y): r(v, 0) for y, v in w30["hdd"].items()}}}, format="num", axis="Heating degree days (Fahrenheit, base 65)", note="Each day adds 65 minus its mean temperature in degrees Fahrenheit; a colder season has more. NOAA's monthly " "summary, converted from its Celsius figure.", provenance="NOAA NCEI, Global Summary of the Month, New York Central Park")) ch_ = city.heat_311_channels() ch_ = ch_[(ch_["year"] >= 2014) & (ch_["year"] <= LATEST)] # complete calendar years only known = ch_[ch_["channel"].isin(["PHONE", "ONLINE", "MOBILE"])] tot = known.groupby("year")["n"].sum() share = known.pivot_table(index="year", columns="channel", values="n", aggfunc="sum").div(tot, axis=0) * 100 latest_year = int(share.index.max()) A.chart("channels", lines( title="Since 2021, most heat complaints have come online or by app rather than by phone", subtitle=f"How 311 heat and hot-water complaints arrived, share of each year's complaints, 2014 to {latest_year}", x=[[str(y), str(y)] for y in share.index], x_ticks=[str(y) for y in share.index if y in (2014, 2018, 2021, latest_year)], series=[["PHONE", "By phone", "Phone"], ["ONLINE", "Online", "Online"], ["MOBILE", "Mobile app", "App"]], cells={"": {c: {str(y): r(v) for y, v in share[c].items()} for c in ["PHONE", "ONLINE", "MOBILE"]}}, format="pct", axis="Share of the year's 311 heat complaints", note="Complaints recorded with another or an unknown channel are left out of the shares.", provenance="NYC 311 Service Requests (NYC Open Data)")) A.results["channels"] = {str(y): {c: r(share.loc[y, c]) for c in ["PHONE", "ONLINE", "MOBILE"]} for y in share.index} places = list(BOROUGHS) + ["New York City"] A.chart("boroughs", lines( title="Per renter household, the Bronx's complaints rose most and Brooklyn's barely moved", subtitle=f"No-heat complaints per 1,000 renter households outside public housing, by borough, heat seasons " f"{label(FIRST)} to {label(LATEST)}", x=xs, x_ticks=ticks, x_axis=SEASON_AXIS, series=[[b, b, b] for b in BOROUGHS] + [["New York City", "New York City", "City"]], cells={"": {b: {str(y): int(round(v)) for y, v in br[b].items()} for b in places}}, format="num", axis="Complaints per 1,000 renter households", note=("No-heat complaints to HPD, October to May, by the borough of the building, divided by the borough's renter " "households outside public housing in the Housing and Vacancy Survey, interpolated between surveys and " "held at the 2023 count after."), provenance=PROV + "; " + PROV_HVS)) # ---- the companion data: one tidy table rows = [] for y, x in s.iterrows(): rows.append({"table": "heat_season", "period": label(y), "group": "all", "complaints": int(x["complaints"]), "complaints_not_duplicates": int(x["nondup"]), "all_heat_hot_water": int(x["all_heat_hot_water"]), "heating_degree_days": r(x["hdd"], 1), "renter_households": r(x["renters"], 0), "expected_at_base_rate": r(x["expected"], 0), "ratio_to_expected": r(x["ratio_rate"], 3), "ratio_per_renter_household": r(x["ratio_household"], 3), "ratio_not_duplicates": r(x["ratio_nondup"], 3), "buildings_with_complaints": int(cz.loc[y, "buildings"]), "top10_share_pct": r(cz.loc[y, "top10_share"], 2), "top10_share_without_100_busiest_pct": r(cz.loc[y, "top10_share_wo_top"], 2), "high_buildings": int(cz.loc[y, "high_buildings"]), "high_share_pct": r(cz.loc[y, "high_share"], 2), "repeat_share_of_all_pct": r(rp.loc[y, "repeat_share_of_all"], 2) if y in rp.index else None}) for y, x in oc.iterrows(): rows.append({"table": "complaint_outcomes", "period": label(y), "group": "first complaints", "complaints": int(x["first_complaints"]), "violations": int(x["violations"]), "inspected": int(x["inspected"]), "violation_rate_inspected_pct": r(x["violation_rate_inspected"], 2), "no_access_pct": r(x["no_access_pct"], 2), "restored_pct": r(x["restored_pct"], 2), "daily_model_ratio": r(dm["ratios"][label(y)], 3)}) for y in bo.index: for b in BOROUGHS: rows.append({"table": "heat_season_borough", "period": label(y), "group": b, "complaints": int(bo.loc[y, b]), "complaints_per_1000_renter_households": r(br.loc[y, b], 1)}) for x in sv.itertuples(): rows.append({"table": "survey_breakdown", "period": str(int(x.year)), "group": x.group, "winter": x.winter, "share_pct": r(x.pct, 2), "moe_pct": r(x.moe, 2), "records": int(x.n), "margin": x.margin, "question": x.question}) for y, x in w30.iterrows(): rows.append({"table": "weather", "period": label(y), "group": "Central Park", "heating_degree_days": r(x["hdd"], 1)}) A.write(data=pd.DataFrame(rows)) print(f"{SLUG}: wrote {A.out}") return 0 if __name__ == "__main__": raise SystemExit(main())