#!/usr/bin/env python3 """02_build_calibration_targets.py -- compute calibration moments from data/raw/. Writes data/processed/calibration_targets.csv with columns target, value, units, period, source_series, method Every row is computed from a file listed in data/raw/MANIFEST.json; nothing is typed in by hand. Where a moment could not be computed from the downloaded data (e.g. 2008/2020 credit-spread peaks, FRED window restriction) the row is kept with value NaN and the reason in `method` so the gap is visible. Run with: python3 -I scripts/02_build_calibration_targets.py Requires pandas, numpy, openpyxl (Fernald xlsx), xlrd (Shiller xls). """ from __future__ import annotations import math from pathlib import Path import numpy as np import pandas as pd ROOT = Path(__file__).resolve().parents[1] RAW = ROOT / "data" / "raw" OUT = ROOT / "data" / "processed" / "calibration_targets.csv" OUT.parent.mkdir(parents=True, exist_ok=True) rows: list[dict] = [] def add(target, value, units, period, source, method): rows.append(dict(target=target, value=(None if value is None or (isinstance(value, float) and math.isnan(value)) else float(value)), units=units, period=period, source_series=source, method=method)) # ---------------------------------------------------------------------------- # helpers # ---------------------------------------------------------------------------- def fred(sid: str) -> pd.Series: df = pd.read_csv(RAW / "fred" / f"{sid}.csv", na_values=["."]) s = pd.Series(df[sid].values, index=pd.to_datetime(df["observation_date"]), name=sid).astype(float) return s.dropna() def annual_mean(s: pd.Series) -> pd.Series: a = s.groupby(s.index.year).mean() a.index.name = "year" return a def hp_filter(y: np.ndarray, lam: float) -> np.ndarray: """Return the cyclical component of y (Hodrick-Prescott, dense solve).""" n = len(y) D = np.zeros((n - 2, n)) for i in range(n - 2): D[i, i], D[i, i + 1], D[i, i + 2] = 1.0, -2.0, 1.0 trend = np.linalg.solve(np.eye(n) + lam * D.T @ D, y) return y - trend def ar1(x: np.ndarray) -> tuple[float, float]: """OLS AR(1) with intercept on a (detrended) series: returns (rho, sd of residual).""" y, xl = x[1:], x[:-1] X = np.column_stack([np.ones_like(xl), xl]) beta, *_ = np.linalg.lstsq(X, y, rcond=None) resid = y - X @ beta return float(beta[1]), float(np.std(resid, ddof=2)) def win(s: pd.Series, y0: int, y1: int) -> pd.Series: return s[(s.index >= y0) & (s.index <= y1)] # ---------------------------------------------------------------------------- # 1. Labor share and capital share # ---------------------------------------------------------------------------- prs = annual_mean(fred("PRS85006173")) # index, 2017 = 100 pwt = fred("LABSHPUSA156NRUG"); pwt.index = pwt.index.year # level (ratio), PWT bea_comp = fred("A4002E1A156NBEA"); bea_comp.index = bea_comp.index.year # percent of GDI bea_wage = fred("W270RE1A156NBEA"); bea_wage.index = bea_wage.index.year last_prs_year = int(prs.index.max()) anchor = float(pwt.loc[2017]) # PWT labour-compensation share in 2017 (index base year) prs_level = prs / 100.0 * anchor for (y0, y1) in [(1947, last_prs_year), (1984, last_prs_year), (2015, last_prs_year)]: add(f"labor_share_index_mean_{y0}_{y1}", win(prs, y0, y1).mean(), "index (2017=100)", f"{y0}-{y1}", "FRED:PRS85006173", "annual average of quarterly BLS nonfarm-business labor-share index") add(f"labor_share_level_mean_{y0}_{y1}", win(prs_level, y0, y1).mean(), "ratio", f"{y0}-{y1}", "FRED:PRS85006173 x FRED:LABSHPUSA156NRUG(2017)", 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.") add(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") add("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)") add("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") add("compensation_share_gdi_latest_2025", bea_comp.loc[2025], "percent of GDI", "2025", "FRED:A4002E1A156NBEA", "latest annual value") add("wage_salary_share_gdi_latest_2025", bea_wage.loc[2025], "percent of GDI", "2025", "FRED:W270RE1A156NBEA", "wages and salaries only (excludes supplements)") add(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)") # ---------------------------------------------------------------------------- # 2. Fernald TFP: alpha (capital share) and TFP process # ---------------------------------------------------------------------------- fern_a = pd.read_excel(RAW / "fernald" / "quarterly_tfp.xlsx", sheet_name="annual") fern_a = fern_a.dropna(subset=["date"]).copy() fern_a["date"] = fern_a["date"].astype(int) fern_a = fern_a.set_index("date") last_f = int(fern_a.dropna(subset=["dtfp_util"]).index.max()) add(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)") add(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") for col, label in [("dtfp_util", "util_adjusted"), ("dtfp", "raw")]: g = fern_a[col].dropna() / 100.0 # annual log change for (y0, y1) in [(1948, last_f), (1984, last_f)]: gg = win(g, y0, y1) 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") 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)") lvl = np.cumsum(gg.values) for lam in (6.25, 100.0): rho, sd = ar1(hp_filter(lvl, lam)) 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") 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)") t = np.arange(len(lvl)) lin = lvl - np.polyval(np.polyfit(t, lvl, 1), t) rho, sd = ar1(lin) 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)") 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)") # quarterly cross-check fern_q = pd.read_excel(RAW / "fernald" / "quarterly_tfp.xlsx", sheet_name="quarterly", header=1) fern_q = fern_q.dropna(subset=["dtfp_util"]).copy() fern_q = fern_q[fern_q["date"].astype(str).str.fullmatch(r"\d{4}:Q[1-4]")].copy() # drop summary rows fern_q["year"] = fern_q["date"].astype(str).str[:4].astype(int) fq = fern_q[(fern_q["year"] >= 1984)] lvl_q = np.cumsum(fq["dtfp_util"].values / 400.0) rho_q, sd_q = ar1(hp_filter(lvl_q, 1600.0)) add(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)") add(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") add(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") # ---------------------------------------------------------------------------- # 3. Growth, investment and consumption shares, depreciation # ---------------------------------------------------------------------------- gdp = annual_mean(fred("GDPC1")); pce = annual_mean(fred("PCECC96")); inv = annual_mean(fred("GPDIC1")) last_full = int(gdp.index.max()) - 1 # last complete year (current year has only part of the quarters) g_ann = np.log(gdp).diff().dropna() for (y0, y1) in [(1948, last_full), (1984, last_full)]: 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") 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") gq = fred("A191RL1Q225SBEA") add(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") for (y0, y1) in [(1984, last_full), (2015, last_full)]: 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)") 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)") K = fred("K1TTOTL1ES000"); K.index = K.index.year M = fred("M1TTOTL1ES000"); M.index = M.index.year dep = (M / K.shift(1)).dropna() add("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") add("depreciation_rate_private_fixed_assets_2024", dep.loc[2024], "per year", "2024", "same", "latest year") add("depreciation_rate_private_fixed_assets_mean_2015_2024", win(dep, 2015, 2024).mean(), "per year", "2015-2024", "same", "recent decade") add("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.") # ---------------------------------------------------------------------------- # 4. Interest rates and inflation # ---------------------------------------------------------------------------- cpi = annual_mean(fred("CPIAUCSL")); infl = (np.log(cpi).diff() * 100).dropna() tb3 = annual_mean(fred("TB3MS")); dff = annual_mean(fred("DFF")); gs10 = annual_mean(fred("DGS10")); tips = annual_mean(fred("DFII10")) last_cpi = int(infl.index.max()) - 1 # last full year for annual-average inflation add(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") rf_real = (tb3 - infl).dropna() add(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)") add(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") add(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-)") add(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") add(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)") last_tips = int(tips.index.max()) - 1 add(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)") add(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") tips_d = fred("DFII10") add("tips_10y_real_yield_latest", tips_d.iloc[-1], "percent", str(tips_d.index[-1].date()), "FRED:DFII10", "latest daily observation") # ---------------------------------------------------------------------------- # 5. Credit spreads, loan standards, bank credit, bank leverage # ---------------------------------------------------------------------------- for sid, lab in [("BAMLC0A0CM", "ig_oas"), ("BAMLH0A0HYM2", "hy_oas")]: s = fred(sid) per = f"{s.index[0].date()}-{s.index[-1].date()}" 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)") add(f"{lab}_max_window", s.max(), "percent", f"{s.idxmax().date()}", f"FRED:{sid}", "max within the window") add(f"{lab}_latest", s.iloc[-1], "percent", str(s.index[-1].date()), f"FRED:{sid}", "latest") 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") 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") sl = fred("DRTSCILM") add("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") add("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") add("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") bc = annual_mean(fred("TOTBKCR")); bcg = np.log(bc).diff().dropna() add(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)") add(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") eq = fred("RALACBW027SBOG"); ta = fred("TLAACBW027SBOG") lev = (eq / ta).dropna(); lev_a = annual_mean(lev) add("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)") add(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") add(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") add("bank_assets_to_equity_leverage_latest", 1 / lev.iloc[-1], "x", str(lev.index[-1].date()), "same", "inverse of equity/assets") # ---------------------------------------------------------------------------- # 6. Wealth distribution (Fed DFA), latest quarter # ---------------------------------------------------------------------------- dfa = pd.read_csv(RAW / "dfa" / "dfa-networth-shares.csv") latest = dfa["Date"].iloc[-1] d = dfa[dfa["Date"] == latest].set_index("Category")["Net worth"] top1 = d["TopPt1"] + d["RemainingTop1"]; top10 = top1 + d["Next9"]; bot50 = d["Bottom50"] add("wealth_share_top0p1_dfa_latest", d["TopPt1"], "percent of household net worth", latest, "DFA dfa-networth-shares.csv", "Fed Distributional Financial Accounts, latest quarter") add("wealth_share_top1_dfa_latest", top1, "percent of household net worth", latest, "DFA dfa-networth-shares.csv", "TopPt1 + RemainingTop1") add("wealth_share_top10_dfa_latest", top10, "percent of household net worth", latest, "DFA dfa-networth-shares.csv", "top1 + Next9") add("wealth_share_next40_dfa_latest", d["Next40"], "percent of household net worth", latest, "DFA dfa-networth-shares.csv", "50th-90th percentile") add("wealth_share_bottom50_dfa_latest", bot50, "percent of household net worth", latest, "DFA dfa-networth-shares.csv", "bottom half") first = dfa["Date"].iloc[0]; d0 = dfa[dfa["Date"] == first].set_index("Category")["Net worth"] add("wealth_share_top1_dfa_1989Q3", d0["TopPt1"] + d0["RemainingTop1"], "percent", first, "DFA dfa-networth-shares.csv", "first quarter available, for comparison") add("wealth_share_bottom50_dfa_1989Q3", d0["Bottom50"], "percent", first, "DFA dfa-networth-shares.csv", "first quarter available") # grouped-data Gini (lower bound: assumes equality within each of the 5 groups) p = np.array([0.50, 0.40, 0.09, 0.009, 0.001]) w = np.array([bot50, d["Next40"], d["Next9"], d["RemainingTop1"], d["TopPt1"]]) / 100.0 L = np.concatenate([[0.0], np.cumsum(w)]) gini_lb = 1.0 - float(np.sum((L[:-1] + L[1:]) * p)) add("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)") # ---------------------------------------------------------------------------- # 7. Equity returns, premium and volatility (Shiller monthly data) # ---------------------------------------------------------------------------- sh = pd.read_excel(RAW / "shiller" / "ie_data.xls", sheet_name="Data", header=None) sh = sh.iloc[8:].copy() sh = sh[pd.to_numeric(sh[0], errors="coerce").notna()] sh = sh.rename(columns={0: "date", 1: "P", 2: "D", 4: "CPI", 6: "GS10"})[["date", "P", "D", "CPI", "GS10"]].astype(float) sh = sh.reset_index(drop=True) assert abs(sh.loc[0, "date"] - 1871.01) < 1e-6 sh["month"] = pd.period_range("1871-01", periods=len(sh), freq="M") sh["year"] = sh["month"].dt.year sh["ret_nom"] = (sh["P"] + sh["D"].ffill() / 12.0) / sh["P"].shift(1) - 1.0 ann = sh.groupby("year").agg(ret=("ret_nom", lambda r: np.prod(1 + r.values) - 1), n=("ret_nom", "size"), cpi_end=("CPI", "last"), gs10=("GS10", "mean")) ann = ann[ann["n"] == 12] ann["infl"] = ann["cpi_end"] / ann["cpi_end"].shift(1) - 1.0 ann["ret_real"] = (1 + ann["ret"]) / (1 + ann["infl"]) - 1.0 ann = ann.dropna(subset=["ret_real"]) last_sh = int(ann.index.max()) for (y0, y1) in [(1872, last_sh), (1926, last_sh), (1984, last_sh)]: a = win(ann, y0, y1) 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") 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") 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)") tb3_full = annual_mean(fred("TB3MS")) joint = ann.join(tb3_full.rename("tb3"), how="inner") for (y0, y1) in [(1934, last_sh), (1948, last_sh), (1984, last_sh)]: j = win(joint, y0, y1) 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") 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") 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") sp = fred("SP500"); sp_a = sp.groupby(sp.index.year).last(); sp_r = np.log(sp_a).diff().dropna() sp_r = sp_r[sp_r.index <= last_full] add(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") add(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") # FRED:BOGZ1FL073164003Q (Z.1 "NYSE Composite Index, Level") is downloaded but NOT used: its values # (e.g. 2,383,423 in 2026Q2) are not an index level in the usual sense and the FRED unit label is "Percent"; # treating it as a price index gave an implausible 30% annual sd. Left out until the concept is clarified. # ---------------------------------------------------------------------------- # write # ---------------------------------------------------------------------------- df = pd.DataFrame(rows, columns=["target", "value", "units", "period", "source_series", "method"]) df.to_csv(OUT, index=False, float_format="%.6g") print(f"wrote {OUT} ({len(df)} rows; {df['value'].isna().sum()} NOT OBTAINED/NOT COMPUTED)") with pd.option_context("display.max_rows", 500, "display.width", 250, "display.max_colwidth", 60): print(df[["target", "value", "units", "period"]].to_string(index=False))