1#!/usr/bin/env python32"""302_build_calibration_targets.py — turn the raw downloads of 01_download_data.py into4 * data/processed/series/*.csv tidy (date,value) time series, one per retrieved series5 * data/processed/calibration_targets.csv one row per calibration target (value, unit, period,6 geography, frequency, source, series_id, url, raw_file,7 retrieved, status, notes)8 * data/processed/bea_io_value_added_2023.csv sector labor shares / intermediate shares (BEA use table)9 * data/processed/epoch_price_performance_fit.csv hardware rows used in the FLOP/s-per-$ trend1011Run: python3 -I scripts/02_build_calibration_targets.py1213Every number written to calibration_targets.csv is computed from a file in data/raw/ listed in14data/raw/manifest.json. If a raw file is missing, the corresponding targets are written with15status=NOT_RETRIEVED and an empty value — nothing is ever filled in from memory.16Raw files are untrusted data: they are only parsed (pandas / openpyxl / xlrd / PyMuPDF), never executed.17"""18from __future__ import annotations1920import csv21import json22import math23import re24import sys25from pathlib import Path2627import numpy as np28import pandas as pd2930ROOT = Path(__file__).resolve().parents[1]31RAW = ROOT / "data" / "raw"32PROC = ROOT / "data" / "processed"33SER = PROC / "series"34SER.mkdir(parents=True, exist_ok=True)35RETRIEVED = "2026-10-27"36AVG_START = 19903738manifest = json.loads((RAW / "manifest.json").read_text()) if (RAW / "manifest.json").exists() else {}39fred_meta = json.loads((RAW / "fred" / "_metadata.json").read_text()) if (RAW / "fred" / "_metadata.json").exists() else {}4041ROWS: list[dict] = []42COLS = ["target_id", "description", "value", "unit", "year_or_period", "geography", "frequency",43 "source", "series_id", "url", "raw_file", "retrieved", "status", "notes"]444546def add(target_id, description, value, unit, period, geography, frequency, source, series_id, url,47 raw_file, notes="", status=None):48 if status is None:49 status = "RETRIEVED" if (value is not None and not (isinstance(value, float) and math.isnan(value))) else "NOT_RETRIEVED"50 if status in ("NOT_RETRIEVED", "INFO"):51 value = ""52 elif isinstance(value, float):53 value = f"{value:.6g}"54 ROWS.append(dict(target_id=target_id, description=description, value=value, unit=unit,55 year_or_period=period, geography=geography, frequency=frequency, source=source,56 series_id=series_id, url=url, raw_file=raw_file, retrieved=RETRIEVED, status=status,57 notes=notes))585960def not_retrieved(target_id, description, url, raw_file, notes, geography="United States", unit="", frequency="", source="", series_id=""):61 add(target_id, description, None, unit, "", geography, frequency, source, series_id, url, raw_file, notes, status="NOT_RETRIEVED")626364def save_series(name: str, s: pd.Series, note: str | None = None) -> None:65 s = s.dropna()66 df = pd.DataFrame({"date": pd.to_datetime(s.index).strftime("%Y-%m-%d"), "value": s.values})67 df.to_csv(SER / f"{name}.csv", index=False)686970# ----------------------------------------------------------------------------71# FRED helpers72# ----------------------------------------------------------------------------73def fred_path(sid: str) -> Path:74 return RAW / "fred" / f"{sid}.csv"757677def fred(sid: str) -> pd.Series | None:78 p = fred_path(sid)79 if not p.exists():80 return None81 df = pd.read_csv(p)82 if list(df.columns) != ["observation_date", sid]:83 raise ValueError(f"unexpected FRED layout in {p}: {df.columns.tolist()}")84 s = pd.Series(pd.to_numeric(df[sid], errors="coerce").values, index=pd.to_datetime(df["observation_date"]), name=sid)85 save_series(f"fred_{sid}", s)86 return s878889def fred_url(sid: str) -> str:90 return f"https://fred.stlouisfed.org/graph/fredgraph.csv?id={sid}"919293def fred_file(sid: str) -> str:94 return f"data/raw/fred/{sid}.csv"959697def fred_desc(sid: str) -> str:98 m = fred_meta.get(sid, {})99 return f"{m.get('title', sid)} [{m.get('units', '?')}; {m.get('frequency', '?')}; vintage last updated {m.get('last_updated', '?')}]"100101102def annual_mean(s: pd.Series, full_years_only: bool = True) -> pd.Series:103 """Calendar-year means. Quarterly series: a year is 'full' with 4 obs, monthly 12, daily >= 240."""104 s = s.dropna()105 g = s.groupby(s.index.year)106 m, n = g.mean(), g.count()107 if full_years_only:108 freq = pd.infer_freq(s.index[:12]) if len(s) >= 12 else None109 spacing = (s.index[1:] - s.index[:-1]).median().days if len(s) > 2 else 365110 # monthly: 11 obs accepted because BLS CPI/average-price series have no October-2025 observation111 need = 4 if 80 <= spacing <= 95 else 11 if 27 <= spacing <= 32 else 240 if spacing <= 5 else 1112 m = m[n >= need]113 m.index = pd.to_datetime([f"{y}-01-01" for y in m.index])114 return m115116117def latest(s: pd.Series):118 s = s.dropna()119 return s.index[-1], float(s.iloc[-1])120121122def period_label(idx, freq: str = "quarterly") -> str:123 """Label an observation date according to the series frequency."""124 ts = pd.Timestamp(idx)125 f = (freq or "").lower()126 if "annual" in f:127 return str(ts.year)128 if "month" in f:129 return ts.strftime("%Y-%m")130 if "dai" in f:131 return ts.strftime("%Y-%m-%d")132 return f"{ts.year}Q{(ts.month - 1) // 3 + 1}"133134135def mean_from(s: pd.Series, start_year: int = AVG_START) -> tuple[float, str]:136 s = s.dropna()137 s = s[s.index.year >= start_year]138 return float(s.mean()), f"{s.index[0].year}-{s.index[-1].year} mean of {len(s)} obs"139140141def ratio_targets(tid, desc, num_id, den_id, unit="ratio", scale=1.0, geography="United States", extra_note="", annual=False):142 """Latest value and AVG_START-latest average of num/den (both FRED). scale multiplies num/den."""143 num, den = fred(num_id), fred(den_id)144 src = "FRED (" + "; ".join(sorted({fred_meta.get(x, {}).get("source", "") for x in (num_id, den_id)} - {""})) + ")"145 url = fred_url(num_id) + " ; " + fred_url(den_id)146 rawf = fred_file(num_id) + " ; " + fred_file(den_id)147 if num is None or den is None:148 not_retrieved(tid + "_latest", desc + ", latest", url, rawf, "FRED file(s) missing", geography, unit, "", src, f"{num_id}/{den_id}")149 not_retrieved(tid + "_avg", desc + f", {AVG_START}-latest average", url, rawf, "FRED file(s) missing", geography, unit, "", src, f"{num_id}/{den_id}")150 return None151 r = (num / den * scale).dropna()152 save_series(f"ratio_{num_id}_over_{den_id}", r)153 freq = "annual" if annual else "quarterly"154 d, v = latest(r)155 note = f"{fred_desc(num_id)} divided by {fred_desc(den_id)}. {extra_note}".strip()156 add(tid + "_latest", desc + ", latest observation", v, unit, period_label(d) if not annual else str(d.year), geography, freq, src, f"{num_id}/{den_id}", url, rawf, note)157 if annual:158 mv, mlab = mean_from(r)159 else:160 a = annual_mean(r)161 mv, mlab = mean_from(a)162 mlab += " (annual means of quarterly ratios)"163 add(tid + "_avg", desc + f", {AVG_START}-latest average", mv, unit, mlab.split(" mean")[0], geography, "annual average", src, f"{num_id}/{den_id}", url, rawf, note + f" Average: {mlab}.")164 # also latest full calendar-year average for quarterly series165 if not annual:166 a = annual_mean(r)167 if len(a):168 add(tid + "_lastfullyear", desc + ", latest full calendar year (mean of quarters)", float(a.iloc[-1]), unit, str(a.index[-1].year), geography, "annual average", src, f"{num_id}/{den_id}", url, rawf, note)169 return r170171172def level_targets(tid, desc, sid, unit, geography="United States", avg=True, transform=None, extra_note="", annual_from_hf=False):173 s = fred(sid)174 src = "FRED (" + fred_meta.get(sid, {}).get("source", "") + ")"175 if s is None:176 not_retrieved(tid + "_latest", desc + ", latest", fred_url(sid), fred_file(sid), "FRED file missing", geography, unit, "", src, sid)177 return None178 if transform is not None:179 s = transform(s)180 d, v = latest(s)181 note = (fred_desc(sid) + ". " + extra_note).strip()182 add(tid + "_latest", desc + ", latest observation", v, unit, period_label(d, fred_meta.get(sid, {}).get("frequency", "")), geography, fred_meta.get(sid, {}).get("frequency", ""), src, sid, fred_url(sid), fred_file(sid), note)183 if annual_from_hf:184 a = annual_mean(s)185 if len(a):186 add(tid + "_lastfullyear", desc + ", latest full calendar-year average", float(a.iloc[-1]), unit, str(a.index[-1].year), geography, "annual average", src, sid, fred_url(sid), fred_file(sid), note)187 if avg:188 mv, mlab = mean_from(a)189 add(tid + "_avg", desc + f", {AVG_START}-latest average of annual means", mv, unit, mlab.split(" mean")[0], geography, "annual average", src, sid, fred_url(sid), fred_file(sid), note + f" Average: {mlab}.")190 elif avg:191 mv, mlab = mean_from(s)192 add(tid + "_avg", desc + f", {AVG_START}-latest average", mv, unit, mlab.split(" mean")[0], geography, fred_meta.get(sid, {}).get("frequency", ""), src, sid, fred_url(sid), fred_file(sid), note + f" Average: {mlab}.")193 return s194195196# ============================================================================197# 1. Labor share198# ============================================================================199level_targets("T01_labor_share_BLS_nonfarm_index", "BLS nonfarm business sector labor share (index 2017=100)", "PRS85006173", "index 2017=100",200 extra_note="Index, not a level; 2017=100. Use for trends only.", annual_from_hf=True)201ratio_targets("T01_labor_share_COE_GDI", "National-accounts labor share: compensation of employees / gross domestic income", "COE", "GDI")202ratio_targets("T01_labor_share_COE_NDI", "National-accounts labor share: compensation of employees / net domestic income", "COE", "W256RC1Q027SBEA",203 extra_note="Net domestic income = GDI less consumption of fixed capital (BEA NIPA table 1.7.5 line 'Net domestic income').")204level_targets("T01_labor_share_PWT_US", "Penn World Table labour-compensation share of GDP, United States", "LABSHPUSA156NRUG", "ratio",205 extra_note="FRED serves PWT 11.0 (1950-2023). The DataverseNL file pwt1001.xlsx (PWT 10.01, to 2019) is kept as cross-check below.")206207# PWT 10.01 spreadsheet cross-check (labsh)208pwt_p = RAW / "pwt" / "pwt1001.xlsx"209if pwt_p.exists():210 pwt = pd.read_excel(pwt_p, sheet_name="Data", usecols=["countrycode", "year", "labsh"])211 for cc, geo, tid in (("USA", "United States", "T01_labor_share_PWT1001_US"), ("CAN", "Canada", "T12_labor_share_PWT1001_CA")):212 s = pwt[pwt.countrycode == cc].dropna(subset=["labsh"]).set_index("year")["labsh"]213 s.index = pd.to_datetime([f"{y}-01-01" for y in s.index])214 save_series(f"pwt1001_labsh_{cc}", s)215 d, v = latest(s)216 mv, mlab = mean_from(s)217 add(tid + "_latest", f"PWT 10.01 labsh (share of labour compensation in GDP at current national prices), {geo}, latest", v, "ratio", str(d.year), geo, "annual",218 "Penn World Table 10.01 (Feenstra, Inklaar & Timmer 2015), DataverseNL doi:10.34894/QT5BCC", "labsh", "https://dataverse.nl/api/access/datafile/354095", "data/raw/pwt/pwt1001.xlsx", "CC BY 4.0.")219 add(tid + "_avg", f"PWT 10.01 labsh, {geo}, {AVG_START}-latest average", mv, "ratio", mlab.split(" mean")[0], geo, "annual",220 "Penn World Table 10.01", "labsh", "https://dataverse.nl/api/access/datafile/354095", "data/raw/pwt/pwt1001.xlsx", mlab)221else:222 not_retrieved("T01_labor_share_PWT1001_US_latest", "PWT 10.01 labsh USA", "https://dataverse.nl/api/access/datafile/354095", "data/raw/pwt/pwt1001.xlsx", "file missing")223224# ============================================================================225# 2. Profits, dividends, GDP226# ============================================================================227ratio_targets("T02_profits_after_tax_noIVA_GDP", "Corporate profits after tax (without IVA and CCAdj) / GDP", "CP", "GDP",228 extra_note="FRED id CP is the series WITHOUT IVA/CCAdj (book profits).")229ratio_targets("T02_profits_after_tax_IVA_CCAdj_GDP", "Corporate profits after tax with IVA and CCAdj / GDP", "CPATAX", "GDP",230 extra_note="FRED id CPATAX is the series WITH IVA and CCAdj (economic profits).")231ratio_targets("T02_net_dividends_GDP", "Net corporate dividends (NIPA, after tax with IVA & CCAdj) / GDP", "DIVIDEND", "GDP")232ratio_targets("T02_net_dividends_annual_GDP", "Net corporate dividend payments (annual NIPA) / GDP (annual)", "B056RC1A027NBEA", "GDPA", annual=True)233level_targets("T02_GDP_annual_nominal", "Nominal GDP, annual", "GDPA", "billions of current US dollars", avg=False)234level_targets("T02_GDP_quarterly_nominal", "Nominal GDP, quarterly SAAR", "GDP", "billions of current US dollars, SAAR", avg=False)235level_targets("T02_GDP_real", "Real GDP, quarterly SAAR", "GDPC1", "billions of chained 2017 dollars, SAAR", avg=False)236# real GDP growth 1990-latest (annual means, log growth)237g = fred("GDPC1")238if g is not None:239 a = annual_mean(g)240 a = a[a.index.year >= AVG_START]241 gr = math.log(a.iloc[-1] / a.iloc[0]) / (a.index[-1].year - a.index[0].year)242 add("T02_real_GDP_growth_avg", f"Average real GDP growth, {a.index[0].year}-{a.index[-1].year} (log, annual means)", gr, "per year",243 f"{a.index[0].year}-{a.index[-1].year}", "United States", "annual", "FRED (BEA)", "GDPC1", fred_url("GDPC1"), fred_file("GDPC1"), fred_desc("GDPC1"))244245# ============================================================================246# 3. Household balance sheet247# ============================================================================248ratio_targets("T03_household_net_worth_GDP", "Household & nonprofit net worth / GDP", "TNWBSHNO", "GDP", unit="ratio (years of GDP)", scale=1e-3,249 extra_note="TNWBSHNO in millions of dollars, GDP in billions SAAR; scaled by 1e-3. Quarterly end-of-period net worth over same-quarter annualized GDP.")250ratio_targets("T03_household_real_estate_GDP", "Household & nonprofit real estate at market value / GDP", "HNOREMQ027S", "GDP", unit="ratio (years of GDP)", scale=1e-3,251 extra_note="Millions over billions SAAR; scaled by 1e-3.")252ratio_targets("T03_owners_equity_real_estate_GDP", "Households' owners' equity in real estate / GDP", "OEHRENWBSHNO", "GDP", unit="ratio (years of GDP)", scale=1e-3,253 extra_note="Millions over billions SAAR; scaled by 1e-3.")254level_targets("T03_owners_equity_share_of_real_estate", "Owners' equity as % of household real estate", "HOEREPHRE", "percent", annual_from_hf=True)255ratio_targets("T03_real_estate_share_of_net_worth", "Household real estate / household net worth", "HNOREMQ027S", "TNWBSHNO", unit="ratio")256257# ============================================================================258# 4. Wealth concentration (DFA)259# ============================================================================260dfa_p = RAW / "fed_dfa" / "unzipped" / "dfa-networth-shares.csv"261dfa_url = "https://www.federalreserve.gov/releases/z1/dataviz/download/zips/dfa.zip"262if dfa_p.exists():263 dfa = pd.read_csv(dfa_p)264 w = dfa.pivot(index="Date", columns="Category", values="Net worth")265 w.index = pd.to_datetime([f"{d[:4]}-{int(d[-1]) * 3 - 2:02d}-01" for d in w.index])266 w = w.sort_index()267 top1 = w["TopPt1"] + w["RemainingTop1"]268 top10 = top1 + w["Next9"]269 bot50 = w["Bottom50"]270 mid = w["Next40"]271 for name, s, desc in (("top1", top1, "Share of household net worth held by the top 1% (TopPt1 + RemainingTop1)"),272 ("top10", top10, "Share of household net worth held by the top 10% (top 1% + Next9)"),273 ("top0p1", w["TopPt1"], "Share of household net worth held by the top 0.1%"),274 ("next40", mid, "Share of household net worth held by the 50th-90th percentiles"),275 ("bottom50", bot50, "Share of household net worth held by the bottom 50%")):276 save_series(f"dfa_networth_share_{name}", s)277 d, v = latest(s)278 add(f"T04_wealth_share_{name}_latest", desc + ", latest quarter", v, "percent", period_label(d), "United States", "quarterly",279 "Federal Reserve Board, Distributional Financial Accounts (dfa-networth-shares.csv, release dated 2026-09-15 in zip)", f"dfa-networth-shares:Net worth:{name}",280 dfa_url, "data/raw/fed_dfa/unzipped/dfa-networth-shares.csv", "Percent of total household net worth; DFA combine SCF and Financial Accounts. Public domain.")281 mv, mlab = mean_from(annual_mean(s))282 add(f"T04_wealth_share_{name}_avg", desc + f", {AVG_START}-latest average", mv, "percent", mlab.split(" mean")[0], "United States", "annual average",283 "Federal Reserve Board, Distributional Financial Accounts", f"dfa-networth-shares:Net worth:{name}", dfa_url, "data/raw/fed_dfa/unzipped/dfa-networth-shares.csv", mlab)284 # corporate-equity concentration (useful for the ownership block)285 eq = dfa.pivot(index="Date", columns="Category", values="Corporate equities and mutual fund shares")286 eq.index = w.index287 eq_top10 = (eq["TopPt1"] + eq["RemainingTop1"] + eq["Next9"]).sort_index()288 save_series("dfa_corporate_equity_share_top10", eq_top10)289 d, v = latest(eq_top10)290 add("T04_corporate_equity_share_top10_latest", "Share of directly/indirectly held corporate equities & mutual fund shares held by the top 10%", v, "percent", period_label(d),291 "United States", "quarterly", "Federal Reserve Board, Distributional Financial Accounts", "dfa-networth-shares:Corporate equities and mutual fund shares:top10",292 dfa_url, "data/raw/fed_dfa/unzipped/dfa-networth-shares.csv", "Excludes equities held through DC pensions (DFA definition).")293else:294 for n in ("top1", "top10", "bottom50"):295 not_retrieved(f"T04_wealth_share_{n}_latest", f"DFA net worth share {n}", dfa_url, "data/raw/fed_dfa/unzipped/dfa-networth-shares.csv", "dfa.zip missing")296# FRED cross-checks of the same DFA numbers297for sid, tid, desc in (("WFRBST01134", "T04_wealth_share_top1_FRED", "Share of net worth held by the top 1% (FRED mirror of DFA)"),298 ("WFRBSN09161", "T04_wealth_share_p90_99_FRED", "Share of net worth held by the 90th-99th percentiles (FRED mirror of DFA)"),299 ("WFRBSN40188", "T04_wealth_share_p50_90_FRED", "Share of net worth held by the 50th-90th percentiles (FRED mirror of DFA)"),300 ("WFRBSB50215", "T04_wealth_share_bottom50_FRED", "Share of net worth held by the bottom 50% (FRED mirror of DFA)")):301 level_targets(tid, desc, sid, "percent", avg=False)302303# ============================================================================304# 5. Investment305# ============================================================================306ratio_targets("T05_GPDI_GDP", "Gross private domestic investment / GDP", "GPDI", "GDP")307ratio_targets("T05_PNFI_GDP", "Private nonresidential fixed investment / GDP", "PNFI", "GDP")308ratio_targets("T05_info_processing_equipment_GDP", "Nonresidential investment in information processing equipment / GDP", "Y034RC1Q027SBEA", "GDP")309ratio_targets("T05_software_GDP", "Nonresidential investment in software (IPP) / GDP", "B985RC1Q027SBEA", "GDP")310ratio_targets("T05_IPE_and_software_GDP", "Private fixed investment in information processing equipment and software / GDP", "A679RC1Q027SBEA", "GDP")311ratio_targets("T05_IPP_GDP", "Nonresidential investment in intellectual property products / GDP", "Y001RC1Q027SBEA", "GDP",312 extra_note="The brief suggested Y001RC1Q027SBEA as 'total nonresidential'; FRED metadata shows it is the IPP component. PNFI is the total.")313314# ============================================================================315# 6. Capital stock and depreciation (BEA Fixed Assets, current cost)316# ============================================================================317ratio_targets("T06_net_stock_private_fixed_assets_GDP", "Current-cost net stock of private fixed assets (total) / GDP", "K1PTOTL1ES000", "GDPA", unit="ratio (years of GDP)", scale=1e-3, annual=True,318 extra_note="End-of-year current-cost net stock (millions) over annual nominal GDP (billions); scaled 1e-3.")319ratio_targets("T06_net_stock_private_nonres_fixed_assets_GDP", "Current-cost net stock of private nonresidential fixed assets / GDP", "K1NTOTL1ES000", "GDPA", unit="ratio (years of GDP)", scale=1e-3, annual=True)320ratio_targets("T06_depreciation_rate_private_total", "Implied depreciation rate: current-cost depreciation / net stock, private fixed assets (total)", "M1PTOTL1ES000", "K1PTOTL1ES000", unit="per year", annual=True,321 extra_note="Both current-cost, BEA Fixed Assets tables 1.1 and 1.3; depreciation during year t over end-of-year-t net stock (not lagged stock).")322ratio_targets("T06_depreciation_rate_private_nonres", "Implied depreciation rate: private nonresidential fixed assets", "M1NTOTL1ES000", "K1NTOTL1ES000", unit="per year", annual=True,323 extra_note="Current-cost depreciation over end-of-year current-cost net stock, nonresidential.")324ratio_targets("T06_CFC_GDP", "Consumption of fixed capital (NIPA, total economy) / GDP", "COFC", "GDP")325ratio_targets("T06_CFC_private_GDP", "Consumption of fixed capital: private / GDP (annual)", "A024RC1A027NBEA", "GDPA", annual=True)326ratio_targets("T06_CFC_private_over_net_stock", "NIPA private CFC / current-cost net stock of private fixed assets", "A024RC1A027NBEA", "K1PTOTL1ES000", unit="per year", scale=1e3, annual=True,327 extra_note="Billions over millions, scaled 1e3. Alternative depreciation-rate measure.")328329# ============================================================================330# 7. Real interest rates331# ============================================================================332level_targets("T07_TIPS10_daily", "10-year TIPS constant-maturity real yield", "DFII10", "percent per year", annual_from_hf=True,333 extra_note="Average over 2003-latest (TIPS data start 2003), not 1990.")334level_targets("T07_TIPS10_monthly", "10-year TIPS constant-maturity real yield (monthly)", "FII10", "percent per year", annual_from_hf=True)335level_targets("T07_real_rate_10y_Cleveland", "10-year real interest rate (Cleveland Fed model: nominal yield minus expected inflation)", "REAINTRATREARAT10Y", "percent per year", annual_from_hf=True)336level_targets("T07_nominal_10y", "10-year Treasury constant-maturity nominal yield (monthly)", "GS10", "percent per year", annual_from_hf=True)337# Treasury.gov cross-check (daily par real yield curve 10 YR)338tr_files = sorted((RAW / "treasury").glob("daily_treasury_real_yield_curve_*.csv"))339if tr_files:340 frames = []341 for f in tr_files:342 d = pd.read_csv(f)343 d.columns = [c.strip() for c in d.columns]344 frames.append(d[["Date", "10 YR"]])345 tr = pd.concat(frames)346 tr["Date"] = pd.to_datetime(tr["Date"], format="%m/%d/%Y")347 s = pd.Series(pd.to_numeric(tr["10 YR"], errors="coerce").values, index=tr["Date"]).sort_index().dropna()348 save_series("treasury_real_yield_10y_daily", s)349 a = annual_mean(s)350 add("T07_TIPS10_treasury_lastfullyear", "10-year par real yield (US Treasury daily real yield curve), latest full calendar-year average", float(a.iloc[-1]), "percent per year", str(a.index[-1].year),351 "United States", "annual average", "US Department of the Treasury, Daily Treasury Par Real Yield Curve Rates", "10 YR",352 "https://home.treasury.gov/resource-center/data-chart-center/interest-rates/daily-treasury-rates.csv/<year>/all?type=daily_treasury_real_yield_curve&field_tdr_date_value=<year>&page&_format=csv",353 "data/raw/treasury/daily_treasury_real_yield_curve_<year>.csv", "Cross-check of DFII10; public domain.")354 mv, mlab = mean_from(a, 2003)355 add("T07_TIPS10_treasury_avg", "10-year par real yield (Treasury), 2003-latest average of annual means", mv, "percent per year", mlab.split(" mean")[0], "United States", "annual average",356 "US Department of the Treasury", "10 YR", "https://home.treasury.gov/resource-center/data-chart-center/interest-rates/", "data/raw/treasury/", mlab)357# Cleveland Fed workbook cross-check358cf_p = RAW / "clevelandfed" / "inflation-expectations.xlsx"359if cf_p.exists():360 cf = pd.read_excel(cf_p, sheet_name="Real Interest Rate")361 cf.columns = [str(c).strip() for c in cf.columns]362 s = pd.Series(pd.to_numeric(cf["Real Rate 10-year"], errors="coerce").values * 100, index=pd.to_datetime(cf["Model Output Date"])).dropna()363 save_series("clevelandfed_real_rate_10y_monthly", s)364 d, v = latest(s)365 add("T07_real_rate_10y_Cleveland_xlsx_latest", "10-year real rate from the Cleveland Fed inflation-expectations workbook (decimal x100)", v, "percent per year", d.strftime("%Y-%m"),366 "United States", "monthly", "Federal Reserve Bank of Cleveland, Inflation Expectations model", "Real Rate 10-year",367 "https://www.clevelandfed.org/-/media/files/webcharts/inflationexpectations/inflation-expectations.xlsx", "data/raw/clevelandfed/inflation-expectations.xlsx",368 "Cross-check of REAINTRATREARAT10Y. Model of Haubrich, Pennacchi & Ritchken (2012).")369 ei = pd.read_excel(cf_p, sheet_name="Ten-year Expected Chart")370 ei.columns = [str(c).strip() for c in ei.columns]371 s2 = pd.Series(pd.to_numeric(ei["10 year Expected Inflation"], errors="coerce").values, index=pd.to_datetime(ei["Model Output Date"])).dropna()372 save_series("clevelandfed_expected_inflation_10y_monthly", s2)373 d, v = latest(s2)374 add("T07_expected_inflation_10y_Cleveland_latest", "10-year expected inflation (Cleveland Fed model)", v, "percent per year", d.strftime("%Y-%m"), "United States", "monthly",375 "Federal Reserve Bank of Cleveland", "10 year Expected Inflation", "https://www.clevelandfed.org/-/media/files/webcharts/inflationexpectations/inflation-expectations.xlsx",376 "data/raw/clevelandfed/inflation-expectations.xlsx", "")377378# ============================================================================379# 8. Equity valuation (Shiller)380# ============================================================================381sh_p = RAW / "shiller" / "ie_data.xls"382sh_url = "http://www.econ.yale.edu/~shiller/data/ie_data.xls"383if sh_p.exists():384 raw = pd.read_excel(sh_p, sheet_name="Data", header=None)385 hdr_row = raw.index[raw.iloc[:, 0].astype(str).str.strip() == "Date"][0]386 df = raw.iloc[hdr_row + 1:, :13].copy()387 df.columns = ["Date", "P", "D", "E", "CPI", "DateFraction", "GS10", "RealP", "RealD", "RealTRPrice", "RealE", "RealTRScaledE", "CAPE"]388 df = df[pd.to_numeric(df["DateFraction"], errors="coerce").notna()]389 for c in df.columns:390 df[c] = pd.to_numeric(df[c], errors="coerce")391 # month from the date fraction (year + (m-0.5)/12)392 yr = df["DateFraction"].astype(float)393 y = np.floor(yr).astype(int)394 m = np.clip(np.round((yr - y) * 12 + 0.5).astype(int), 1, 12)395 idx = pd.to_datetime([f"{a}-{b:02d}-01" for a, b in zip(y, m)])396 df.index = idx397 cape = df["CAPE"].dropna()398 dp = (df["D"] / df["P"]).dropna()399 ep = (df["E"] / df["P"]).dropna()400 for name, s in (("cape", cape), ("dividend_yield", dp), ("earnings_yield", ep), ("sp_price", df["P"]), ("sp_dividend", df["D"]), ("sp_earnings", df["E"]), ("cpi", df["CPI"]), ("gs10", df["GS10"])):401 save_series(f"shiller_{name}", s)402 vint = f"file ends {df.index[-1].strftime('%Y-%m')} (the Yale copy retrieved on {RETRIEVED} is not more recent)"403 for tid, s, desc, unit in (("T08_CAPE", cape, "Shiller cyclically adjusted P/E (CAPE, P/E10)", "ratio"),404 ("T08_dividend_yield", dp, "S&P composite dividend yield D/P (trailing 12-month dividends / price)", "ratio"),405 ("T08_earnings_yield", ep, "S&P composite earnings yield E/P (trailing 12-month earnings / price)", "ratio")):406 d, v = latest(s)407 add(tid + "_latest", desc + ", latest observation in file", v, unit, d.strftime("%Y-%m"), "United States", "monthly",408 "Robert J. Shiller, online data for 'Irrational Exuberance' (ie_data.xls)", "ie_data.xls Data sheet", sh_url, "data/raw/shiller/ie_data.xls",409 f"VINTAGE CAVEAT: {vint}. D and E are 12-month trailing sums interpolated monthly by Shiller; P is the monthly average of daily closes.")410 a = annual_mean(s)411 mv, mlab = mean_from(a)412 add(tid + "_avg", desc + f", {AVG_START}-latest average of annual means", mv, unit, mlab.split(" mean")[0], "United States", "annual average",413 "Robert J. Shiller, online data (ie_data.xls)", "ie_data.xls Data sheet", sh_url, "data/raw/shiller/ie_data.xls", mlab + f". {vint}.")414 # long-run (1881-latest) CAPE mean as context415 add("T08_CAPE_longrun_avg", "Shiller CAPE, full-sample mean of annual means", float(annual_mean(cape).mean()), "ratio",416 f"{annual_mean(cape).index[0].year}-{annual_mean(cape).index[-1].year}", "United States", "annual average", "Robert J. Shiller, online data (ie_data.xls)", "ie_data.xls Data sheet", sh_url, "data/raw/shiller/ie_data.xls", vint)417else:418 for t in ("T08_CAPE_latest", "T08_dividend_yield_latest", "T08_earnings_yield_latest"):419 not_retrieved(t, t, sh_url, "data/raw/shiller/ie_data.xls", "ie_data.xls missing")420421# ============================================================================422# 9. Housing423# ============================================================================424f31_p = RAW / "shiller" / "Fig3-1.xls"425f31_url = "http://www.econ.yale.edu/~shiller/data/Fig3-1.xls"426if f31_p.exists():427 raw = pd.read_excel(f31_p, header=None)428 hdr_row = raw.index[raw.iloc[:, 0].astype(str).str.strip() == "Date"][0]429 dates = pd.to_numeric(raw.iloc[hdr_row + 1:, 0], errors="coerce")430 vals = pd.to_numeric(raw.iloc[hdr_row + 1:, 1], errors="coerce")431 ok = dates.notna() & vals.notna()432 yrs = dates[ok].astype(float).values433 yy = np.floor(yrs + 1e-9).astype(int)434 mm = np.clip(np.round((yrs - yy) * 12 + 1).astype(int), 1, 12) # annual rows (integer years) -> January435 idx = pd.to_datetime([f"{a}-{b:02d}-01" for a, b in zip(yy, mm)])436 rhpi = pd.Series(vals[ok].values, index=idx).sort_index()437 rhpi = rhpi[~rhpi.index.duplicated()]438 save_series("shiller_real_home_price_index", rhpi)439 d, v = latest(rhpi)440 add("T09_real_home_price_index_latest", "Shiller US real home price index (1890=100), latest observation in file", v, "index 1890=100", d.strftime("%Y-%m"), "United States", "monthly (annual before 1953)",441 "Robert J. Shiller, Irrational Exuberance Fig. 3.1 data (Fig3-1.xls)", "Fig3-1.xls col B", f31_url, "data/raw/shiller/Fig3-1.xls",442 f"VINTAGE CAVEAT: file ends {d.strftime('%Y-%m')}. Nominal index = Case-Shiller national (S&P/CoreLogic) spliced to Grebler/PHCPI; deflated by CPI by Shiller.")443 a = annual_mean(rhpi, full_years_only=False)444 for (y0, tid) in ((AVG_START, "T09_real_home_price_growth_1990"), (1890, "T09_real_home_price_growth_1890")):445 aa = a[a.index.year >= y0]446 gr = math.log(aa.iloc[-1] / aa.iloc[0]) / (aa.index[-1].year - aa.index[0].year)447 add(tid, f"Average real home price growth (log, annual means), {aa.index[0].year}-{aa.index[-1].year}", gr, "per year", f"{aa.index[0].year}-{aa.index[-1].year}", "United States", "annual",448 "Robert J. Shiller, Fig3-1.xls", "Fig3-1.xls col B", f31_url, "data/raw/shiller/Fig3-1.xls", "Computed from annual means of the real index.")449else:450 not_retrieved("T09_real_home_price_index_latest", "Shiller real home price index", f31_url, "data/raw/shiller/Fig3-1.xls", "file missing")451452# Rent/price proxy: CPI rent of primary residence over Case-Shiller national index, both rebased to Jan-2000 = 100453hp, rent, cpi = fred("CSUSHPISA"), fred("CUSR0000SEHA"), fred("CPIAUCSL")454if hp is not None and rent is not None:455 base = pd.Timestamp("2000-01-01")456 rp = (rent / rent[base]) / (hp / hp[base]) * 100457 rp = rp.dropna()458 save_series("rent_price_proxy_index_2000eq100", rp)459 d, v = latest(rp)460 note = ("PROXY ONLY: index of (CPI rent of primary residence) / (Case-Shiller national HPI), both rebased to Jan 2000 = 100. "461 "It tracks the change in the rent-price ratio relative to 2000, not its level; no key-less level series of the rent-price ratio was available.")462 add("T09_rent_price_ratio_proxy_latest", "Rent-price ratio proxy index (2000=100), latest", v, "index Jan-2000=100", d.strftime("%Y-%m"), "United States", "monthly",463 "FRED (BLS CUSR0000SEHA; S&P Dow Jones Indices CSUSHPISA)", "CUSR0000SEHA/CSUSHPISA", fred_url("CUSR0000SEHA") + " ; " + fred_url("CSUSHPISA"),464 fred_file("CUSR0000SEHA") + " ; " + fred_file("CSUSHPISA"), note)465 a = annual_mean(rp)466 mv, mlab = mean_from(a)467 add("T09_rent_price_ratio_proxy_avg", f"Rent-price ratio proxy index (2000=100), {AVG_START}-latest average", mv, "index Jan-2000=100", mlab.split(" mean")[0], "United States", "annual average",468 "FRED (BLS; S&P Dow Jones Indices)", "CUSR0000SEHA/CSUSHPISA", fred_url("CUSR0000SEHA") + " ; " + fred_url("CSUSHPISA"), fred_file("CUSR0000SEHA") + " ; " + fred_file("CSUSHPISA"), note + " " + mlab)469 # rent inflation vs house price inflation 1990-latest470 for sid, s, tid, desc in (("CUSR0000SEHA", rent, "T09_rent_cpi_growth_1990", "Average growth of CPI rent of primary residence"),471 ("CSUSHPISA", hp, "T09_case_shiller_growth_1990", "Average growth of Case-Shiller national home price index (nominal)")):472 a = annual_mean(s)473 a = a[a.index.year >= AVG_START]474 gr = math.log(a.iloc[-1] / a.iloc[0]) / (a.index[-1].year - a.index[0].year)475 add(tid, f"{desc} (log, annual means), {a.index[0].year}-{a.index[-1].year}", gr, "per year", f"{a.index[0].year}-{a.index[-1].year}", "United States", "annual",476 "FRED (" + fred_meta.get(sid, {}).get("source", "") + ")", sid, fred_url(sid), fred_file(sid), fred_desc(sid))477 if cpi is not None:478 real_hp = (hp / cpi).dropna()479 a = annual_mean(real_hp)480 a = a[a.index.year >= AVG_START]481 gr = math.log(a.iloc[-1] / a.iloc[0]) / (a.index[-1].year - a.index[0].year)482 add("T09_case_shiller_real_growth_1990", f"Average real growth of Case-Shiller national HPI deflated by CPI-U (log, annual means), {a.index[0].year}-{a.index[-1].year}", gr, "per year",483 f"{a.index[0].year}-{a.index[-1].year}", "United States", "annual", "FRED (S&P Dow Jones Indices; BLS)", "CSUSHPISA/CPIAUCSL", fred_url("CSUSHPISA") + " ; " + fred_url("CPIAUCSL"),484 fred_file("CSUSHPISA") + " ; " + fred_file("CPIAUCSL"), "")485else:486 not_retrieved("T09_rent_price_ratio_proxy_latest", "Rent-price proxy", fred_url("CUSR0000SEHA"), fred_file("CUSR0000SEHA"), "FRED files missing")487488# ============================================================================489# 10. Energy: data centers (LBNL) and total US electricity (EIA)490# ============================================================================491lb_p = RAW / "lbnl" / "lbnl-2024-united-states-data-center-energy-usage-report.pdf"492lb_entry = manifest.get("lbnl:report_pdf", {})493lb_url = lb_entry.get("url", "https://eta-publications.lbl.gov/sites/default/files/2024-12/lbnl-2024-united-states-data-center-energy-usage-report.pdf")494lb_src = "Shehabi, A., Newkirk, A., Smith, S.J., et al. (2024), 2024 United States Data Center Energy Usage Report, LBNL-2001637, Lawrence Berkeley National Laboratory (DOI 10.71468/P1WC7Q)"495if lb_p.exists():496 import fitz # PyMuPDF497 doc = fitz.open(lb_p)498 pages = {i + 1: doc[i].get_text() for i in range(doc.page_count)}499500 def find(pattern, want_pages):501 for p in want_pages:502 t = " ".join(pages[p].split())503 m = re.search(pattern, t)504 if m:505 return m, p506 return None, None507508 def printed_page(p):509 m = re.search(r"Energy Usage Report (\d+)", pages[p].split("\n")[0] + " " + pages[p][:120])510 return m.group(1) if m else "?"511512 m, p = find(r"reaching (\d+) TWh by 2023, representing ([\d.]+)% of total U\.S\. electricity consumption", [7, 53])513 if m:514 add("T10_datacenter_electricity_2023_TWh", "US data center electricity use, 2023", float(m.group(1)), "TWh", "2023", "United States", "annual", lb_src, "text",515 lb_url, "data/raw/lbnl/lbnl-2024-united-states-data-center-energy-usage-report.pdf",516 f"Quote (PDF page {p}, printed page {printed_page(p)}): '...{m.group(0)}'. Includes servers, storage, network and infrastructure (PUE), bottom-up model.")517 add("T10_datacenter_share_of_US_electricity_2023", "US data center electricity use as share of total US electricity consumption, 2023", float(m.group(2)), "percent", "2023", "United States", "annual",518 lb_src, "text", lb_url, "data/raw/lbnl/lbnl-2024-united-states-data-center-energy-usage-report.pdf", f"Quote (PDF page {p}, printed page {printed_page(p)}): '{m.group(0)}'.")519 m, p = find(r"by 2018 data centers consumed about (\d+) TWh, representing ([\d.]+)% of total (?:annual )?U\.S\. electricity consumption", [7, 53])520 if m:521 add("T10_datacenter_electricity_2018_TWh", "US data center electricity use, 2018", float(m.group(1)), "TWh", "2018", "United States", "annual", lb_src, "text", lb_url,522 "data/raw/lbnl/lbnl-2024-united-states-data-center-energy-usage-report.pdf", f"Quote (PDF page {p}, printed page {printed_page(p)}): '{m.group(0)}'.")523 add("T10_datacenter_share_of_US_electricity_2018", "US data center share of total US electricity consumption, 2018", float(m.group(2)), "percent", "2018", "United States", "annual", lb_src, "text", lb_url,524 "data/raw/lbnl/lbnl-2024-united-states-data-center-energy-usage-report.pdf", f"Quote (PDF page {p}, printed page {printed_page(p)}).")525 m, p = find(r"low and high end (?:of roughly|representing about) (\d+) and (\d+) TWh in 2028", [8, 53])526 if m:527 add("T10_datacenter_electricity_2028_low_TWh", "US data center electricity use, 2028 projection, low scenario", float(m.group(1)), "TWh", "2028", "United States", "annual (projection)", lb_src, "text", lb_url,528 "data/raw/lbnl/lbnl-2024-united-states-data-center-energy-usage-report.pdf", f"Quote (PDF page {p}, printed page {printed_page(p)}): '{m.group(0)}'. Scenario range, not a forecast.")529 add("T10_datacenter_electricity_2028_high_TWh", "US data center electricity use, 2028 projection, high scenario", float(m.group(2)), "TWh", "2028", "United States", "annual (projection)", lb_src, "text", lb_url,530 "data/raw/lbnl/lbnl-2024-united-states-data-center-energy-usage-report.pdf", f"Quote (PDF page {p}, printed page {printed_page(p)}).")531 m, p = find(r"([\d.]+)% to ([\d.]+)% of total U\.S\. electricity consumption", [8, 53])532 if m:533 add("T10_datacenter_share_2028_low", "US data center share of total US electricity consumption, 2028 projection, low", float(m.group(1)), "percent", "2028", "United States", "annual (projection)", lb_src, "text", lb_url,534 "data/raw/lbnl/lbnl-2024-united-states-data-center-energy-usage-report.pdf", f"Quote (PDF page {p}, printed page {printed_page(p)}): '{m.group(0)}'.")535 add("T10_datacenter_share_2028_high", "US data center share of total US electricity consumption, 2028 projection, high", float(m.group(2)), "percent", "2028", "United States", "annual (projection)", lb_src, "text", lb_url,536 "data/raw/lbnl/lbnl-2024-united-states-data-center-energy-usage-report.pdf", f"Quote (PDF page {p}, printed page {printed_page(p)}).")537else:538 not_retrieved("T10_datacenter_electricity_2023_TWh", "US data center electricity use 2023", lb_url, "data/raw/lbnl/", "PDF missing (LBNL host returned HTTP 202 challenge page; eScholarship mirror not retrieved)")539540eia_p = RAW / "eia" / "epa_01_01.xlsx"541eia_url = "https://www.eia.gov/electricity/annual/xls/epa_01_01.xlsx"542if eia_p.exists():543 t = pd.read_excel(eia_p, header=None)544 hdr = t.index[t.iloc[:, 0].astype(str).str.strip() == "Fuel"][0]545 years = [str(t.iloc[hdr, 2]).replace("Year ", ""), str(t.iloc[hdr, 3]).replace("Year ", "")]546 row = t.index[(t.iloc[:, 0].astype(str).str.strip() == "All Energy Sources") & (t.iloc[:, 1].astype(str).str.contains("Utility Scale"))][0]547 for j, y in enumerate(years):548 add(f"T10_US_net_generation_{y}_TWh", f"US net electricity generation, all sectors, utility-scale facilities, {y}", float(t.iloc[row, 2 + j]) / 1e3, "TWh", y, "United States", "annual",549 "EIA, Electric Power Annual / Electric Power Monthly Table 1.1 (Total Electric Power Industry Summary Statistics)", "All Energy Sources, Utility Scale Facilities, Net Generation (thousand MWh)", eia_url,550 "data/raw/eia/epa_01_01.xlsx", "Thousand megawatthours converted to TWh (/1000). Public domain. Excludes small-scale solar (reported separately in the table).")551 row2 = t.index[(t.iloc[:, 0].astype(str).str.strip() == "All Sectors")]552 if len(row2):553 r2 = row2[0]554 hdr2 = t.index[(t.iloc[:, 0].astype(str).str.strip() == "Sector")][0]555 years2 = [str(t.iloc[hdr2, 1]).replace("Year ", ""), str(t.iloc[hdr2, 2]).replace("Year ", "")]556 for j, y in enumerate(years2):557 v = float(t.iloc[r2, 1 + j]) / 1e3558 add(f"T10_US_retail_sales_{y}_TWh", f"US retail sales of electricity to ultimate customers, all sectors, {y}", v, "TWh", y, "United States", "annual",559 "EIA, Electric Power Annual Table 1.1 (Sales of Electricity to Ultimate Customers, thousand MWh)", "All Sectors retail sales", eia_url, "data/raw/eia/epa_01_01.xlsx",560 "Thousand MWh / 1000. Retail sales approximate end-use consumption (exclude direct use and losses).")561else:562 not_retrieved("T10_US_net_generation_latest_TWh", "US net generation", eia_url, "data/raw/eia/epa_01_01.xlsx", "file missing")563564# ============================================================================565# 11. Compute price-performance (Epoch AI hardware dataset)566# ============================================================================567ep_p = RAW / "epoch" / "ml_hardware.csv"568ep_url = "https://epoch.ai/data/ml_hardware.csv"569if ep_p.exists():570 e = pd.read_csv(ep_p)571 e["Release date"] = pd.to_datetime(e["Release date"], errors="coerce")572 e["pp"] = pd.to_numeric(e["Price-performance"], errors="coerce")573 e["price"] = pd.to_numeric(e["Release price (USD)"], errors="coerce")574 e["maxperf"] = pd.to_numeric(e["Max performance"], errors="coerce")575 sub = e.dropna(subset=["Release date", "pp"]).copy()576 sub["t"] = sub["Release date"].dt.year + (sub["Release date"].dt.dayofyear - 1) / 365.25577 sub["log2pp"] = np.log2(sub["pp"])578 sub[["Hardware name", "Manufacturer", "Type", "Release date", "price", "maxperf", "pp", "t", "log2pp"]].sort_values("Release date").to_csv(PROC / "epoch_price_performance_fit.csv", index=False)579580 def fit(df):581 X = np.vstack([np.ones(len(df)), df["t"].values]).T582 b, res, *_ = np.linalg.lstsq(X, df["log2pp"].values, rcond=None)583 slope = b[1]584 yhat = X @ b585 r2 = 1 - ((df["log2pp"] - yhat) ** 2).sum() / ((df["log2pp"] - df["log2pp"].mean()) ** 2).sum()586 return slope, 1 / slope, r2587588 for tid, df_, lab in (("T11_compute_price_performance_doubling_all", sub, "all hardware with a price and price-performance entry"),589 ("T11_compute_price_performance_doubling_2015plus", sub[sub["t"] >= 2015], "hardware released 2015 or later")):590 if len(df_) >= 5:591 slope, dbl, r2 = fit(df_)592 add(tid, f"Doubling time of FLOP/s per US dollar (Epoch 'Price-performance' column), OLS of log2(price-performance) on release year, {lab}", dbl, "years per doubling",593 f"{int(df_['t'].min())}-{int(df_['t'].max())}", "global (hardware market)", "cross-section of hardware releases",594 "Epoch AI, 'Machine Learning Hardware' dataset (https://epoch.ai/data/machine-learning-hardware), CC BY 4.0", "Price-performance; Release date", ep_url, "data/raw/epoch/ml_hardware.csv",595 f"n={len(df_)}, slope={slope:.3f} doublings/yr (x{2**slope:.2f}/yr), R2={r2:.2f}. Price-performance = Epoch's 'Max performance' FLOP/s (best available precision) over release price; our own OLS, not Epoch's published estimate. Last modified field max: {e['Last modified'].max() if 'Last modified' in e else '?'}.")596 else:597 not_retrieved(tid, "Epoch price-performance doubling time", ep_url, "data/raw/epoch/ml_hardware.csv", "too few rows")598 add("T11_compute_price_performance_caveat", "Caveat for T11: Epoch's 'Max performance' uses the highest-throughput precision available for each chip (FP32 for pre-2017 GPUs, tensor FP16/BF16 or lower for recent accelerators), so part of the trend reflects precision changes, not only hardware progress", float("nan"), "", "", "global", "", "Epoch AI", "Max performance", ep_url, "data/raw/epoch/ml_hardware.csv", "Informational row; no value.", status="INFO")599 add("T11_compute_price_performance_n", "Number of hardware entries with release price and price-performance", float(len(sub)), "count", f"{int(sub['t'].min())}-{int(sub['t'].max())}", "global", "", "Epoch AI", "Price-performance", ep_url, "data/raw/epoch/ml_hardware.csv", "")600else:601 not_retrieved("T11_compute_price_performance_doubling_all", "Epoch price-performance doubling time", ep_url, "data/raw/epoch/ml_hardware.csv", "file missing")602603# ============================================================================604# 12. Canada605# ============================================================================606level_targets("T12_labor_share_PWT_CA", "Penn World Table labour-compensation share of GDP, Canada", "LABSHPCAA156NRUG", "ratio", geography="Canada")607sc1 = RAW / "statcan" / "unzipped" / "36100103.csv"608sc2 = RAW / "statcan" / "unzipped" / "36100580.csv"609sc_lic = "Statistics Canada Open Licence (https://www.statcan.gc.ca/en/reference/licence)."610if sc1.exists():611 d = pd.read_csv(sc1, usecols=["REF_DATE", "GEO", "Seasonal adjustment", "Estimates", "VALUE", "SCALAR_FACTOR"], low_memory=False)612 d = d[(d["GEO"] == "Canada") & (d["Seasonal adjustment"] == "Seasonally adjusted at annual rates")]613 piv = d.pivot_table(index="REF_DATE", columns="Estimates", values="VALUE")614 piv.index = pd.to_datetime(piv.index + "-01")615 coe, gdp = piv["Compensation of employees"], piv["Gross domestic product at market prices"]616 ls = (coe / gdp).dropna()617 save_series("statcan_36100103_COE_over_GDP_quarterly", ls)618 save_series("statcan_36100103_GDP_market_prices_SAAR", gdp)619 d_, v = latest(ls)620 add("T12_labor_share_COE_GDP_CA_latest", "Canada: compensation of employees / GDP at market prices (income-based, SAAR)", v, "ratio", period_label(d_), "Canada", "quarterly",621 "Statistics Canada, Table 36-10-0103-01 (GDP income-based, quarterly)", "Compensation of employees / Gross domestic product at market prices", "https://www150.statcan.gc.ca/n1/tbl/csv/36100103-eng.zip",622 "data/raw/statcan/unzipped/36100103.csv", "Both in millions of dollars SAAR. " + sc_lic)623 mv, mlab = mean_from(annual_mean(ls))624 add("T12_labor_share_COE_GDP_CA_avg", f"Canada: compensation of employees / GDP, {AVG_START}-latest average", mv, "ratio", mlab.split(" mean")[0], "Canada", "annual average",625 "Statistics Canada, Table 36-10-0103-01", "COE/GDP", "https://www150.statcan.gc.ca/n1/tbl/csv/36100103-eng.zip", "data/raw/statcan/unzipped/36100103.csv", mlab + ". " + sc_lic)626 if sc2.exists():627 b = pd.read_csv(sc2, usecols=["REF_DATE", "Sectors", "Valuation", "Categories", "VALUE"], low_memory=False)628 b = b[(b["Sectors"] == "Households and non-profit institutions serving households") & (b["Valuation"] == "Market value") & (b["Categories"] == "Net worth")]629 nw = pd.Series(b["VALUE"].values, index=pd.to_datetime(b["REF_DATE"] + "-01")).sort_index()630 save_series("statcan_36100580_household_NPISH_net_worth_market_value", nw)631 r = (nw / gdp).dropna()632 save_series("statcan_household_net_worth_over_GDP", r)633 d_, v = latest(r)634 add("T12_household_net_worth_GDP_CA_latest", "Canada: households & NPISH net worth (market value) / GDP at market prices (SAAR)", v, "ratio (years of GDP)", period_label(d_), "Canada", "quarterly",635 "Statistics Canada, Tables 36-10-0580-01 (national balance sheet) and 36-10-0103-01 (GDP)", "Net worth (Households and NPISH, market value) / GDP", "https://www150.statcan.gc.ca/n1/tbl/csv/36100580-eng.zip ; https://www150.statcan.gc.ca/n1/tbl/csv/36100103-eng.zip",636 "data/raw/statcan/unzipped/36100580.csv ; data/raw/statcan/unzipped/36100103.csv", "Both in millions of CAD; end-of-quarter net worth over same-quarter annualized GDP. " + sc_lic)637 mv, mlab = mean_from(annual_mean(r))638 add("T12_household_net_worth_GDP_CA_avg", f"Canada: household net worth / GDP, {AVG_START}-latest average", mv, "ratio (years of GDP)", mlab.split(" mean")[0], "Canada", "annual average",639 "Statistics Canada, Tables 36-10-0580-01 and 36-10-0103-01", "Net worth / GDP", "https://www150.statcan.gc.ca/n1/tbl/csv/36100580-eng.zip", "data/raw/statcan/unzipped/36100580.csv", mlab + ". " + sc_lic)640 else:641 not_retrieved("T12_household_net_worth_GDP_CA_latest", "Canada household net worth / GDP", "https://www150.statcan.gc.ca/n1/tbl/csv/36100580-eng.zip", "data/raw/statcan/unzipped/36100580.csv", "file missing", geography="Canada")642else:643 not_retrieved("T12_labor_share_COE_GDP_CA_latest", "Canada COE/GDP", "https://www150.statcan.gc.ca/n1/tbl/csv/36100103-eng.zip", "data/raw/statcan/unzipped/36100103.csv", "file missing", geography="Canada")644645# ============================================================================646# 13. Electricity price647# ============================================================================648level_targets("T13_electricity_price_kWh", "Average price of electricity per kWh, US city average (BLS average price data)", "APU000072610", "US dollars per kWh", annual_from_hf=True)649650# ============================================================================651# BEA Input-Output: 2023 summary use table -> sector labor shares, intermediate shares652# ============================================================================653bea_p = RAW / "bea_io" / "unzipped" / "IOUse_After_Redefinitions_PRO_1997-2023_Summary.xlsx"654bea_url = "https://apps.bea.gov/industry/iTables%20Static%20Files/AllTablesIO.zip"655if bea_p.exists():656 u = pd.read_excel(bea_p, sheet_name="2023", header=None)657 code_row = u.index[u.iloc[:, 1].astype(str).str.strip() == "Commodities/Industries"][0]658 name_row = code_row + 1659 codes = u.iloc[code_row, 2:].astype(str).str.strip()660 names = u.iloc[name_row, 2:].astype(str).str.strip()661 # industry columns = those before the first final-use column (codes starting with 'F')662 ind_mask = ~codes.str.startswith("F") & (codes != "nan") & ~names.str.contains("Total", na=False)663 first_label = u.iloc[:, 1].astype(str).str.strip()664665 def row_vals(label):666 i = u.index[first_label == label][0]667 return pd.to_numeric(u.iloc[i, 2:], errors="coerce")668669 V = {k: row_vals(k) for k in ("Compensation of employees", "Taxes on production and imports, less subsidies", "Gross operating surplus", "Total Value Added", "Total Intermediate", "Total Industry Output")}670 out = pd.DataFrame({"industry_code": codes.values, "industry_name": names.values,671 "compensation_of_employees": V["Compensation of employees"].values,672 "taxes_less_subsidies": V["Taxes on production and imports, less subsidies"].values,673 "gross_operating_surplus": V["Gross operating surplus"].values,674 "total_value_added": V["Total Value Added"].values,675 "total_intermediate": V["Total Intermediate"].values,676 "total_industry_output": V["Total Industry Output"].values})[ind_mask.values]677 out = out.dropna(subset=["total_value_added"])678 out["labor_share_of_value_added"] = out["compensation_of_employees"] / out["total_value_added"]679 out["labor_share_of_factor_cost_VA"] = out["compensation_of_employees"] / (out["compensation_of_employees"] + out["gross_operating_surplus"])680 out["intermediate_share_of_gross_output"] = out["total_intermediate"] / out["total_industry_output"]681 out["value_added_share_of_total"] = out["total_value_added"] / out["total_value_added"].sum()682 out.to_csv(PROC / "bea_io_value_added_2023.csv", index=False)683 tot_va, tot_coe, tot_gos = out["total_value_added"].sum(), out["compensation_of_employees"].sum(), out["gross_operating_surplus"].sum()684 note = "Summary-level use table after redefinitions, producers' prices, millions of current dollars, 2023 (file dated 2024-09-06 inside AllTablesIO.zip). Sum over all 71 summary industries (private + government). Public domain."685 add("T01_labor_share_BEA_IO_2023_VA", "Compensation of employees / total value added, all industries, BEA use table 2023", tot_coe / tot_va, "ratio", "2023", "United States", "annual",686 "BEA Input-Output Accounts, IOUse_After_Redefinitions_PRO_1997-2023_Summary.xlsx", "V001/Total Value Added", bea_url, "data/raw/bea_io/unzipped/IOUse_After_Redefinitions_PRO_1997-2023_Summary.xlsx", note)687 add("T01_labor_share_BEA_IO_2023_factor_cost", "Compensation / (compensation + gross operating surplus), all industries, BEA use table 2023", tot_coe / (tot_coe + tot_gos), "ratio", "2023", "United States", "annual",688 "BEA Input-Output Accounts (summary use table 2023)", "V001/(V001+V003)", bea_url, "data/raw/bea_io/unzipped/IOUse_After_Redefinitions_PRO_1997-2023_Summary.xlsx", note + " Gross operating surplus includes CFC and proprietors' income.")689 add("T01_intermediate_share_BEA_IO_2023", "Total intermediate inputs / total industry output, all industries, BEA use table 2023", out["total_intermediate"].sum() / out["total_industry_output"].sum(), "ratio", "2023", "United States", "annual",690 "BEA Input-Output Accounts (summary use table 2023)", "Total Intermediate/Total Industry Output", bea_url, "data/raw/bea_io/unzipped/IOUse_After_Redefinitions_PRO_1997-2023_Summary.xlsx", note + " Per-industry values in data/processed/bea_io_value_added_2023.csv.")691 # a few sectors of direct interest to the model692 for code, tid in (("22", "T01_labor_share_BEA_IO_2023_utilities"), ("5415", "T01_labor_share_BEA_IO_2023_computer_systems_design"), ("513", "T01_labor_share_BEA_IO_2023_publishing_incl_software"),693 ("334", "T01_labor_share_BEA_IO_2023_computer_electronic_products"), ("531", "T01_labor_share_BEA_IO_2023_real_estate"), ("23", "T01_labor_share_BEA_IO_2023_construction"),694 ("5412OP", "T01_labor_share_BEA_IO_2023_professional_services"), ("621", "T01_labor_share_BEA_IO_2023_ambulatory_health")):695 r = out[out["industry_code"] == code]696 if len(r):697 r = r.iloc[0]698 add(tid, f"Labor share of value added, BEA summary industry {code} ({r['industry_name']}), 2023", float(r["labor_share_of_value_added"]), "ratio", "2023", "United States", "annual",699 "BEA Input-Output Accounts (summary use table 2023)", f"V001/Total Value Added, column {code}", bea_url, "data/raw/bea_io/unzipped/IOUse_After_Redefinitions_PRO_1997-2023_Summary.xlsx",700 f"Compensation {r['compensation_of_employees']:.0f} / VA {r['total_value_added']:.0f} ($mn).")701else:702 not_retrieved("T01_labor_share_BEA_IO_2023_VA", "BEA IO labor share 2023", bea_url, "data/raw/bea_io/", "AllTablesIO.zip missing")703704# ============================================================================705# BLS API cross-checks (independent retrieval path for PRS85006173, APU000072610, CUSR0000SEHA)706# ============================================================================707for sid, desc in (("PRS85006173", "BLS nonfarm business labor share index (BLS API v1)"), ("APU000072610", "Electricity price per kWh (BLS API v1)"), ("CUSR0000SEHA", "CPI rent of primary residence (BLS API v1)")):708 files = sorted((RAW / "bls").glob(f"{sid}_*.json"))709 if not files:710 continue711 recs = []712 for f in files:713 js = json.loads(f.read_text())714 for srs in js.get("Results", {}).get("series", []):715 for o in srs.get("data", []):716 p = o["period"]717 if p.startswith("M") and p != "M13":718 dt = f"{o['year']}-{int(p[1:]):02d}-01"719 elif p.startswith("Q") and p != "Q05":720 dt = f"{o['year']}-{(int(p[1:]) - 1) * 3 + 1:02d}-01"721 else:722 continue723 try:724 recs.append((dt, float(o["value"])))725 except ValueError: # BLS uses '-' for missing726 continue727 s = pd.Series(dict(recs))728 s.index = pd.to_datetime(s.index)729 s = s.sort_index()730 save_series(f"bls_api_{sid}", s)731 d, v = latest(s)732 add(f"T00_crosscheck_BLS_API_{sid}_latest", desc + ", latest", v, fred_meta.get(sid, {}).get("units", ""), period_label(d, fred_meta.get(sid, {}).get("frequency", "")), "United States", fred_meta.get(sid, {}).get("frequency", ""),733 "U.S. Bureau of Labor Statistics public API v1 (key-less)", sid, "https://api.bls.gov/publicAPI/v1/timeseries/data/ (POST, 10-year windows)", f"data/raw/bls/{sid}_<years>.json",734 "Independent retrieval path; should equal the FRED value for the same date.")735736# ----------------------------------------------------------------------------737PROC.mkdir(parents=True, exist_ok=True)738with open(PROC / "calibration_targets.csv", "w", newline="") as f:739 w = csv.DictWriter(f, fieldnames=COLS)740 w.writeheader()741 for r in ROWS:742 w.writerow(r)743n_ok = sum(r["status"] == "RETRIEVED" for r in ROWS)744n_nr = sum(r["status"] == "NOT_RETRIEVED" for r in ROWS)745print(f"calibration_targets.csv: {len(ROWS)} rows, {n_ok} RETRIEVED, {n_nr} NOT_RETRIEVED")746for r in ROWS:747 if r["status"] == "NOT_RETRIEVED":748 print(" NOT_RETRIEVED:", r["target_id"], "-", r["notes"])749print(f"series written: {len(list(SER.glob('*.csv')))}")750