1#!/usr/bin/env python32"""02_build_calibration_targets.py -- compute calibration moments from data/raw/.34Writes data/processed/calibration_targets.csv with columns5 target, value, units, period, source_series, method6Every row is computed from a file listed in data/raw/MANIFEST.json; nothing is7typed in by hand. Where a moment could not be computed from the downloaded8data (e.g. 2008/2020 credit-spread peaks, FRED window restriction) the row is9kept with value NaN and the reason in `method` so the gap is visible.1011Run with: python3 -I scripts/02_build_calibration_targets.py12Requires pandas, numpy, openpyxl (Fernald xlsx), xlrd (Shiller xls).13"""14from __future__ import annotations1516import math17from pathlib import Path1819import numpy as np20import pandas as pd2122ROOT = Path(__file__).resolve().parents[1]23RAW = ROOT / "data" / "raw"24OUT = ROOT / "data" / "processed" / "calibration_targets.csv"25OUT.parent.mkdir(parents=True, exist_ok=True)2627rows: list[dict] = []282930def add(target, value, units, period, source, method):31 rows.append(dict(target=target, value=(None if value is None or (isinstance(value, float) and math.isnan(value)) else float(value)),32 units=units, period=period, source_series=source, method=method))333435# ----------------------------------------------------------------------------36# helpers37# ----------------------------------------------------------------------------38def fred(sid: str) -> pd.Series:39 df = pd.read_csv(RAW / "fred" / f"{sid}.csv", na_values=["."])40 s = pd.Series(df[sid].values, index=pd.to_datetime(df["observation_date"]), name=sid).astype(float)41 return s.dropna()424344def annual_mean(s: pd.Series) -> pd.Series:45 a = s.groupby(s.index.year).mean()46 a.index.name = "year"47 return a484950def hp_filter(y: np.ndarray, lam: float) -> np.ndarray:51 """Return the cyclical component of y (Hodrick-Prescott, dense solve)."""52 n = len(y)53 D = np.zeros((n - 2, n))54 for i in range(n - 2):55 D[i, i], D[i, i + 1], D[i, i + 2] = 1.0, -2.0, 1.056 trend = np.linalg.solve(np.eye(n) + lam * D.T @ D, y)57 return y - trend585960def ar1(x: np.ndarray) -> tuple[float, float]:61 """OLS AR(1) with intercept on a (detrended) series: returns (rho, sd of residual)."""62 y, xl = x[1:], x[:-1]63 X = np.column_stack([np.ones_like(xl), xl])64 beta, *_ = np.linalg.lstsq(X, y, rcond=None)65 resid = y - X @ beta66 return float(beta[1]), float(np.std(resid, ddof=2))676869def win(s: pd.Series, y0: int, y1: int) -> pd.Series:70 return s[(s.index >= y0) & (s.index <= y1)]717273# ----------------------------------------------------------------------------74# 1. Labor share and capital share75# ----------------------------------------------------------------------------76prs = annual_mean(fred("PRS85006173")) # index, 2017 = 10077pwt = fred("LABSHPUSA156NRUG"); pwt.index = pwt.index.year # level (ratio), PWT78bea_comp = fred("A4002E1A156NBEA"); bea_comp.index = bea_comp.index.year # percent of GDI79bea_wage = fred("W270RE1A156NBEA"); bea_wage.index = bea_wage.index.year80last_prs_year = int(prs.index.max())8182anchor = float(pwt.loc[2017]) # PWT labour-compensation share in 2017 (index base year)83prs_level = prs / 100.0 * anchor8485for (y0, y1) in [(1947, last_prs_year), (1984, last_prs_year), (2015, last_prs_year)]:86 add(f"labor_share_index_mean_{y0}_{y1}", win(prs, y0, y1).mean(), "index (2017=100)", f"{y0}-{y1}",87 "FRED:PRS85006173", "annual average of quarterly BLS nonfarm-business labor-share index")88 add(f"labor_share_level_mean_{y0}_{y1}", win(prs_level, y0, y1).mean(), "ratio", f"{y0}-{y1}",89 "FRED:PRS85006173 x FRED:LABSHPUSA156NRUG(2017)",90 f"BLS index rescaled: level_t = index_t/100 x PWT labsh(2017)={anchor:.4f}. ASSUMPTION: the PWT 2017 economy-wide labour-compensation share is used as the level of the BLS nonfarm-business share in its base year; the two concepts differ (PWT includes imputed self-employment labour income), so the rescaled level is indicative only.")91add(f"labor_share_level_latest_{last_prs_year}", prs_level.loc[last_prs_year], "ratio", str(last_prs_year), "FRED:PRS85006173 rescaled", "same rescaling, latest full/partial year average")92add("labor_share_pwt_mean_1984_2023", win(pwt, 1984, 2023).mean(), "ratio", "1984-2023", "FRED:LABSHPUSA156NRUG", "Penn World Table labsh, annual mean (level, no rescaling)")93add("compensation_share_gdi_mean_1984_2025", win(bea_comp, 1984, 2025).mean(), "percent of GDI", "1984-2025", "FRED:A4002E1A156NBEA", "BEA compensation of employees / gross domestic income, annual mean")94add("compensation_share_gdi_latest_2025", bea_comp.loc[2025], "percent of GDI", "2025", "FRED:A4002E1A156NBEA", "latest annual value")95add("wage_salary_share_gdi_latest_2025", bea_wage.loc[2025], "percent of GDI", "2025", "FRED:W270RE1A156NBEA", "wages and salaries only (excludes supplements)")96add(f"capital_share_1_minus_labor_level_1984_{last_prs_year}", 1 - win(prs_level, 1984, last_prs_year).mean(), "ratio", f"1984-{last_prs_year}", "derived", "1 - rescaled BLS labor-share level (gross capital share incl. depreciation)")9798# ----------------------------------------------------------------------------99# 2. Fernald TFP: alpha (capital share) and TFP process100# ----------------------------------------------------------------------------101fern_a = pd.read_excel(RAW / "fernald" / "quarterly_tfp.xlsx", sheet_name="annual")102fern_a = fern_a.dropna(subset=["date"]).copy()103fern_a["date"] = fern_a["date"].astype(int)104fern_a = fern_a.set_index("date")105last_f = int(fern_a.dropna(subset=["dtfp_util"]).index.max())106add(f"capital_share_fernald_alpha_mean_1984_{last_f}", win(fern_a["alpha"], 1984, last_f).mean(), "ratio", f"1984-{last_f}", "Fernald quarterly_tfp.xlsx (annual, col alpha)", "annual mean of Fernald's capital share (business sector)")107add(f"capital_share_fernald_alpha_mean_1948_{last_f}", win(fern_a["alpha"], 1948, last_f).mean(), "ratio", f"1948-{last_f}", "Fernald quarterly_tfp.xlsx (annual, col alpha)", "annual mean")108109for col, label in [("dtfp_util", "util_adjusted"), ("dtfp", "raw")]:110 g = fern_a[col].dropna() / 100.0 # annual log change111 for (y0, y1) in [(1948, last_f), (1984, last_f)]:112 gg = win(g, y0, y1)113 add(f"tfp_{label}_growth_mean_{y0}_{y1}", gg.mean(), "log change per year", f"{y0}-{y1}", f"Fernald annual {col}", "mean of 100*dlog TFP / 100")114 add(f"tfp_{label}_growth_sd_{y0}_{y1}", gg.std(ddof=1), "log change per year", f"{y0}-{y1}", f"Fernald annual {col}", "first-difference method: sd of annual log TFP growth (= innovation sd if log TFP is a random walk with drift)")115 lvl = np.cumsum(gg.values)116 for lam in (6.25, 100.0):117 rho, sd = ar1(hp_filter(lvl, lam))118 add(f"tfp_{label}_ar1_rho_hp{lam:g}_{y0}_{y1}", rho, "per year", f"{y0}-{y1}", f"Fernald annual {col}", f"HP filter (lambda={lam:g}) of cumulated log TFP, OLS AR(1) with intercept on the cyclical component")119 add(f"tfp_{label}_ar1_sd_hp{lam:g}_{y0}_{y1}", sd, "log units per year", f"{y0}-{y1}", f"Fernald annual {col}", "residual sd of the same AR(1)")120 t = np.arange(len(lvl))121 lin = lvl - np.polyval(np.polyfit(t, lvl, 1), t)122 rho, sd = ar1(lin)123 add(f"tfp_{label}_ar1_rho_lintrend_{y0}_{y1}", rho, "per year", f"{y0}-{y1}", f"Fernald annual {col}", "linear detrending of cumulated log TFP, OLS AR(1)")124 add(f"tfp_{label}_ar1_sd_lintrend_{y0}_{y1}", sd, "log units per year", f"{y0}-{y1}", f"Fernald annual {col}", "residual sd of the same AR(1)")125126# quarterly cross-check127fern_q = pd.read_excel(RAW / "fernald" / "quarterly_tfp.xlsx", sheet_name="quarterly", header=1)128fern_q = fern_q.dropna(subset=["dtfp_util"]).copy()129fern_q = fern_q[fern_q["date"].astype(str).str.fullmatch(r"\d{4}:Q[1-4]")].copy() # drop summary rows130fern_q["year"] = fern_q["date"].astype(str).str[:4].astype(int)131fq = fern_q[(fern_q["year"] >= 1984)]132lvl_q = np.cumsum(fq["dtfp_util"].values / 400.0)133rho_q, sd_q = ar1(hp_filter(lvl_q, 1600.0))134add(f"tfp_util_adjusted_ar1_rho_quarterly_hp1600_1984_{int(fq['year'].max())}", rho_q, "per quarter", f"1984Q1-{fq['date'].iloc[-1]}", "Fernald quarterly dtfp_util", "HP(1600) cyclical component of cumulated quarterly log TFP, AR(1)")135add(f"tfp_util_adjusted_ar1_rho_quarterly_annualised_1984_{int(fq['year'].max())}", rho_q ** 4, "per year", f"1984Q1-{fq['date'].iloc[-1]}", "Fernald quarterly dtfp_util", "rho_q^4")136add(f"tfp_util_adjusted_ar1_sd_quarterly_hp1600_1984_{int(fq['year'].max())}", sd_q, "log units per quarter", f"1984Q1-{fq['date'].iloc[-1]}", "Fernald quarterly dtfp_util", "residual sd")137138# ----------------------------------------------------------------------------139# 3. Growth, investment and consumption shares, depreciation140# ----------------------------------------------------------------------------141gdp = annual_mean(fred("GDPC1")); pce = annual_mean(fred("PCECC96")); inv = annual_mean(fred("GPDIC1"))142last_full = int(gdp.index.max()) - 1 # last complete year (current year has only part of the quarters)143g_ann = np.log(gdp).diff().dropna()144for (y0, y1) in [(1948, last_full), (1984, last_full)]:145 add(f"real_gdp_growth_mean_{y0}_{y1}", win(g_ann, y0, y1).mean(), "log change per year", f"{y0}-{y1}", "FRED:GDPC1", "log change of annual-average real GDP, mean")146 add(f"real_gdp_growth_sd_{y0}_{y1}", win(g_ann, y0, y1).std(ddof=1), "log change per year", f"{y0}-{y1}", "FRED:GDPC1", "sd of annual log growth")147gq = fred("A191RL1Q225SBEA")148add(f"real_gdp_growth_saar_mean_1984_{last_full}", gq[(gq.index.year >= 1984) & (gq.index.year <= last_full)].mean(), "percent, annualised", f"1984Q1-{last_full}Q4", "FRED:A191RL1Q225SBEA", "mean of BEA quarterly annualised growth rates")149for (y0, y1) in [(1984, last_full), (2015, last_full)]:150 add(f"investment_gdp_ratio_mean_{y0}_{y1}", win(inv / gdp, y0, y1).mean(), "ratio", f"{y0}-{y1}", "FRED:GPDIC1 / FRED:GDPC1", "gross private domestic investment over GDP, both chained-2017$ (ratio of chained aggregates is approximate; nominal shares differ slightly)")151 add(f"consumption_gdp_ratio_mean_{y0}_{y1}", win(pce / gdp, y0, y1).mean(), "ratio", f"{y0}-{y1}", "FRED:PCECC96 / FRED:GDPC1", "personal consumption over GDP, chained-2017$ (same caveat)")152153K = fred("K1TTOTL1ES000"); K.index = K.index.year154M = fred("M1TTOTL1ES000"); M.index = M.index.year155dep = (M / K.shift(1)).dropna()156add("depreciation_rate_private_fixed_assets_mean_1984_2024", win(dep, 1984, 2024).mean(), "per year", "1984-2024", "FRED:M1TTOTL1ES000 / FRED:K1TTOTL1ES000(t-1)", "BEA current-cost depreciation of private fixed assets divided by previous year-end current-cost net stock")157add("depreciation_rate_private_fixed_assets_2024", dep.loc[2024], "per year", "2024", "same", "latest year")158add("depreciation_rate_private_fixed_assets_mean_2015_2024", win(dep, 2015, 2024).mean(), "per year", "2015-2024", "same", "recent decade")159add("capital_output_ratio_private_fixed_assets_2024", None, "ratio", "2024", "FRED:K1TTOTL1ES000", "NOT COMPUTED here: K is current-cost (nominal) while GDPC1 is chained 2017$; nominal GDP not downloaded. Modeling team may add GDPA if needed.")160161# ----------------------------------------------------------------------------162# 4. Interest rates and inflation163# ----------------------------------------------------------------------------164cpi = annual_mean(fred("CPIAUCSL")); infl = (np.log(cpi).diff() * 100).dropna()165tb3 = annual_mean(fred("TB3MS")); dff = annual_mean(fred("DFF")); gs10 = annual_mean(fred("DGS10")); tips = annual_mean(fred("DFII10"))166last_cpi = int(infl.index.max()) - 1 # last full year for annual-average inflation167add(f"cpi_inflation_mean_1984_{last_cpi}", win(infl, 1984, last_cpi).mean(), "percent per year", f"1984-{last_cpi}", "FRED:CPIAUCSL", "log change of annual-average CPI")168rf_real = (tb3 - infl).dropna()169add(f"real_tbill_rate_expost_mean_1984_{last_cpi}", win(rf_real, 1984, last_cpi).mean(), "percent per year", f"1984-{last_cpi}", "FRED:TB3MS - FRED:CPIAUCSL", "annual-average 3-month T-bill rate minus same-year CPI inflation (ex post)")170add(f"real_tbill_rate_expost_sd_1984_{last_cpi}", win(rf_real, 1984, last_cpi).std(ddof=1), "percent per year", f"1984-{last_cpi}", "same", "sd of annual ex-post real T-bill rate")171add(f"real_tbill_rate_expost_mean_1948_{last_cpi}", win(rf_real, 1948, last_cpi).mean(), "percent per year", f"1948-{last_cpi}", "same", "post-war mean (TB3MS starts 1934, CPI 1947; window 1948-)")172add(f"real_fedfunds_rate_expost_mean_1984_{last_cpi}", win(dff - infl, 1984, last_cpi).mean(), "percent per year", f"1984-{last_cpi}", "FRED:DFF - FRED:CPIAUCSL", "annual-average effective federal funds rate minus CPI inflation")173add(f"real_10y_treasury_expost_mean_1984_{last_cpi}", win(gs10 - infl, 1984, last_cpi).mean(), "percent per year", f"1984-{last_cpi}", "FRED:DGS10 - FRED:CPIAUCSL", "annual-average 10y nominal yield minus same-year CPI inflation (crude ex post)")174last_tips = int(tips.index.max()) - 1175add(f"tips_10y_real_yield_mean_2003_{last_tips}", win(tips, 2003, last_tips).mean(), "percent per year", f"2003-{last_tips}", "FRED:DFII10", "mean of annual averages of daily 10-year TIPS yield (ex ante real)")176add(f"tips_10y_real_yield_mean_2015_{last_tips}", win(tips, 2015, last_tips).mean(), "percent per year", f"2015-{last_tips}", "FRED:DFII10", "recent decade")177tips_d = fred("DFII10")178add("tips_10y_real_yield_latest", tips_d.iloc[-1], "percent", str(tips_d.index[-1].date()), "FRED:DFII10", "latest daily observation")179180# ----------------------------------------------------------------------------181# 5. Credit spreads, loan standards, bank credit, bank leverage182# ----------------------------------------------------------------------------183for sid, lab in [("BAMLC0A0CM", "ig_oas"), ("BAMLH0A0HYM2", "hy_oas")]:184 s = fred(sid)185 per = f"{s.index[0].date()}-{s.index[-1].date()}"186 add(f"{lab}_mean_window", s.mean(), "percent", per, f"FRED:{sid}", "mean over the only downloadable window (FRED restricts ICE BofA series to 3 years since April 2026)")187 add(f"{lab}_max_window", s.max(), "percent", f"{s.idxmax().date()}", f"FRED:{sid}", "max within the window")188 add(f"{lab}_latest", s.iloc[-1], "percent", str(s.index[-1].date()), f"FRED:{sid}", "latest")189 add(f"{lab}_peak_2008", None, "percent", "2008", f"FRED:{sid}", "NOT OBTAINED: pre-2023 history not downloadable from FRED (ICE licence); Fed FSR Nov-2020 chart reads ~2000 bp HY / ~680 bp BBB in late 2008 -- see ledger")190 add(f"{lab}_peak_2020", None, "percent", "2020", f"FRED:{sid}", "NOT OBTAINED: same restriction; Fed FSR Nov-2020 chart reads HY ~800 bp, BBB ~400 bp in March 2020 -- see ledger")191192sl = fred("DRTSCILM")193add("sloos_ci_tightening_mean_1990_2025", sl[(sl.index.year >= 1990) & (sl.index.year <= 2025)].mean(), "net percent of banks tightening", "1990Q2-2025Q4", "FRED:DRTSCILM", "mean")194add("sloos_ci_tightening_peak_2008_2009", sl[(sl.index.year >= 2008) & (sl.index.year <= 2009)].max(), "net percent", str(sl[(sl.index.year >= 2008) & (sl.index.year <= 2009)].idxmax().date()), "FRED:DRTSCILM", "max in 2008-2009")195add("sloos_ci_tightening_peak_2020", sl[sl.index.year == 2020].max(), "net percent", str(sl[sl.index.year == 2020].idxmax().date()), "FRED:DRTSCILM", "max in 2020")196bc = annual_mean(fred("TOTBKCR")); bcg = np.log(bc).diff().dropna()197add(f"bank_credit_nominal_growth_mean_1984_{last_full}", win(bcg, 1984, last_full).mean(), "log change per year", f"1984-{last_full}", "FRED:TOTBKCR", "annual-average bank credit, log growth (nominal)")198add(f"bank_credit_nominal_growth_sd_1984_{last_full}", win(bcg, 1984, last_full).std(ddof=1), "log change per year", f"1984-{last_full}", "FRED:TOTBKCR", "sd")199eq = fred("RALACBW027SBOG"); ta = fred("TLAACBW027SBOG")200lev = (eq / ta).dropna(); lev_a = annual_mean(lev)201add("bank_equity_to_assets_latest", lev.iloc[-1], "ratio", str(lev.index[-1].date()), "FRED:RALACBW027SBOG / FRED:TLAACBW027SBOG", "H.8 residual (assets less liabilities) over total assets, all commercial banks (book, not regulatory)")202add(f"bank_equity_to_assets_mean_2010_{last_full}", win(lev_a, 2010, last_full).mean(), "ratio", f"2010-{last_full}", "same", "mean of annual averages")203add(f"bank_equity_to_assets_mean_1984_{last_full}", win(lev_a, 1984, last_full).mean(), "ratio", f"1984-{last_full}", "same", "mean of annual averages")204add("bank_assets_to_equity_leverage_latest", 1 / lev.iloc[-1], "x", str(lev.index[-1].date()), "same", "inverse of equity/assets")205206# ----------------------------------------------------------------------------207# 6. Wealth distribution (Fed DFA), latest quarter208# ----------------------------------------------------------------------------209dfa = pd.read_csv(RAW / "dfa" / "dfa-networth-shares.csv")210latest = dfa["Date"].iloc[-1]211d = dfa[dfa["Date"] == latest].set_index("Category")["Net worth"]212top1 = d["TopPt1"] + d["RemainingTop1"]; top10 = top1 + d["Next9"]; bot50 = d["Bottom50"]213add("wealth_share_top0p1_dfa_latest", d["TopPt1"], "percent of household net worth", latest, "DFA dfa-networth-shares.csv", "Fed Distributional Financial Accounts, latest quarter")214add("wealth_share_top1_dfa_latest", top1, "percent of household net worth", latest, "DFA dfa-networth-shares.csv", "TopPt1 + RemainingTop1")215add("wealth_share_top10_dfa_latest", top10, "percent of household net worth", latest, "DFA dfa-networth-shares.csv", "top1 + Next9")216add("wealth_share_next40_dfa_latest", d["Next40"], "percent of household net worth", latest, "DFA dfa-networth-shares.csv", "50th-90th percentile")217add("wealth_share_bottom50_dfa_latest", bot50, "percent of household net worth", latest, "DFA dfa-networth-shares.csv", "bottom half")218first = dfa["Date"].iloc[0]; d0 = dfa[dfa["Date"] == first].set_index("Category")["Net worth"]219add("wealth_share_top1_dfa_1989Q3", d0["TopPt1"] + d0["RemainingTop1"], "percent", first, "DFA dfa-networth-shares.csv", "first quarter available, for comparison")220add("wealth_share_bottom50_dfa_1989Q3", d0["Bottom50"], "percent", first, "DFA dfa-networth-shares.csv", "first quarter available")221# grouped-data Gini (lower bound: assumes equality within each of the 5 groups)222p = np.array([0.50, 0.40, 0.09, 0.009, 0.001])223w = np.array([bot50, d["Next40"], d["Next9"], d["RemainingTop1"], d["TopPt1"]]) / 100.0224L = np.concatenate([[0.0], np.cumsum(w)])225gini_lb = 1.0 - float(np.sum((L[:-1] + L[1:]) * p))226add("wealth_gini_dfa_grouped_lower_bound_latest", gini_lb, "Gini", latest, "DFA dfa-networth-shares.csv", "Gini from 5 grouped shares (bottom50/next40/next9/next0.9/top0.1); within-group equality assumed, so this is a LOWER BOUND on the true Gini (SCF-based estimates ~0.85)")227228# ----------------------------------------------------------------------------229# 7. Equity returns, premium and volatility (Shiller monthly data)230# ----------------------------------------------------------------------------231sh = pd.read_excel(RAW / "shiller" / "ie_data.xls", sheet_name="Data", header=None)232sh = sh.iloc[8:].copy()233sh = sh[pd.to_numeric(sh[0], errors="coerce").notna()]234sh = sh.rename(columns={0: "date", 1: "P", 2: "D", 4: "CPI", 6: "GS10"})[["date", "P", "D", "CPI", "GS10"]].astype(float)235sh = sh.reset_index(drop=True)236assert abs(sh.loc[0, "date"] - 1871.01) < 1e-6237sh["month"] = pd.period_range("1871-01", periods=len(sh), freq="M")238sh["year"] = sh["month"].dt.year239sh["ret_nom"] = (sh["P"] + sh["D"].ffill() / 12.0) / sh["P"].shift(1) - 1.0240ann = sh.groupby("year").agg(ret=("ret_nom", lambda r: np.prod(1 + r.values) - 1), n=("ret_nom", "size"),241 cpi_end=("CPI", "last"), gs10=("GS10", "mean"))242ann = ann[ann["n"] == 12]243ann["infl"] = ann["cpi_end"] / ann["cpi_end"].shift(1) - 1.0244ann["ret_real"] = (1 + ann["ret"]) / (1 + ann["infl"]) - 1.0245ann = ann.dropna(subset=["ret_real"])246last_sh = int(ann.index.max())247for (y0, y1) in [(1872, last_sh), (1926, last_sh), (1984, last_sh)]:248 a = win(ann, y0, y1)249 add(f"equity_real_return_mean_{y0}_{y1}", a["ret_real"].mean() * 100, "percent per year", f"{y0}-{y1}", "Shiller ie_data.xls (P, D, CPI)", "S&P composite total return (monthly P + D/12, compounded over calendar year), deflated by Dec/Dec CPI, arithmetic mean")250 add(f"equity_real_return_sd_{y0}_{y1}", a["ret_real"].std(ddof=1) * 100, "percent per year", f"{y0}-{y1}", "Shiller ie_data.xls", "sd of annual real total return")251 add(f"equity_premium_over_gs10_{y0}_{y1}", (a["ret"] - a["gs10"] / 100).mean() * 100, "percent per year", f"{y0}-{y1}", "Shiller ie_data.xls (ret, GS10)", "nominal equity total return minus annual-average 10-year Treasury yield (NB: yield, not bond holding-period return)")252tb3_full = annual_mean(fred("TB3MS"))253joint = ann.join(tb3_full.rename("tb3"), how="inner")254for (y0, y1) in [(1934, last_sh), (1948, last_sh), (1984, last_sh)]:255 j = win(joint, y0, y1)256 add(f"equity_premium_over_tbill_{y0}_{y1}", (j["ret"] - j["tb3"] / 100).mean() * 100, "percent per year", f"{y0}-{y1}", "Shiller ie_data.xls + FRED:TB3MS", "nominal equity total return minus annual-average 3-month T-bill rate, arithmetic mean")257 add(f"equity_premium_over_tbill_sd_{y0}_{y1}", (j["ret"] - j["tb3"] / 100).std(ddof=1) * 100, "percent per year", f"{y0}-{y1}", "same", "sd of annual excess return")258 add(f"sharpe_ratio_equity_{y0}_{y1}", (j["ret"] - j["tb3"] / 100).mean() / (j["ret"] - j["tb3"] / 100).std(ddof=1), "ratio", f"{y0}-{y1}", "same", "mean / sd of annual excess return")259260sp = fred("SP500"); sp_a = sp.groupby(sp.index.year).last(); sp_r = np.log(sp_a).diff().dropna()261sp_r = sp_r[sp_r.index <= last_full]262add(f"sp500_price_log_return_mean_{int(sp_r.index.min())}_{int(sp_r.index.max())}", sp_r.mean() * 100, "percent per year", f"{int(sp_r.index.min())}-{int(sp_r.index.max())}", "FRED:SP500", "year-end to year-end log price return (no dividends), 10-year FRED window")263add(f"sp500_price_log_return_sd_{int(sp_r.index.min())}_{int(sp_r.index.max())}", sp_r.std(ddof=1) * 100, "percent per year", f"{int(sp_r.index.min())}-{int(sp_r.index.max())}", "FRED:SP500", "sd")264# FRED:BOGZ1FL073164003Q (Z.1 "NYSE Composite Index, Level") is downloaded but NOT used: its values265# (e.g. 2,383,423 in 2026Q2) are not an index level in the usual sense and the FRED unit label is "Percent";266# treating it as a price index gave an implausible 30% annual sd. Left out until the concept is clarified.267268# ----------------------------------------------------------------------------269# write270# ----------------------------------------------------------------------------271df = pd.DataFrame(rows, columns=["target", "value", "units", "period", "source_series", "method"])272df.to_csv(OUT, index=False, float_format="%.6g")273print(f"wrote {OUT} ({len(df)} rows; {df['value'].isna().sum()} NOT OBTAINED/NOT COMPUTED)")274with pd.option_context("display.max_rows", 500, "display.width", 250, "display.max_colwidth", 60):275 print(df[["target", "value", "units", "period"]].to_string(index=False))276