"""Who is leaving New York City, who is arriving, and how that changed over ten years — every number and figure in the article, with the working. py src/article_migration.py This file is the article's notebook. The article keeps to what was found; the how and the why live here, next to the code that computes them. ====================================================================================== 1. The question, and the data that answer it ====================================================================================== Every year more people leave New York City for the rest of the country than arrive from it, and every year the gap is filled, or not, by people arriving from abroad and by births. The published numbers are nets. This article opens the nets up: how many people actually arrived and left each year from 2015 to 2024, who they were, where they came from and went to, and where the flows stand in September 2026. Five sources, and the article says which one each number comes from. * **The published mobility table.** Census table B07001 (residence one year ago, for the population one year and over) for New York city from the one-year ACS: 2019 from the Bureau's sequence-format summary file (sequence 18), 2021-2023 from the table-based summary files, 2024 from the Census Reporter mirror. It gives five rows: same house, moved within the county, from another county in the state, from another state, from abroad. The article reproduces each before going anywhere the table does not. What the table cannot give is anyone who *left*. * **The census microdata, two ways.** Arrivals: the 2015-2024 one-year PUMS files for the metro study area (data/interim/acs1/metro_person_.csv.gz), cut to residents of the five boroughs, whose residence-one-year-ago fields say where each person came from. Leavers: the *national* one-year person files, from which extract_migration.py pulled every record with a five-borough migration PUMA a year earlier who now lives elsewhere in the United States (data/interim/migration/leavers_.csv.gz). Both carry the 80 replicate weights, and both are joined to their housing records for household income, tenure and children. Every count, share and margin labelled "census" is a weighted tabulation of those records. The 2020 file carries the Bureau's experimental weights and is flagged wherever it appears. * **The Census Bureau's population estimates** (Vintage 2020 for July 2010 to July 2020, Vintage 2025 for July 2020 to July 2025): births, deaths, net domestic and net international migration for the five counties, summed. Administrative estimates built from tax returns, Medicare enrolment and immigration data, not a survey; they carry the picture to July 2025 and are the numbers the City quotes. * **IRS county-to-county migration data**, filing years 2015-16 to 2022-23: returns, exemptions (roughly people) and adjusted gross income for every county pair, so the income of the households that filed from the city one year and from somewhere else the next, and where they went. Filers only, and the latest year is 2022-23. * **Two series that run into 2026.** The Department of Homeless Services' daily shelter census to September 2026, the administrative shadow of the 2022-2024 arrivals; and the Current Population Survey basic monthly files, June 2024 to August 2026, for the foreign-born share of the city's residents (data/interim/cps_nativity_nyc.csv.gz, built by extract_cps_nativity.py). ====================================================================================== 2. Universe and definitions ====================================================================================== * **Residents**: people living in the five boroughs, one year old and over, in the survey year, including group quarters (dormitories, nursing homes, shelters), which is the universe of table B07001. Households only where a household variable (income, tenure, children) is used. * **Flows**, from the residence-one-year-ago items. A resident either stayed in the same house, moved within the city (previous migration PUMA in the five boroughs), or arrived: from the metro suburbs (a migration PUMA in the Long Island, Lower Hudson Valley, northern New Jersey or southwestern Connecticut part of the study area, using the PUMA-to-migration-PUMA composition files), from elsewhere in New York State, from elsewhere in New Jersey or Connecticut, from another state, or from abroad (outside the fifty states and DC, so including Puerto Rico, as the published table does). A leaver lived in the five boroughs a year ago and now lives elsewhere in the United States; destinations are the same regions plus Florida, Pennsylvania, California, Texas, the rest of the Northeast, South, Midwest and West. People who left the country are in no survey. * **Net domestic migration** (census): arrivals from elsewhere in the United States minus leavers to elsewhere in the United States, computed within one year's file with its replicate weights. The Bureau's population-estimates figure of the same name is a different measurement (tax and Medicare records, a July-to-July year). * **Out-migration rate**: leavers divided by everyone who lived in the city a year ago and is still in a US survey (stayers, within-city movers and leavers). It omits people who died or left the country in the year. * **Household income**: the household's total income, in 2024 dollars (the Bureau's within-year factor, then the New York-metro CPI). For a mover it is the income of the household they live in *now*, after the move. * **Degree**: bachelor's or higher, among people 25 and over. **Employed**: 16 and over. **Race/ethnicity**: Hispanic of any race; otherwise White, Black, Asian, Other. **With children**: someone under 18 lives in the household. **Owner**: the household owns its home, with or without a mortgage. ====================================================================================== 3. Margins ====================================================================================== PUMS ships 80 replicate weights per person. For any statistic θ computed with the full weight, the same statistic is recomputed with each replicate weight r, and Var(θ) = (4 / 80) · Σ_r (θ_r − θ)² Margins are 90% (1.645 SE), the Bureau's convention. Medians are weighted medians without a margin. Differences and nets are computed replicate by replicate, so the margin on a net flow accounts for both sides. The population estimates, the IRS data and the shelter census carry no sampling margin. The CPS public file has no replicate weights; its margins are approximate, from the spread of the monthly estimates. Outputs (all regenerable; nothing is edited by hand): output/articles/who_is_leaving_and_who_is_arriving/results.json every number quoted in the prose output/articles/who_is_leaving_and_who_is_arriving/pums_direct.csv census-side estimates with margins output/articles/who_is_leaving_and_who_is_arriving/reproduction.csv the published-table reproduction output/articles/who_is_leaving_and_who_is_arriving/charts.json cubes behind the interactive charts output/articles/who_is_leaving_and_who_is_arriving/fig*.svg static twins of the same figures output/articles/who_is_leaving_and_who_is_arriving/data.csv companion dataset (re-randomised ids) site/public/articles/who-is-leaving-and-who-is-arriving/ the same, copied for the website """ from __future__ import annotations import argparse import json import shutil from pathlib import Path import numpy as np import pandas as pd from _paths import P, INTERIM, OUTPUT, RAW, REFERENCE import article_migration_charts as CH from article_commute import weighted_median SLUG = "who-is-leaving-and-who-is-arriving" OUT = OUTPUT / "articles" / "who_is_leaving_and_who_is_arriving" SITE_OUT = P.root / "site" / "public" / "articles" / SLUG MIG_DIR = INTERIM / "migration" Z90 = 1.645 MIN_RECORDS = 30 YEARS = list(range(2015, 2025)) FOCUS = (2015, 2019, 2024) BOROUGHS = ["Manhattan", "Brooklyn", "Queens", "Bronx", "Staten Island"] NYC_MIGPUMA = {"2010": {3700: "Bronx", 3800: "Manhattan", 3900: "Staten Island", 4000: "Brooklyn", 4100: "Queens"}, "2020": {4100: "Manhattan", 4200: "Bronx", 4300: "Brooklyn", 4400: "Queens", 4500: "Staten Island"}} STATE_NAMES = {1: "Alabama", 2: "Alaska", 4: "Arizona", 5: "Arkansas", 6: "California", 8: "Colorado", 9: "Connecticut", 10: "Delaware", 11: "District of Columbia", 12: "Florida", 13: "Georgia", 15: "Hawaii", 16: "Idaho", 17: "Illinois", 18: "Indiana", 19: "Iowa", 20: "Kansas", 21: "Kentucky", 22: "Louisiana", 23: "Maine", 24: "Maryland", 25: "Massachusetts", 26: "Michigan", 27: "Minnesota", 28: "Mississippi", 29: "Missouri", 30: "Montana", 31: "Nebraska", 32: "Nevada", 33: "New Hampshire", 34: "New Jersey", 35: "New Mexico", 36: "New York", 37: "North Carolina", 38: "North Dakota", 39: "Ohio", 40: "Oklahoma", 41: "Oregon", 42: "Pennsylvania", 44: "Rhode Island", 45: "South Carolina", 46: "South Dakota", 47: "Tennessee", 48: "Texas", 49: "Utah", 50: "Vermont", 51: "Virginia", 53: "Washington", 54: "West Virginia", 55: "Wisconsin", 56: "Wyoming"} NORTHEAST = {9, 23, 25, 33, 34, 36, 42, 44, 50} MIDWEST = {17, 18, 19, 20, 26, 27, 29, 31, 38, 39, 46, 55} SOUTH = {1, 5, 10, 11, 12, 13, 21, 22, 24, 28, 37, 40, 45, 47, 48, 51, 54} WEST = {2, 4, 6, 8, 15, 16, 30, 32, 35, 41, 49, 53, 56} SUBURB_REGIONS = {"long_island": "Long Island", "lower_hudson": "Lower Hudson Valley", "north_jersey": "Northern New Jersey", "southwest_ct": "Southwestern Connecticut"} ORIGINS = [("abroad", "From abroad"), ("other_state", "From another state"), ("suburbs", "From the metro suburbs"), ("rest_ny", "From elsewhere in New York State"), ("rest_nj_ct", "From elsewhere in New Jersey or Connecticut")] DESTINATIONS = [("suburbs", "To the metro suburbs"), ("rest_ny", "To elsewhere in New York State"), ("rest_nj_ct", "To elsewhere in New Jersey or Connecticut"), ("florida", "To Florida"), ("pennsylvania", "To Pennsylvania"), ("other_northeast", "To the rest of the Northeast"), ("california", "To California"), ("texas", "To Texas"), ("other_south", "To the rest of the South"), ("midwest", "To the Midwest"), ("other_west", "To the rest of the West")] SUBURB_DESTS = [("long_island", "Long Island"), ("lower_hudson", "Lower Hudson Valley"), ("north_jersey", "Northern New Jersey"), ("southwest_ct", "Southwestern Connecticut")] STATUSES = [("stayed", "Did not move"), ("within", "Moved within the city"), ("arrived_domestic", "Arrived from elsewhere in the US"), ("arrived_abroad", "Arrived from abroad"), ("left", "Left for elsewhere in the US")] AGE_BANDS = [("1_17", "Under 18", 1, 18), ("18_24", "18–24", 18, 25), ("25_34", "25–34", 25, 35), ("35_44", "35–44", 35, 45), ("45_64", "45–64", 45, 65), ("65_plus", "65 and over", 65, 200)] INCOME_BANDS = [("lt25k", "Under $25k", -np.inf, 25_000), ("25k_50k", "$25k–50k", 25_000, 50_000), ("50k_100k", "$50k–100k", 50_000, 100_000), ("100k_150k", "$100k–150k", 100_000, 150_000), ("150k_200k", "$150k–200k", 150_000, 200_000), ("200k_300k", "$200k–300k", 200_000, 300_000), ("300k_plus", "$300k and over", 300_000, np.inf)] RACES = ["White", "Black", "Hispanic", "Asian", "Other"] FACTORS = {"age": ("Age", [(k, l) for k, l, *_ in AGE_BANDS]), "income": ("Household income now, after any move (2024 dollars)", [(k, l) for k, l, *_ in INCOME_BANDS]), "degree": ("Degree (25 and over)", [("ba_plus", "Bachelor's or higher"), ("no_degree", "No bachelor's")]), "race": ("Race or ethnicity", [(r, r) for r in RACES]), "children": ("Children in the household now", [("with_children", "Lives with children"), ("no_children", "No children in the household")]), "nativity": ("Nativity", [("us_born", "US-born"), ("foreign_born", "Foreign-born")]), "tenure": ("Tenure of the home now", [("owner", "Lives in an owned home"), ("renter", "Lives in a rented home")]), "borough": ("Borough", [(b, b if b != "Bronx" else "The Bronx") for b in BOROUGHS])} WORLD_REGIONS = [("caribbean_central", "Caribbean, Mexico and Central America"), ("south_america", "South America"), ("asia", "Asia"), ("europe", "Europe"), ("africa", "Africa"), ("other_world", "Canada, Oceania and elsewhere")] TRAITS = [("under_18", "Under 18"), ("age_18_34", "Aged 18 to 34"), ("age_65_plus", "65 and over"), ("ba_plus", "Bachelor's degree or higher (25+)"), ("employed", "Employed (16+)"), ("income_150k_plus", "Household income now $150k or more"), ("income_lt50k", "Household income now under $50k"), ("foreign_born", "Foreign-born"), ("noncitizen", "Not a US citizen"), ("white", "White"), ("black", "Black"), ("hispanic", "Hispanic"), ("asian", "Asian"), ("with_children", "Lives with children now"), ("owner", "Lives in an owned home now"), ("student", "Enrolled in college or school (18+)")] INK = "#6f6e69"; GRID = "#d5d4cd"; BLUE = "#2a78d6"; ORANGE = "#eb6834"; GREEN = "#4a9c6d"; TICK = "#1a1a1a"; FIG_W = 6.4 SF_2019 = RAW / "acs_sf" / "NewYork_All_Geographies_2019_1yr" TABLES_2024 = REFERENCE / "tables" / "acs_mobility_nyc_us_2024_1yr.json" NYC_COUNTIES = {"005": "Bronx", "047": "Brooklyn", "061": "Manhattan", "081": "Queens", "085": "Staten Island"} SUBURB_COUNTIES_IRS = {"36059": "long_island", "36103": "long_island", "36119": "lower_hudson", "36087": "lower_hudson", "36079": "lower_hudson", "36071": "lower_hudson", "34003": "north_jersey", "34017": "north_jersey", "34013": "north_jersey", "34039": "north_jersey", "34031": "north_jersey", "34023": "north_jersey", "34025": "north_jersey", "34027": "north_jersey", "34035": "north_jersey", "34037": "north_jersey", "34019": "north_jersey", "34029": "north_jersey", "09001": "southwest_ct"} # ------------------------------------------------------------------ helpers def cpi_deflators() -> dict[int, float]: d = pd.read_csv(REFERENCE / "cpi_annual.csv"); d = d[d["name"] == "cpi_ny_all_items"] return dict(zip(d["year"].astype(int), d["deflator_to_2024"].astype(float))) def vintage(year: int) -> str: return "2010" if year <= 2021 else "2020" def _geo(v: str) -> pd.DataFrame: f = "geography_pumas_2010.csv" if v == "2010" else "geography_pumas.csv" g = pd.read_csv(REFERENCE / f, dtype={"state_fips": str, "puma": str}) g["puma_geoid"] = g["state_fips"].str.zfill(2) + g["puma"].str.zfill(5) g["in_nyc"] = pd.to_numeric(g["in_nyc"]).astype(int) return g[["puma_geoid", "in_nyc", "borough", "region"]] def migpuma_regions(v: str) -> dict[tuple[int, int], str]: """(state, migration PUMA) -> study-area region, from the composition files.""" f = REFERENCE / ("puma_migpuma_2010.xls" if v == "2010" else "puma_migpuma_2020.xls") x = pd.read_excel(f, dtype=str, header=None) x = x[x.iloc[:, 0].str.fullmatch(r"\d{2}", na=False)].iloc[:, :4] x.columns = ["st", "puma", "migst", "migpuma"] x["puma_geoid"] = x["st"].str.zfill(2) + x["puma"].str.zfill(5) g = _geo(v); m = x.merge(g, on="puma_geoid", how="inner") m["key"] = list(zip(m["migst"].astype(int), m["migpuma"].astype(int))) out = {} for key, grp in m.groupby("key"): # A migration PUMA that mixes city and suburb, or suburb and outside, takes its majority region. reg = grp["region"].mode().iloc[0] out[key] = reg return out def race_label(rac1p, hisp) -> np.ndarray: r = pd.to_numeric(rac1p, errors="coerce").to_numpy(); h = pd.to_numeric(hisp, errors="coerce").to_numpy() return np.where(h > 1, "Hispanic", np.where(r == 1, "White", np.where(r == 2, "Black", np.where(r == 6, "Asian", "Other")))) def age_band(a) -> np.ndarray: v = np.asarray(a, float) return np.select([v < 18, v < 25, v < 35, v < 45, v < 65], ["1_17", "18_24", "25_34", "35_44", "45_64"], "65_plus") def income_band(v) -> np.ndarray: v = np.asarray(v, float) out = np.select([v < 25_000, v < 50_000, v < 100_000, v < 150_000, v < 200_000, v < 300_000], ["lt25k", "25k_50k", "50k_100k", "100k_150k", "150k_200k", "200k_300k"], "300k_plus") return np.where(np.isnan(v), "na", out) def world_region_last_year(migsp) -> np.ndarray: """Where an arrival from abroad lived a year ago, from the migration state/country code: Puerto Rico and the island areas, then the same regions as place of birth.""" p = np.asarray(migsp, float) return np.select([p < 100, (p >= 100) & (p < 200), (p >= 200) & (p < 300), (p >= 300) & (p < 360), (p >= 360) & (p < 400), (p >= 400) & (p < 500)], ["puerto_rico", "europe", "asia", "caribbean_central", "south_america", "africa"], "other_world") def migsp_names() -> dict[int, str]: d = pd.read_csv(REFERENCE / "PUMS_Data_Dictionary_2020-2024.csv", header=None, dtype=str, names=list("abcdefg"), on_bad_lines="skip") d = d[(d["a"] == "VAL") & (d["b"] == "MIGSP") & d["e"].str.fullmatch(r"\d+")] return {int(k): v for k, v in zip(d["e"], d["g"])} def world_region(pobp) -> np.ndarray: p = np.asarray(pobp, float) return np.select([p < 100, (p >= 100) & (p < 200), (p >= 200) & (p < 300), (p >= 300) & (p < 360), (p >= 360) & (p < 400), (p >= 400) & (p < 500)], ["us", "europe", "asia", "caribbean_central", "south_america", "africa"], "other_world") def state_group(st) -> np.ndarray: s = np.asarray(st, float) return np.select([s == 36, np.isin(s, [34, 9]), s == 12, s == 42, s == 6, s == 48, np.isin(s, list(NORTHEAST)), np.isin(s, list(SOUTH)), np.isin(s, list(MIDWEST)), np.isin(s, list(WEST))], ["rest_ny", "rest_nj_ct", "florida", "pennsylvania", "california", "texas", "other_northeast", "other_south", "midwest", "other_west"], "other_west") def pobp_names() -> dict[int, str]: d = pd.read_csv(REFERENCE / "PUMS_Data_Dictionary_2020-2024.csv", header=None, dtype=str, names=list("abcdefg"), on_bad_lines="skip") d = d[(d["a"] == "VAL") & (d["b"] == "POBP")] return {int(k): v for k, v in zip(d["e"], d["g"])} # ------------------------------------------------------------------ loading PERSON_COLS = ["SERIALNO", "PUMA", "STATE", "PWGTP", "AGEP", "SEX", "RAC1P", "HISP", "NATIVITY", "CIT", "POBP", "YOEP", "SCHL", "SCH", "ESR", "COW", "PERNP", "PINCP", "ADJINC", "MIG", "MIGSP", "MIGPUMA", "RELSHIPP", "JWTRNS", "POWSP", "POWPUMA"] HH_COLS = ["SERIALNO", "HINCP", "TEN", "NP", "NOC", "HHT", "MV", "TYPEHUGQ", "R18", "R65", "WGTP"] def _derive(h: pd.DataFrame, year: int, defl: float) -> pd.DataFrame: h["age_band"] = age_band(h["AGEP"]) h["race"] = race_label(h["RAC1P"], h["HISP"]) h["foreign_born"] = pd.to_numeric(h["NATIVITY"], errors="coerce") == 2 cit = pd.to_numeric(h["CIT"], errors="coerce") h["noncitizen"] = cit == 5; h["naturalized"] = cit == 4 h["adult25"] = h["AGEP"] >= 25; h["adult16"] = h["AGEP"] >= 16 h["ba_plus"] = (pd.to_numeric(h["SCHL"], errors="coerce") >= 21) & h["adult25"] h["employed"] = pd.to_numeric(h["ESR"], errors="coerce").isin([1, 2, 4, 5]) & h["adult16"] h["student"] = pd.to_numeric(h["SCH"], errors="coerce").isin([2, 3]) & (h["AGEP"] >= 18) adj = pd.to_numeric(h["ADJINC"], errors="coerce") / 1e6 hinc = pd.to_numeric(h["HINCP"], errors="coerce") gq = pd.to_numeric(h["TYPEHUGQ"], errors="coerce") != 1 h["hh_income"] = np.where(gq | hinc.isna(), np.nan, hinc * adj * defl) h["income_band"] = income_band(h["hh_income"]) h["household"] = ~gq & hinc.notna() ten = pd.to_numeric(h["TEN"], errors="coerce") h["owner"] = ten.isin([1, 2]); h["renter"] = ten.isin([3, 4]) h["with_children"] = pd.to_numeric(h["R18"], errors="coerce") == 1 h["world_region"] = world_region(h["POBP"]) h["from_region"] = np.where(h["status"] == "arrived_abroad", world_region_last_year(pd.to_numeric(h["MIGSP"], errors="coerce")), "") h["wfh"] = pd.to_numeric(h["JWTRNS"], errors="coerce") == 11 h["year"] = year return h def load_year(year: int) -> tuple[pd.DataFrame, np.ndarray]: """One frame per year: NYC residents (stayed, moved within, arrived) stacked with the leavers found in the national file, with replicate weights in the same order.""" v = vintage(year); geo = _geo(v); codes = NYC_MIGPUMA[v]; code_of = {b: c for c, b in codes.items()} regions = migpuma_regions(v); defl = cpi_deflators().get(year, 1.0) # Residents. p = pd.read_csv(INTERIM / "acs1" / f"metro_person_{year}.csv.gz", usecols=PERSON_COLS, dtype={"STATE": str, "PUMA": str, "SERIALNO": str, "POWPUMA": str}) rw = np.load(INTERIM / "acs1" / f"metro_person_{year}_repwts.npy", mmap_mode="r") assert len(rw) == len(p) p["puma_geoid"] = p["STATE"].str.zfill(2) + p["PUMA"].str.zfill(5) p = p.merge(geo, on="puma_geoid", how="left") keep = ((p["in_nyc"] == 1) & p["MIG"].notna()).to_numpy() r = p[keep].reset_index(drop=True); rw_r = np.asarray(rw[keep], dtype=np.float32) hh = pd.read_csv(INTERIM / "acs1" / f"metro_housing_{year}.csv.gz", usecols=lambda c: c in HH_COLS, dtype={"SERIALNO": str}) r = r.merge(hh, on="SERIALNO", how="left") mig = pd.to_numeric(r["MIG"], errors="coerce"); sp = pd.to_numeric(r["MIGSP"], errors="coerce"); mp = pd.to_numeric(r["MIGPUMA"], errors="coerce") own_code = r["borough"].map(code_of) origin_region = pd.Series([regions.get((int(s), int(m))) if s == s and m == m else None for s, m in zip(sp, mp)], index=r.index) us_state = sp.between(1, 56) r["status"] = np.select([mig == 1, (mig == 3) & (sp == 36) & mp.isin(list(codes)), (mig == 2) | ((mig == 3) & ~us_state)], ["stayed", "within", "arrived_abroad"], "arrived_domestic") r["origin"] = np.select([r["status"] == "arrived_abroad", r["status"] != "arrived_domestic", origin_region.isin(list(SUBURB_REGIONS)), sp == 36, sp.isin([34, 9])], ["abroad", "", "suburbs", "rest_ny", "rest_nj_ct"], "other_state") r["origin_region"] = np.where(r["origin"] == "suburbs", origin_region.fillna(""), "") r["same_county"] = (mig == 3) & (sp == 36) & (mp == own_code) r["diff_county_same_state"] = (mig == 3) & (sp == 36) & (mp != own_code) r["destination"] = ""; r["dest_region"] = ""; r["origin_borough"] = np.where(r["status"].isin(["stayed", "within"]), r["borough"], "") r.loc[r["status"] == "within", "origin_borough"] = mp[r["status"] == "within"].map(codes) r["now_borough"] = r["borough"]; r["works_in_nyc"] = np.nan r = _derive(r, year, defl) # Leavers. l = pd.read_csv(MIG_DIR / f"leavers_{year}.csv.gz", dtype={"STATE": str, "PUMA": str, "SERIALNO": str, "POWPUMA": str}) rw_l = np.load(MIG_DIR / f"leavers_{year}_repwts.npy") assert len(rw_l) == len(l) keep = (~l["now_nyc"]).to_numpy() l = l[keep].reset_index(drop=True); rw_l = np.asarray(rw_l[keep], dtype=np.float32) lhh = pd.read_csv(MIG_DIR / f"leavers_{year}_hh.csv.gz", usecols=lambda c: c in HH_COLS, dtype={"SERIALNO": str}) l = l.merge(lhh, on="SERIALNO", how="left") l["puma_geoid"] = l["STATE"].str.zfill(2) + l["PUMA"].str.zfill(5) l = l.merge(geo, on="puma_geoid", how="left") st = pd.to_numeric(l["STATE"], errors="coerce") l["status"] = "left"; l["origin"] = ""; l["origin_region"] = ""; l["same_county"] = False; l["diff_county_same_state"] = False suburb = l["region"].isin(list(SUBURB_REGIONS)) l["destination"] = np.where(suburb, "suburbs", state_group(st)) l["dest_region"] = np.where(suburb, l["region"].fillna(""), "") l["dest_state"] = st.map(STATE_NAMES) l["now_borough"] = ""; l["borough"] = l["origin_borough"] pow_state = pd.to_numeric(l["POWSP"], errors="coerce"); pow_puma = pd.to_numeric(l["POWPUMA"], errors="coerce") l["works_in_nyc"] = np.where(pd.to_numeric(l["ESR"], errors="coerce").isin([1, 2, 4, 5]) & pow_state.notna(), (pow_state == 36) & pow_puma.isin(list(codes)), np.nan) l = _derive(l, year, defl) cols = sorted(set(r.columns) & set(l.columns)) h = pd.concat([r[cols], l[cols + ["dest_state"]]], ignore_index=True) h["resident"] = h["status"] != "left" return h, np.concatenate([rw_r, rw_l]) 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["PWGTP"].to_numpy(dtype=float) @staticmethod def _moe(full, reps): return float(Z90 * np.sqrt(4 / 80 * ((reps - full) ** 2).sum())) def count_reps(self, mask): mask = np.asarray(mask); return float(self.w[mask].sum()), self.rw[mask].sum(0).astype(float) def count(self, mask) -> tuple[float, float]: full, reps = self.count_reps(mask); return full, self._moe(full, reps) def diff(self, mask_a, mask_b) -> tuple[float, float]: """a minus b, replicate by replicate.""" fa, ra = self.count_reps(mask_a); fb, rb = self.count_reps(mask_b) return fa - fb, self._moe(fa - fb, ra - rb) 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) if w.sum() == 0: return float("nan"), float("nan") full = (w * x).sum() / w.sum(); reps = (rw * x[:, None]).sum(0) / rw.sum(0) return float(full), self._moe(full, reps) def median(self, mask, col: str) -> float: mask = np.asarray(mask); v = self.h.loc[mask, col].to_numpy(dtype=float); ok = ~np.isnan(v) if ok.sum() < MIN_RECORDS: return float("nan") return weighted_median(v[ok], self.w[mask][ok]) # ------------------------------------------------------------------ published tables def published_tables() -> dict: """B07001 as published for New York city: the five residence-one-year-ago rows.""" cells = {"total": 1, "same_house": 17, "same_county": 33, "diff_county_same_state": 49, "different_state": 65, "abroad": 81} out = {"table_note": "B07001: population 1 year and over by residence 1 year ago; cells 1, 17, 33, 49, 65 and 81 are the totals of each row."} g = pd.read_csv(SF_2019 / "g20191ny.csv", header=None, dtype=str, encoding="latin-1"); logrec = dict(zip(g[48], g[4])) e = pd.read_csv(SF_2019 / "e20191ny0018000.txt", header=None, dtype=str, low_memory=False) m = pd.read_csv(SF_2019 / "m20191ny0018000.txt", header=None, dtype=str, low_memory=False) lr = logrec["16000US3651000"]; er = e[e[5] == lr].iloc[0]; mr = m[m[5] == lr].iloc[0]; start = 7 - 1 out["2019"] = {k: {"estimate": float(er[start + c - 1]), "moe90": float(mr[start + c - 1])} for k, c in cells.items()} for y in (2021, 2022, 2023): d = pd.read_csv(REFERENCE / "tables" / f"acsdt1y{y}-b07001.dat", sep="|", dtype=str).set_index("GEO_ID") r = d.loc["1600000US3651000"] out[str(y)] = {k: {"estimate": float(r[f"B07001_E{c:03d}"]), "moe90": float(r[f"B07001_M{c:03d}"])} for k, c in cells.items()} j = json.loads(TABLES_2024.read_text(encoding="utf-8")); b = j["data"]["16000US3651000"]["B07001"] out["2024"] = {k: {"estimate": b["estimate"][f"B07001{c:03d}"], "moe90": b["error"][f"B07001{c:03d}"]} for k, c in cells.items()} bu = j["data"]["01000US"]["B07001"] out["2024_us"] = {k: {"estimate": bu["estimate"][f"B07001{c:03d}"], "moe90": bu["error"][f"B07001{c:03d}"]} for k, c in cells.items()} return out # ------------------------------------------------------------------ census estimates def year_estimates(cs: Census) -> tuple[dict, list[dict]]: h = cs.h; year = int(h["year"].iloc[0]); rows = [] status = h["status"].to_numpy(); res = h["resident"].to_numpy(); origin = h["origin"].to_numpy(); dest = h["destination"].to_numpy() oreg = h["origin_region"].to_numpy(); dreg = h["dest_region"].to_numpy(); boro = h["borough"].to_numpy() stayed = status == "stayed"; within = status == "within"; arr_dom = status == "arrived_domestic"; arr_abr = status == "arrived_abroad"; left = status == "left" arrived = arr_dom | arr_abr; lived_here = stayed | within | left def put(key, label, value, moe, n, kind): rows.append({"year": year, "key": key, "label": label, "value": value, "moe90": moe, "n_records": int(n), "kind": kind}) out: dict = {"year": year, "n_records": {"residents": int(res.sum()), "leavers": int(left.sum())}, "experimental_weights": year == 2020} pop, popm = cs.count(res); out["population"] = {"estimate": pop, "moe90": popm} put("population", "NYC residents 1 and over", pop, popm, res.sum(), "count") out["flows"] = {} for k, label in STATUSES: mk = status == k; c, cm = cs.count(mk) out["flows"][k] = {"label": label, "count": c, "moe90": cm, "n_records": int(mk.sum()), "share_of_population": c / pop} put(f"flow_{k}", label, c, cm, mk.sum(), "count") c, cm = cs.count(arrived); out["flows"]["arrived_all"] = {"label": "Arrived from anywhere", "count": c, "moe90": cm, "n_records": int(arrived.sum())} net, netm = cs.diff(arr_dom, left); out["net_domestic"] = {"estimate": net, "moe90": netm} put("net_domestic", "Net domestic migration (arrivals from elsewhere in the US minus leavers)", net, netm, (arr_dom | left).sum(), "count") net_all, net_allm = cs.diff(arrived, left); out["net_all_observed"] = {"estimate": net_all, "moe90": net_allm} rate, ratem = cs.share(lived_here, left); out["out_rate"] = {"estimate": rate, "moe90": ratem} put("out_rate", "Share of last year's residents who left for elsewhere in the US", rate, ratem, lived_here.sum(), "share") out["in_rate"] = {"domestic": cs.share(res, arr_dom)[0], "abroad": cs.share(res, arr_abr)[0], "any": cs.share(res, arrived)[0]} out["within_rate"] = cs.share(res, within)[0] # Origins and destinations. out["origins"] = {} for k, label in ORIGINS: mk = arrived & (origin == k); c, cm = cs.count(mk) out["origins"][k] = {"label": label, "count": c, "moe90": cm, "n_records": int(mk.sum()), "share_of_arrivals": c / out["flows"]["arrived_all"]["count"]} put(f"origin_{k}", label, c, cm, mk.sum(), "count") out["origins_suburbs"] = {k: dict(zip(("count", "moe90"), cs.count(arrived & (oreg == k))), label=l) for k, l in SUBURB_DESTS} out["destinations"] = {} for k, label in DESTINATIONS: mk = left & (dest == k); c, cm = cs.count(mk) out["destinations"][k] = {"label": label, "count": c, "moe90": cm, "n_records": int(mk.sum()), "share_of_leavers": c / out["flows"]["left"]["count"]} put(f"destination_{k}", label, c, cm, mk.sum(), "count") out["destinations_suburbs"] = {k: dict(zip(("count", "moe90"), cs.count(left & (dreg == k))), label=l) for k, l in SUBURB_DESTS} st = h["dest_state"].to_numpy(dtype=object) top = pd.Series(cs.w[left]).groupby(st[left]).sum().sort_values(ascending=False) out["destination_states"] = [{"state": s, "count": float(v), "share_of_leavers": float(v / out["flows"]["left"]["count"])} for s, v in top.head(12).items()] # Exchanges: net with each region. out["exchange"] = {} for k, label in [("suburbs", "The metro suburbs"), ("rest_ny", "Elsewhere in New York State"), ("rest_nj_ct", "Elsewhere in New Jersey and Connecticut"), ("other_state", "Other states")]: a = arrived & (origin == k); l_ = left & ((dest == k) if k != "other_state" else ~np.isin(dest, ["suburbs", "rest_ny", "rest_nj_ct", ""])) n_, nm = cs.diff(a, l_) out["exchange"][k] = {"label": label, "arrived": cs.count(a)[0], "left": cs.count(l_)[0], "net": n_, "net_moe90": nm} # Boroughs: arrivals by current borough, leavers by origin borough. out["boroughs"] = {} for b in BOROUGHS: mb = boro == b; bp = cs.count(mb & res)[0] n_, nm = cs.diff(mb & arr_dom, mb & left); r_, rm = cs.share(mb & lived_here, left) out["boroughs"][b] = {"population": bp, "arrived_domestic": cs.count(mb & arr_dom)[0], "arrived_abroad": cs.count(mb & arr_abr)[0], "left": cs.count(mb & left)[0], "within": cs.count(mb & within)[0], "net_domestic": n_, "net_domestic_moe90": nm, "out_rate": r_, "out_rate_moe90": rm, "arrived_abroad_share": cs.count(mb & arr_abr)[0] / bp, "n_records": int(mb.sum())} put(f"out_rate_{b}", f"Out-migration rate, {b}", r_, rm, (mb & lived_here).sum(), "share") # Groups: who arrives, who leaves, net and rates by characteristic. levels = {"age": h["age_band"].to_numpy(), "income": h["income_band"].to_numpy(), "race": h["race"].to_numpy(), "degree": np.where(h["adult25"], np.where(h["ba_plus"], "ba_plus", "no_degree"), "na"), "children": np.where(h["household"], np.where(h["with_children"], "with_children", "no_children"), "na"), "nativity": np.where(h["foreign_born"], "foreign_born", "us_born"), "tenure": np.where(h["household"], np.where(h["owner"], "owner", "renter"), "na"), "borough": boro} out["groups"] = {} for fk, (flabel, opts) in FACTORS.items(): lv = levels[fk]; out["groups"][fk] = {"label": flabel, "levels": {}} for k, label in opts: mk = lv == k; n_, nm = cs.diff(mk & arr_dom, mk & left); r_, rm = cs.share(mk & lived_here, left) ab, abm = cs.count(mk & arr_abr); ad, adm = cs.count(mk & arr_dom); lf, lfm = cs.count(mk & left) out["groups"][fk]["levels"][k] = {"label": label, "population": cs.count(mk & res)[0], "arrived_domestic": ad, "arrived_domestic_moe90": adm, "arrived_abroad": ab, "arrived_abroad_moe90": abm, "left": lf, "left_moe90": lfm, "net_domestic": n_, "net_domestic_moe90": nm, "out_rate": r_, "out_rate_moe90": rm, "in_rate_domestic": cs.share(mk & res, arr_dom)[0], "in_rate_abroad": cs.share(mk & res, arr_abr)[0], "n_records": int(mk.sum()), "n_left": int((mk & left).sum()), "n_arrived_domestic": int((mk & arr_dom).sum())} put(f"net_{fk}_{k}", f"Net domestic migration, {label}", n_, nm, (mk & (arr_dom | left)).sum(), "count") put(f"out_rate_{fk}_{k}", f"Out-migration rate, {label}", r_, rm, (mk & lived_here).sum(), "share") # Profiles: what each flow looks like. traits = {"under_18": h["AGEP"] < 18, "age_18_34": h["AGEP"].between(18, 34), "age_65_plus": h["AGEP"] >= 65, "ba_plus": h["ba_plus"], "employed": h["employed"], "income_150k_plus": h["hh_income"] >= 150_000, "income_lt50k": h["hh_income"] < 50_000, "foreign_born": h["foreign_born"], "noncitizen": h["noncitizen"], "white": h["race"] == "White", "black": h["race"] == "Black", "hispanic": h["race"] == "Hispanic", "asian": h["race"] == "Asian", "with_children": h["with_children"], "owner": h["owner"], "student": h["student"]} bases = {"ba_plus": h["adult25"].to_numpy(), "employed": h["adult16"].to_numpy(), "income_150k_plus": h["household"].to_numpy(), "income_lt50k": h["household"].to_numpy(), "with_children": h["household"].to_numpy(), "owner": h["household"].to_numpy(), "student": (h["AGEP"] >= 18).to_numpy()} out["profile"] = {} for k, label in STATUSES + [("all_residents", "All residents")]: mk = res if k == "all_residents" else (status == k); prof = {"label": label, "n_records": int(mk.sum())} for t, (tl, _) in zip([x[0] for x in TRAITS], TRAITS): base = mk & bases.get(t, np.ones(len(h), bool)); v, vm = cs.share(base, traits[t].to_numpy()) prof[t] = {"share": v, "moe90": vm, "label": tl} prof["median_age"] = cs.median(mk, "AGEP"); prof["median_hh_income"] = cs.median(mk & h["household"].to_numpy(), "hh_income") prof["median_hh_income_adults"] = cs.median(mk & h["household"].to_numpy() & h["adult25"].to_numpy(), "hh_income") out["profile"][k] = prof # Arrivals from other states and from abroad share little; give each origin a profile too. out["origin_profile"] = {} for k, label in ORIGINS: mk = arrived & (origin == k) if mk.sum() < MIN_RECORDS: continue out["origin_profile"][k] = {"label": label, "n_records": int(mk.sum()), "median_age": cs.median(mk, "AGEP"), "median_hh_income": cs.median(mk & h["household"].to_numpy(), "hh_income"), "ba_plus": cs.share(mk & h["adult25"].to_numpy(), h["ba_plus"].to_numpy())[0], "age_18_34": cs.share(mk, h["AGEP"].between(18, 34).to_numpy())[0], "under_18": cs.share(mk, (h["AGEP"] < 18).to_numpy())[0], "employed": cs.share(mk & h["adult16"].to_numpy(), h["employed"].to_numpy())[0], "student": cs.share(mk & (h["AGEP"] >= 18).to_numpy(), h["student"].to_numpy())[0]} out["destination_profile"] = {} for k, label in DESTINATIONS: mk = left & (dest == k) if mk.sum() < MIN_RECORDS: continue out["destination_profile"][k] = {"label": label, "n_records": int(mk.sum()), "median_age": cs.median(mk, "AGEP"), "median_hh_income": cs.median(mk & h["household"].to_numpy(), "hh_income"), "ba_plus": cs.share(mk & h["adult25"].to_numpy(), h["ba_plus"].to_numpy())[0], "with_children": cs.share(mk & h["household"].to_numpy(), h["with_children"].to_numpy())[0], "owner": cs.share(mk & h["household"].to_numpy(), h["owner"].to_numpy())[0], "under_18": cs.share(mk, (h["AGEP"] < 18).to_numpy())[0], "income_150k_plus": cs.share(mk & h["household"].to_numpy(), (h["hh_income"] >= 150_000).to_numpy())[0]} # Leavers: still working in the city. wk = left & h["works_in_nyc"].notna().to_numpy() out["leavers_work"] = {"works_in_nyc": cs.share(wk, h["works_in_nyc"].fillna(False).to_numpy())[0], "wfh": cs.share(wk, h["wfh"].to_numpy())[0], "n_records": int(wk.sum()), "suburbs_works_in_nyc": cs.share(wk & (dest == "suburbs"), h["works_in_nyc"].fillna(False).to_numpy())[0]} # Arrivals from abroad: where they were born, citizenship, age. wr = h["world_region"].to_numpy(); fr = h["from_region"].to_numpy() out["abroad"] = {"regions": {}, "from_regions": {}, "citizenship": {}, "top_countries": [], "top_countries_lived": []} for k, label in WORLD_REGIONS + [("puerto_rico", "Puerto Rico and the island areas")]: mk = arr_abr & (fr == k); c, cm = cs.count(mk) out["abroad"]["from_regions"][k] = {"label": label, "count": c, "moe90": cm, "share": c / out["flows"]["arrived_abroad"]["count"], "n_records": int(mk.sum())} mnames = migsp_names(); msp = pd.to_numeric(h["MIGSP"], errors="coerce").to_numpy(); pob_ = pd.to_numeric(h["POBP"], errors="coerce").to_numpy() topm = pd.Series(cs.w[arr_abr & (msp >= 100)]).groupby(msp[arr_abr & (msp >= 100)]).sum().sort_values(ascending=False).head(10) out["abroad"]["top_countries_lived"] = [{"code": int(c), "country": mnames.get(int(c), str(int(c))), "count": float(v), "n_records": int((arr_abr & (msp == c)).sum())} for c, v in topm.items()] out["abroad"]["lived_in_birth_country_share"] = cs.share(arr_abr & (msp >= 100) & (pob_ >= 100), msp == pob_)[0] for k, label in WORLD_REGIONS: mk = arr_abr & (wr == k); c, cm = cs.count(mk) out["abroad"]["regions"][k] = {"label": label, "count": c, "moe90": cm, "share": c / out["flows"]["arrived_abroad"]["count"], "n_records": int(mk.sum())} us_born_abroad = arr_abr & (wr == "us") out["abroad"]["regions"]["us"] = {"label": "US-born returning", "count": cs.count(us_born_abroad)[0], "moe90": cs.count(us_born_abroad)[1], "share": cs.count(us_born_abroad)[0] / out["flows"]["arrived_abroad"]["count"], "n_records": int(us_born_abroad.sum())} for k, mk in [("us_citizen_born", arr_abr & ~h["foreign_born"].to_numpy()), ("naturalized", arr_abr & h["naturalized"].to_numpy()), ("noncitizen", arr_abr & h["noncitizen"].to_numpy())]: c, cm = cs.count(mk); out["abroad"]["citizenship"][k] = {"count": c, "moe90": cm, "share": c / out["flows"]["arrived_abroad"]["count"]} names = pobp_names(); pob = pd.to_numeric(h["POBP"], errors="coerce").to_numpy() top = pd.Series(cs.w[arr_abr & (pob >= 100)]).groupby(pob[arr_abr & (pob >= 100)]).sum().sort_values(ascending=False).head(10) out["abroad"]["top_countries"] = [{"code": int(c), "country": names.get(int(c), str(int(c))), "count": float(v), "n_records": int((arr_abr & (pob == c)).sum())} for c, v in top.items()] out["abroad"]["student_share_18_plus"] = cs.share(arr_abr & (h["AGEP"] >= 18).to_numpy(), h["student"].to_numpy())[0] out["abroad"]["employed_share_16_plus"] = cs.share(arr_abr & h["adult16"].to_numpy(), h["employed"].to_numpy())[0] out["foreign_born_share"] = cs.share(res, h["foreign_born"].to_numpy()) return out, rows # ------------------------------------------------------------------ reproduction def reproduction(pub: dict, years: dict, censuses: dict) -> pd.DataFrame: rows = [] for y in ("2019", "2021", "2022", "2023", "2024"): if y not in pub or int(y) not in censuses: continue cs = censuses[int(y)]; h = cs.h; res = h["resident"].to_numpy(); status = h["status"].to_numpy() ours = {"total": cs.count(res), "same_house": cs.count(status == "stayed"), "same_county": cs.count(h["same_county"].to_numpy()), "diff_county_same_state": cs.count(h["diff_county_same_state"].to_numpy()), "different_state": cs.count((status == "arrived_domestic") & ~h["same_county"].to_numpy() & ~h["diff_county_same_state"].to_numpy()), "abroad": cs.count(status == "arrived_abroad")} labels = {"total": "Residents 1 and over", "same_house": "Same house one year ago", "same_county": "Moved within the same county (borough)", "diff_county_same_state": "Moved from a different county in New York State (other boroughs included)", "different_state": "Moved from a different state", "abroad": "Moved from abroad"} for k, lab in labels.items(): # The Bureau's own test for a difference between two estimates: the margin on # the difference is the root of the summed squared margins. The published # table comes from the full ACS sample and the PUMS is a subsample of it, so # the two are not independent, but this is the published guidance. p = pub[y][k]; o, om = ours[k]; band = float(np.sqrt(p["moe90"] ** 2 + om ** 2)) rows.append({"statistic": f"{lab}, NYC, {y}", "published": p["estimate"], "published_moe90": p["moe90"], "reproduced": o, "reproduced_moe90": om, "difference": o - p["estimate"], "difference_moe90": band, "within_margin": bool(abs(o - p["estimate"]) <= band), "kind": "count"}) return pd.DataFrame(rows) # ------------------------------------------------------------------ administrative series def components() -> dict: """Births, deaths and net migration for the five counties, July to July, two vintages spliced.""" def load(path, years, vint): d = pd.read_csv(path, encoding="latin-1", dtype={"STATE": str, "COUNTY": str}) n = d[(d["STATE"] == "36") & d["COUNTY"].isin(NYC_COUNTIES)] out = [] for y in years: row = {"year": y, "vintage": vint, "population": int(n[f"POPESTIMATE{y}"].sum()), "births": int(n[f"BIRTHS{y}"].sum()), "deaths": int(n[f"DEATHS{y}"].sum()), "international": int(n[f"INTERNATIONALMIG{y}"].sum()), "domestic": int(n[f"DOMESTICMIG{y}"].sum())} row["natural"] = row["births"] - row["deaths"]; row["net_migration"] = row["international"] + row["domestic"] row["change"] = int(n[f"NPOPCHG{y}"].sum()) if f"NPOPCHG{y}" in n.columns else int(n[f"NPOPCHG_{y}"].sum()) row["boroughs"] = {NYC_COUNTIES[c]: {"population": int(r[f"POPESTIMATE{y}"]), "international": int(r[f"INTERNATIONALMIG{y}"]), "domestic": int(r[f"DOMESTICMIG{y}"])} for c, r in n.set_index("COUNTY").iterrows()} out.append(row) return out v20 = load(RAW / "popest" / "co-est2020-alldata.csv", range(2011, 2021), "2020") v25 = load(RAW / "popest" / "co-est2025-alldata.csv", range(2021, 2026), "2025") v25_2020_partial = load(RAW / "popest" / "co-est2025-alldata.csv", [2020], "2025")[0] d24 = pd.read_csv(RAW / "popest" / "co-est2024-alldata.csv", encoding="latin-1", dtype={"STATE": str, "COUNTY": str}) n24 = d24[(d24["STATE"] == "36") & d24["COUNTY"].isin(NYC_COUNTIES)] v24 = {y: {"international": int(n24[f"INTERNATIONALMIG{y}"].sum()), "domestic": int(n24[f"DOMESTICMIG{y}"].sum()), "population": int(n24[f"POPESTIMATE{y}"].sum())} for y in range(2021, 2025)} return {"note": "July-to-July years. 2011-2020 from Vintage 2020 (built on the 2010 census), 2021-2025 from Vintage 2025 (built on the 2020 census); the two vintages do not splice: " "the 2020 census counted about 8.80 million residents (the Vintage 2025 estimates base) against the 8,253,213 the Vintage 2020 series carried for July 2020. The 2020 row of Vintage 2025 covers April to July 2020 only and is not used.", "years": v20 + v25, "vintage2025_april_july_2020": v25_2020_partial, "vintage2024_overlap": v24, "estimates_base_2020": int(pd.read_csv(RAW / "popest" / "co-est2025-alldata.csv", encoding="latin-1", dtype={"STATE": str, "COUNTY": str}).query("STATE == '36' and COUNTY in @NYC_COUNTIES")["ESTIMATESBASE2020"].sum())} def irs_flows() -> dict: """Returns, individuals and AGI moving into and out of the five counties, by filing-year pair.""" pairs = [("1516", "2015–16"), ("1617", "2016–17"), ("1718", "2017–18"), ("1819", "2018–19"), ("1920", "2019–20"), ("2021", "2020–21"), ("2122", "2021–22"), ("2223", "2022–23")] out = {"note": "IRS Statistics of Income county-to-county migration: filers whose address county changed between two filing years; n1 returns, n2 exemptions (roughly people), AGI in thousands of dollars. " "Flows between the five boroughs are removed so the totals are for the city. Filers only: non-filers, and most people arriving from abroad, are not seen.", "years": []} nyc = set(NYC_COUNTIES) for code, label in pairs: o = pd.read_csv(RAW / "irs_migration" / f"countyoutflow{code}.csv", dtype={"y1_statefips": str, "y1_countyfips": str, "y2_statefips": str, "y2_countyfips": str}, encoding="latin-1") i = pd.read_csv(RAW / "irs_migration" / f"countyinflow{code}.csv", dtype={"y1_statefips": str, "y1_countyfips": str, "y2_statefips": str, "y2_countyfips": str}, encoding="latin-1") for d in (o, i): for c in ("n1", "n2", "agi"): d[c] = pd.to_numeric(d[c], errors="coerce").fillna(0) # The 2020-21 release writes its FIPS codes unpadded ("36", "5"); pad every year the same way. for c in ("y1_statefips", "y2_statefips"): d[c] = d[c].astype(str).str.strip().str.zfill(2) for c in ("y1_countyfips", "y2_countyfips"): d[c] = d[c].astype(str).str.strip().str.zfill(3) so = o[(o["y1_statefips"] == "36") & o["y1_countyfips"].isin(nyc)]; si = i[(i["y2_statefips"] == "36") & i["y2_countyfips"].isin(nyc)] tot_o = so[(so["y2_statefips"] == "96")][["n1", "n2", "agi"]].sum(); tot_i = si[(si["y1_statefips"] == "96")][["n1", "n2", "agi"]].sum() intra_o = so[(so["y2_statefips"] == "36") & so["y2_countyfips"].isin(nyc) & (so["y2_countyfips"] != so["y1_countyfips"])][["n1", "n2", "agi"]].sum() intra_i = si[(si["y1_statefips"] == "36") & si["y1_countyfips"].isin(nyc) & (si["y1_countyfips"] != si["y2_countyfips"])][["n1", "n2", "agi"]].sum() non = so[(so["y2_statefips"] == "36") & (so["y2_countyfips"] == so["y1_countyfips"])][["n1", "n2", "agi"]].sum() foreign_o = so[so["y2_statefips"] == "98"][["n1", "n2", "agi"]].sum(); foreign_i = si[si["y1_statefips"] == "98"][["n1", "n2", "agi"]].sum() out_ = tot_o - intra_o; in_ = tot_i - intra_i # Destinations and origins by county, outside the city. # State codes 57-59 are the IRS's "other flows" aggregates (suppressed small pairs), not states. dests = so[so["y2_statefips"].str.match(r"^(0[1-9]|[1-4][0-9]|5[0-6])$") & ~((so["y2_statefips"] == "36") & so["y2_countyfips"].isin(nyc))].copy() dests["fips"] = dests["y2_statefips"] + dests["y2_countyfips"] origs = si[si["y1_statefips"].str.match(r"^(0[1-9]|[1-4][0-9]|5[0-6])$") & ~((si["y1_statefips"] == "36") & si["y1_countyfips"].isin(nyc))].copy() origs["fips"] = origs["y1_statefips"] + origs["y1_countyfips"] def by_state(d, col): g = d.groupby(col)[["n1", "n2", "agi"]].sum().sort_values("n2", ascending=False) return [{"state": STATE_NAMES.get(int(s), s), "returns": int(r["n1"]), "individuals": int(r["n2"]), "agi_thousands": int(r["agi"])} for s, r in g.head(10).iterrows()] def suburbs(d): m = d["fips"].map(SUBURB_COUNTIES_IRS); g = d[m.notna()].groupby(m[m.notna()])[["n1", "n2", "agi"]].sum() return {k: {"returns": int(r["n1"]), "individuals": int(r["n2"]), "agi_thousands": int(r["agi"])} for k, r in g.iterrows()} row = {"pair": code, "label": label, "out": {"returns": int(out_["n1"]), "individuals": int(out_["n2"]), "agi_thousands": int(out_["agi"]), "agi_per_return": float(out_["agi"] * 1000 / out_["n1"])}, "in": {"returns": int(in_["n1"]), "individuals": int(in_["n2"]), "agi_thousands": int(in_["agi"]), "agi_per_return": float(in_["agi"] * 1000 / in_["n1"])}, "non_migrants": {"returns": int(non["n1"]), "individuals": int(non["n2"]), "agi_thousands": int(non["agi"]), "agi_per_return": float(non["agi"] * 1000 / non["n1"])}, "foreign": {"out_returns": int(foreign_o["n1"]), "in_returns": int(foreign_i["n1"])}, "net_individuals": int(in_["n2"] - out_["n2"]), "net_returns": int(in_["n1"] - out_["n1"]), "net_agi_thousands": int(in_["agi"] - out_["agi"]), "destination_states": by_state(dests, "y2_statefips"), "origin_states": by_state(origs, "y1_statefips"), "suburbs_out": suburbs(dests), "suburbs_in": suburbs(origs)} row["out"]["individuals_per_return"] = row["out"]["individuals"] / row["out"]["returns"] row["other_flows_out_individuals"] = int(so[so["y2_statefips"].isin(["57", "58", "59"])]["n2"].sum()) row["other_flows_in_individuals"] = int(si[si["y1_statefips"].isin(["57", "58", "59"])]["n2"].sum()) out["years"].append(row) return out def shelter_series() -> dict: d = pd.read_csv(RAW / "homeless" / "nyc_dhs_daily_report.csv"); d["date"] = pd.to_datetime(d["Date of Census"]); d = d.sort_values("date") d["month"] = d["date"].dt.to_period("M") cols = {"Total Individuals in Shelter": "total", "Total Single Adults in Shelter": "single_adults", "Total Individuals in Families with Children in Shelter ": "families_with_children", "Individuals in Adult Families in Shelter": "adult_families", "Total Children in Shelter": "children"} cols = {k: v for k, v in cols.items() if k in d.columns} m = d.groupby("month")[list(cols)].mean().rename(columns=cols).round() m = m[m.index >= pd.Period("2021-01", "M")] months = [{"month": str(k), **{c: float(r[c]) for c in m.columns}, "days": int((d["month"] == k).sum())} for k, r in m.iterrows()] y = d.groupby(d["date"].dt.year)[list(cols)].mean().rename(columns=cols).round() peak = m["total"].idxmax() return {"note": "Daily census of the Department of Homeless Services shelter system, averaged by month. It counts everyone in DHS shelters and does not separate asylum seekers; " "the humanitarian centres run by other agencies are not in it.", "months": months, "last_date": d["date"].max().strftime("%Y-%m-%d"), "annual": {int(k): {c: float(r[c]) for c in y.columns} for k, r in y.iterrows()}, "peak_month": str(peak), "peak_total": float(m.loc[peak, "total"]), "latest_month": str(m.index[-1]), "latest_total": float(m["total"].iloc[-1]), "early_2022": float(m.loc[pd.Period("2022-03", "M"), "total"])} def cps_nativity() -> dict | None: f = INTERIM / "cps_nativity_summary.json" if not f.exists(): return None s = json.loads(f.read_text(encoding="utf-8")); months = s["months"] periods = {"2024_jun_dec": ("June–December 2024", lambda m: m["label"] < "2025"), "2025_jan_aug": ("January–August 2025", lambda m: "2025-01" <= m["label"] <= "2025-08"), "2025_sep_dec": ("September–December 2025", lambda m: "2025-09" <= m["label"] <= "2025-12"), "2026_jan_aug": ("January–August 2026", lambda m: m["label"] >= "2026")} out = {"months": months, "periods": {}, "first_month": months[0]["label"], "last_month": months[-1]["label"], "notes": s.get("notes")} for k, (label, f_) in periods.items(): ms = [m for m in months if f_(m)] if not ms: continue block = {"label": label, "months": len(ms), "n": int(sum(m["n"] for m in ms))} for stat in ("foreign_born_share", "noncitizen_share", "foreign_born", "population"): v = np.array([m[stat] for m in ms], float); block[stat] = float(v.mean()); block[stat + "_moe90"] = float(Z90 * v.std(ddof=1) / np.sqrt(len(v))) if len(v) > 1 else float("nan") v = np.array([m["national"]["foreign_born_share"] for m in ms], float); block["us_foreign_born_share"] = float(v.mean()) out["periods"][k] = block return out # ------------------------------------------------------------------ figures def _mpl(): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"font.family": "sans-serif", "font.size": 10, "axes.edgecolor": GRID, "axes.labelcolor": INK, "xtick.color": INK, "ytick.color": INK, "text.color": INK, "svg.fonttype": "none"}) return plt def _style(ax, title, sub): ax.set_title(title, loc="left", fontsize=12, color="#1a1a1a", pad=22) ax.annotate(sub, xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left") for s in ("top", "right", "left"): ax.spines[s].set_visible(False) ax.grid(axis="x", color=GRID, linewidth=0.8); ax.set_axisbelow(True); ax.tick_params(length=0) def _save(fig, path: Path): path.parent.mkdir(parents=True, exist_ok=True) fig.savefig(path, format="svg", bbox_inches="tight"); (path.parent / "pdf").mkdir(exist_ok=True) fig.savefig(path.parent / "pdf" / path.with_suffix(".pdf").name, format="pdf", bbox_inches="tight") def _grouped(rows, path, title, sub, labels, fmt="{:,.0f}", height=None, xlabel=None, colors=(BLUE, GREEN, ORANGE)): """rows: (label, [(value, moe), ...]) with one value per series label; negative values allowed.""" plt = _mpl(); n = len(rows); k = len(labels) fig, ax = plt.subplots(figsize=(FIG_W, height or (1.3 + 0.28 * k * n)), dpi=100); fig.patch.set_alpha(0) y = np.arange(n)[::-1]; hgt = 0.8 / k vals = [v for r in rows for v, m in r[1] if v == v]; moes = [m if m == m else 0 for r in rows for v, m in r[1]] top = max(abs(v) + m for v, m in zip(vals, moes)) if vals else 1; neg = min(vals) < 0 if vals else False for j, lab in enumerate(labels): v = np.array([r[1][j][0] for r in rows], float); m = np.array([r[1][j][1] if r[1][j][1] == r[1][j][1] else 0 for r in rows], float) yy = y + ((k - 1) / 2 - j) * hgt ax.barh(yy, v, height=hgt * 0.92, color=colors[j % len(colors)], label=lab) ax.errorbar(v, yy, xerr=m, fmt="none", ecolor=TICK, elinewidth=1.1) for yi, vi, mi in zip(yy, v, m): if vi != vi: continue x = vi + mi + top * 0.015 if vi >= 0 else vi - mi - top * 0.015 ax.text(x, yi, fmt.format(vi), va="center", ha="left" if vi >= 0 else "right", fontsize=8, color=INK) ax.set_yticks(y); ax.set_yticklabels([r[0] for r in rows], fontsize=9) ax.set_xlim(-top * 1.25 if neg else 0, top * 1.25) if neg: ax.axvline(0, color=TICK, linewidth=0.8) if xlabel: ax.set_xlabel(xlabel, fontsize=9, color=INK) _style(ax, title, sub); ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.06 if not xlabel else -0.12), ncol=min(k, 3), frameon=False, fontsize=9, labelcolor=INK) fig.tight_layout(); _save(fig, path); plt.close(fig) def _triple(years: dict, getter, keys, path, title, sub, fmt="{:,.0f}", height=None, xlabel=None): rows = [(label, [getter(years[str(y)], k) for y in FOCUS]) for k, label in keys] _grouped(rows, path, title, sub, [str(y) for y in FOCUS], fmt=fmt, height=height, xlabel=xlabel) def fig_flows(years: dict, path): plt = _mpl(); ys = [years[str(y)] for y in YEARS if str(y) in years] fig, ax = plt.subplots(figsize=(FIG_W, 3.6), dpi=100); fig.patch.set_alpha(0) x = [r["year"] for r in ys] series = [("arrived_abroad", "Arrived from abroad", GREEN), ("arrived_domestic", "Arrived from elsewhere in the US", BLUE), ("left", "Left for elsewhere in the US", ORANGE)] for k, lab, col in series: v = np.array([r["flows"][k]["count"] / 1000 for r in ys]); m = np.array([r["flows"][k]["moe90"] / 1000 for r in ys]) ax.plot(x, v, color=col, linewidth=2, marker="o", markersize=4, label=lab); ax.fill_between(x, v - m, v + m, color=col, alpha=0.15, linewidth=0) net = np.array([r["net_domestic"]["estimate"] / 1000 for r in ys]) ax.plot(x, net, color=TICK, linewidth=1.4, linestyle="--", marker="o", markersize=3, label="Net domestic (arrivals minus leavers)") ax.axhline(0, color=GRID, linewidth=0.8) ax.set_xticks(x); ax.set_xticklabels([str(xi) if xi != 2020 else "2020*" for xi in x], fontsize=8.5) ax.set_yticks([-200, -100, 0, 100, 200, 300]); ax.set_yticklabels(["−200k", "−100k", "0", "100k", "200k", "300k"]) ax.grid(axis="y", color=GRID, linewidth=0.8); ax.grid(axis="x", visible=False) for s in ("top", "right"): ax.spines[s].set_visible(False) ax.set_title("Who arrives, who leaves, year by year", loc="left", fontsize=12, color="#1a1a1a", pad=22) ax.annotate("People a year, thousands (census, one-year files; band: 90% margin; *experimental weights)", xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left") ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.1), ncol=2, frameon=False, fontsize=8.5, labelcolor=INK) fig.tight_layout(); _save(fig, path); plt.close(fig) SHORT = {"abroad": "Abroad", "other_state": "Another state", "suburbs": "The metro suburbs", "rest_ny": "Elsewhere in NY State", "rest_nj_ct": "Elsewhere in NJ or CT", "florida": "Florida", "pennsylvania": "Pennsylvania", "other_northeast": "Rest of the Northeast", "california": "California", "texas": "Texas", "other_south": "Rest of the South", "midwest": "Midwest", "other_west": "Rest of the West"} def fig_origins(years, path): _triple(years, lambda r, k: (r["origins"][k]["count"] / 1000, r["origins"][k]["moe90"] / 1000), [(k, SHORT[k]) for k, _ in ORIGINS], path, "Where the arrivals come from", "People who moved to NYC in the past year, by where they lived a year earlier (census, thousands)", fmt="{:,.0f}k", height=4.6) def fig_destinations(years, path): _triple(years, lambda r, k: (r["destinations"][k]["count"] / 1000, r["destinations"][k]["moe90"] / 1000), [(k, SHORT[k]) for k, _ in DESTINATIONS], path, "Where the leavers go", "People who lived in NYC a year earlier and now live elsewhere in the US, by destination (census, thousands)", fmt="{:,.0f}k", height=7.2) def fig_net_age(years, path): keys = [(k, l) for k, l in FACTORS["age"][1]] _triple(years, lambda r, k: (r["groups"]["age"]["levels"][k]["net_domestic"] / 1000, r["groups"]["age"]["levels"][k]["net_domestic_moe90"] / 1000), keys, path, "Who the city gains and loses to the rest of the country", "Net domestic migration by age: arrivals from elsewhere in the US minus leavers (census, thousands)", fmt="{:+,.0f}k", height=5) def fig_net_income(years, path): keys = [(k, l) for k, l in FACTORS["income"][1]] _triple(years, lambda r, k: (r["groups"]["income"]["levels"][k]["net_domestic"] / 1000, r["groups"]["income"]["levels"][k]["net_domestic_moe90"] / 1000), keys, path, "Net domestic migration by household income", "Arrivals from elsewhere in the US minus leavers, by the income of the household they live in now, 2024 dollars (census, thousands)", fmt="{:+,.0f}k", height=5.4) def fig_rates(years, path): keys = [(k, l) for k, l in FACTORS["age"][1]] _triple(years, lambda r, k: (r["groups"]["age"]["levels"][k]["out_rate"] * 100, r["groups"]["age"]["levels"][k]["out_rate_moe90"] * 100), keys, path, "Who leaves, as a share of who was here", "Share of last year's residents who now live elsewhere in the US, by age (census)", fmt="{:.1f}%", height=5) def fig_profile(y24, path): keys = [(t, l) for t, l in TRAITS if t in ("age_18_34", "ba_plus", "income_150k_plus", "income_lt50k", "foreign_born", "with_children", "owner", "under_18")] rows = [(l, [(y24["profile"][s][t]["share"] * 100, y24["profile"][s][t]["moe90"] * 100) for s in ("arrived_domestic", "left", "stayed")]) for t, l in keys] _grouped(rows, path, "Arrivals, leavers and everyone else, 2024", "Share of each group with the trait (census, 2024 one-year file)", ["Arrived from elsewhere in the US", "Left for elsewhere in the US", "Did not move"], fmt="{:.0f}%", height=6.4) def fig_abroad(years, path): keys = WORLD_REGIONS + [("puerto_rico", "Puerto Rico and the island areas")] _triple(years, lambda r, k: (r["abroad"]["from_regions"][k]["count"] / 1000, r["abroad"]["from_regions"][k]["moe90"] / 1000), keys, path, "Arrivals from abroad, by where they lived a year earlier", "People who lived outside the fifty states a year earlier, by the country or region they lived in (census, thousands)", fmt="{:,.0f}k", height=5.4) def fig_components(comp: dict, path): plt = _mpl(); ys = comp["years"] fig, ax = plt.subplots(figsize=(FIG_W, 3.6), dpi=100); fig.patch.set_alpha(0) x = np.arange(len(ys)); w = 0.27 ax.bar(x - w, [r["natural"] / 1000 for r in ys], width=w, color="#c9a227", label="Births minus deaths") ax.bar(x, [r["international"] / 1000 for r in ys], width=w, color=GREEN, label="Net international migration") ax.bar(x + w, [r["domestic"] / 1000 for r in ys], width=w, color=ORANGE, label="Net domestic migration") ax.plot(x, [r["change"] / 1000 for r in ys], color=TICK, linewidth=1.4, marker="o", markersize=3, label="Total change") ax.axhline(0, color=TICK, linewidth=0.8); ax.axvline(9.5, color=GRID, linewidth=1, linestyle=":") ax.text(9.6, 220, "Vintage 2025 →", fontsize=8, color=INK); ax.text(9.4, 220, "← Vintage 2020", fontsize=8, color=INK, ha="right") ax.set_xticks(x); ax.set_xticklabels([str(r["year"]) for r in ys], fontsize=8, rotation=45) ax.set_yticks([-300, -200, -100, 0, 100, 200]); ax.set_yticklabels(["−300k", "−200k", "−100k", "0", "100k", "200k"]) ax.grid(axis="y", color=GRID, linewidth=0.8); ax.grid(axis="x", visible=False) for s in ("top", "right"): ax.spines[s].set_visible(False) ax.set_title("The Census Bureau's ledger, July to July", loc="left", fontsize=12, color="#1a1a1a", pad=22) ax.annotate("Components of population change, NYC, thousands (population estimates; administrative, no margins)", xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left") ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.18), ncol=2, frameon=False, fontsize=8.5, labelcolor=INK) fig.tight_layout(); _save(fig, path); plt.close(fig) def fig_irs(irs: dict, path): plt = _mpl(); ys = irs["years"] fig, ax = plt.subplots(figsize=(FIG_W, 3.2), dpi=100); fig.patch.set_alpha(0) x = np.arange(len(ys)) for k, lab, col in (("out", "Filed from the city, then from elsewhere", ORANGE), ("in", "Filed from elsewhere, then from the city", BLUE), ("non_migrants", "Filed from the city both years", INK)): ax.plot(x, [r[k]["agi_per_return"] / 1000 for r in ys], color=col, linewidth=2, marker="o", markersize=4, label=lab) ax.set_xticks(x); ax.set_xticklabels([r["label"] for r in ys], fontsize=8.5, rotation=30) ax.set_ylim(0, None); ax.yaxis.set_major_formatter(lambda v, _: f"${v:,.0f}k") ax.grid(axis="y", color=GRID, linewidth=0.8); ax.grid(axis="x", visible=False) for s in ("top", "right"): ax.spines[s].set_visible(False) ax.set_title("What the tax returns say about income", loc="left", fontsize=12, color="#1a1a1a", pad=22) ax.annotate("Average adjusted gross income per return, by whether the return moved (IRS county migration data; nominal dollars)", xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left") ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.22), ncol=1, frameon=False, fontsize=8.5, labelcolor=INK) fig.tight_layout(); _save(fig, path); plt.close(fig) def fig_shelter(sh: dict, cps: dict | None, path): plt = _mpl(); ms = sh["months"] fig, ax = plt.subplots(figsize=(FIG_W, 3.2), dpi=100); fig.patch.set_alpha(0) x = np.arange(len(ms)); ax.bar(x, [m["total"] / 1000 for m in ms], color=BLUE, width=0.85, label="People in DHS shelters, monthly average") ticks = [i for i, m in enumerate(ms) if m["month"].endswith("-01") or i == 0] ax.set_xticks(ticks); ax.set_xticklabels([ms[i]["month"][:4] for i in ticks], fontsize=8.5) ax.set_ylim(0, 100); ax.set_yticks([0, 25, 50, 75, 100]); ax.set_yticklabels(["0", "25k", "50k", "75k", "100k"]) ax.grid(axis="y", color=GRID, linewidth=0.8); ax.grid(axis="x", visible=False) for s in ("top", "right"): ax.spines[s].set_visible(False) ax.set_title("The shelter census, to September 2026", loc="left", fontsize=12, color="#1a1a1a", pad=22) ax.annotate(f"People in the Department of Homeless Services shelter system, monthly average, January 2021 to {sh['last_date']} (administrative)", xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left") fig.tight_layout(); _save(fig, path); plt.close(fig) def fig_cps_nativity(cps: dict, path): plt = _mpl(); ms = cps["months"] fig, ax = plt.subplots(figsize=(FIG_W, 3.2), dpi=100); fig.patch.set_alpha(0) x = np.arange(len(ms)); v = np.array([m["foreign_born_share"] * 100 for m in ms]) ax.bar(x, v, color=BLUE, width=0.85, label="Foreign-born share of NYC residents") us = np.array([m["national"]["foreign_born_share"] * 100 for m in ms]); ax.plot(x, us, color=ORANGE, linewidth=1.6, label="United States") ax.set_xticks(x[::3]); ax.set_xticklabels([ms[i]["label"] for i in range(0, len(ms), 3)], fontsize=8.5) ax.set_ylim(0, 50); ax.set_yticks([0, 10, 20, 30, 40, 50]); ax.set_yticklabels(["0%", "10%", "20%", "30%", "40%", "50%"]) ax.grid(axis="y", color=GRID, linewidth=0.8); ax.grid(axis="x", visible=False) for sp in ("top", "right"): ax.spines[sp].set_visible(False) ax.set_title("The foreign-born share, month by month, to 2026", loc="left", fontsize=12, color="#1a1a1a", pad=22) ax.annotate("Residents of the five boroughs born outside the US (Current Population Survey, monthly; no October 2025 file)", xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left") ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.1), ncol=2, frameon=False, fontsize=9, labelcolor=INK) fig.tight_layout(); _save(fig, path); plt.close(fig) # ------------------------------------------------------------------- companion data def companion_dataset(frames: list[pd.DataFrame], rng: np.random.Generator) -> pd.DataFrame: parts = [] for h in frames: parts.append(pd.DataFrame({ "year": h["year"], "weight": h["PWGTP"].astype(int), "status": h["status"], "borough": h["borough"], "origin": h["origin"].replace("", np.nan), "origin_region": h["origin_region"].replace("", np.nan), "destination": h["destination"].replace("", np.nan), "destination_region": h["dest_region"].replace("", np.nan), "destination_state": h["dest_state"], "age": h["AGEP"], "age_band": h["age_band"], "sex": h["SEX"].map({1: "male", 2: "female"}), "race_ethnicity": h["race"], "foreign_born": h["foreign_born"], "noncitizen": h["noncitizen"], "region_of_birth": h["world_region"], "region_lived_in_a_year_ago": h["from_region"].replace("", np.nan), "bachelors_or_higher_25plus": h["ba_plus"].where(h["adult25"]), "employed_16plus": h["employed"].where(h["adult16"]), "student_18plus": h["student"].where(h["AGEP"] >= 18), "household_income_2024": h["hh_income"].round(0), "income_band": h["income_band"].replace("na", np.nan), "owner": h["owner"].where(h["household"]), "with_children": h["with_children"].where(h["household"]), "works_in_nyc": h["works_in_nyc"], "works_from_home": h["wfh"], })) df = pd.concat(parts, ignore_index=True) # Booleans as 0/1 and integer ids: the Worker serves assets up to 25 MiB and the # three years of records are 217,000 rows. for c in df.columns: if df[c].dtype == bool or df[c].dropna().isin([True, False]).all() and df[c].notna().any() and c not in ("year", "weight", "age"): df[c] = df[c].map({True: 1, False: 0}).astype("Int64") df.insert(0, "row_id", rng.permutation(len(df)) + 1) return df.sort_values("row_id").reset_index(drop=True) # ------------------------------------------------------------------------- main def main(argv=None) -> int: ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--seed", type=int, default=20240101) ap.add_argument("--years", nargs="*", type=int, default=YEARS) args = ap.parse_args(argv) rng = np.random.default_rng(args.seed); OUT.mkdir(parents=True, exist_ok=True) print("published tables ...") pub = published_tables() for y in ("2019", "2024"): p = pub[y]; print(f" {y}: residents {p['total']['estimate']:,.0f}, from another state {p['different_state']['estimate']:,.0f}, from abroad {p['abroad']['estimate']:,.0f}") years: dict[int, dict] = {}; pums_rows = []; frames = {}; censuses = {} for y in args.years: try: h, rw = load_year(y) except FileNotFoundError as e: print(f" {y}: not on disk ({e}), skipped"); continue cs = Census(h, rw); res, rows = year_estimates(cs); years[y] = res; pums_rows += rows; censuses[y] = cs if y in FOCUS: frames[y] = h f = res["flows"] print(f" {y}: residents {res['population']['estimate']:,.0f}; arrived abroad {f['arrived_abroad']['count']:,.0f}, arrived US {f['arrived_domestic']['count']:,.0f}, " f"left {f['left']['count']:,.0f} (n {f['left']['n_records']:,}); net domestic {res['net_domestic']['estimate']:+,.0f} ± {res['net_domestic']['moe90']:,.0f}; out-rate {res['out_rate']['estimate']:.2%}") print("reproduction ...") repro = reproduction(pub, years, censuses) print(repro[["statistic", "published", "reproduced", "within_margin"]].to_string(index=False)) print("population estimates, IRS, shelter, CPS ...") comp = components(); irs = irs_flows(); sh = shelter_series(); cps = cps_nativity() for r in comp["years"][-3:]: print(f" popest {r['year']}: international {r['international']:+,}, domestic {r['domestic']:+,}, change {r['change']:+,}") for r in irs["years"][-2:]: print(f" IRS {r['label']}: out {r['out']['individuals']:,} people (${r['out']['agi_per_return']:,.0f}/return), in {r['in']['individuals']:,} (${r['in']['agi_per_return']:,.0f}), stayers ${r['non_migrants']['agi_per_return']:,.0f}") print(f" shelter: peak {sh['peak_total']:,.0f} in {sh['peak_month']}, latest {sh['latest_total']:,.0f} in {sh['latest_month']}") if cps: for k, b in cps["periods"].items(): print(f" CPS {b['label']}: foreign-born {b['foreign_born_share']:.1%} ± {b['foreign_born_share_moe90']:.1%}, non-citizen {b['noncitizen_share']:.1%}") y15, y19, y24 = (years.get(y) for y in FOCUS) change = {} if y15 and y24: for k, _ in STATUSES: change[f"{k}_2015_2024"] = y24["flows"][k]["count"] - y15["flows"][k]["count"] change["net_domestic_2015_2024"] = y24["net_domestic"]["estimate"] - y15["net_domestic"]["estimate"] change["out_rate_points_2015_2024"] = (y24["out_rate"]["estimate"] - y15["out_rate"]["estimate"]) * 100 res = {"article": SLUG, "as_of": "2026-09-11", "definitions": {"universe": "NYC residents 1 and over (household and group quarters) in each one-year file; leavers are residents of the five boroughs a year earlier now living elsewhere in the US, from the national file", "dollar_year": 2024, "deflator": "New York-metro CPI-U (all items) to 2024 dollars after the Bureau's within-year factor", "census_margin": "90% (successive difference replication, 80 replicate weights)", "origins": dict(ORIGINS), "destinations": dict(DESTINATIONS), "statuses": dict(STATUSES), "factors": {k: v[0] for k, v in FACTORS.items()}, "suburbs": SUBURB_REGIONS, "abroad": "outside the fifty states and DC one year ago, so Puerto Rico and the island areas count as abroad, as in table B07001"}, "published": pub, "years": {str(y): r for y, r in years.items()}, "change": change, "components": comp, "irs": irs, "shelter": sh, "cps_nativity": cps, "context": {"dcp_july_2026": {"population_july_2025": 8_585_000, "population_july_2024_revised": 8_597_000, "population_july_2024_vintage_2024": 8_478_000, "census_2020": 8_804_000, "source": "dcp_population_estimates_july_2026"}, "asylum_seekers_mayor_july_2025": {"through_care_since_spring_2022": 237_000, "peak_in_care": 69_000, "peak_month": "2024-01", "in_care_july_2025": 37_000, "weekly_arrivals_peak_may_2023": 4_000, "weekly_arrivals_july_2025": 100, "source": "nyc_mayor_asylum_arrival_center_closure_2025"}, "osc_taxpayer_migration_2026": {"scope": "New York State", "tax_year": 2024, "part_year_filers": 256_164, "moved_in": 121_251, "left": 134_913, "net": -13_662, "married_100k_500k_net_loss": 8_200, "over_500k_left_share": 0.01, "largest_net_loss_year": 2020, "largest_net_loss": -112_458, "source": "osc_taxpayer_migration_2026"}}, "reproduction": repro.to_dict(orient="records"), "seed": args.seed} print("figures ...") ys = res["years"] fig_flows(ys, OUT / "fig1_flows.svg"); fig_origins(ys, OUT / "fig2_origins.svg"); fig_destinations(ys, OUT / "fig3_destinations.svg") fig_net_age(ys, OUT / "fig4_net_age.svg"); fig_net_income(ys, OUT / "fig5_net_income.svg"); fig_rates(ys, OUT / "fig6_rates.svg") fig_profile(ys["2024"], OUT / "fig7_profile.svg"); fig_abroad(ys, OUT / "fig8_abroad.svg"); fig_components(comp, OUT / "fig9_components.svg") fig_irs(irs, OUT / "fig10_irs.svg"); fig_shelter(sh, cps, OUT / "fig11_shelter.svg") if cps: fig_cps_nativity(cps, OUT / "fig12_cps_nativity.svg") print("interactive chart cubes ...") charts = CH.build(censuses, ys, comp, irs, sh, cps) CH.write(OUT / "charts.json", charts) print(f" {len(charts['charts'])} charts, {sum(len(c['cells']) for c in charts['charts'].values())} filter combinations, {(OUT / 'charts.json').stat().st_size / 1024:.0f} KB") companion_dataset([frames[y] for y in FOCUS if y in frames], rng).to_csv(OUT / "data.csv", index=False) pd.DataFrame(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") 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") print(f"outputs -> {OUT.relative_to(P.root)} and {SITE_OUT.relative_to(P.root)}") return 0 if __name__ == "__main__": raise SystemExit(main())