1"""F1-conf-01: 12-month CPI inflation from two ALFRED vintages.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json4from pathlib import Path5import pandas as pd67D, O = Path("data"), Path("out")8O.mkdir(exist_ok=True)9F_OLD, F_NEW, F_ASM = "cpi_vintage_2016-07-15.csv", "cpi_vintage_2026-10-28.csv", "cpi_assembled.csv"10old = pd.read_csv(D / F_OLD, dtype={"date": str})11new = pd.read_csv(D / F_NEW, dtype={"date": str})12asm = pd.read_csv(D / F_ASM, dtype={"date": str, "vintage": str})13for df in (old, new, asm):14 df["date"] = pd.to_datetime(df["date"])15issues = []1617# --- diagnose the hand-assembled file (used only for checking, not for estimation)18exp = pd.concat([old[old.date < "2014-01-01"].assign(vintage="2016-07-15"),19 new.assign(vintage="2026-10-28")], ignore_index=True)20full = pd.DataFrame({"date": pd.date_range("1947-01-01", "2026-09-01", freq="MS")})21missing_dates = sorted(set(full.date) - set(asm.date))22if missing_dates:23 issues.append({"issue": f"cpi_assembled.csv has no row for {', '.join(d.strftime('%Y-%m') for d in missing_dates)} "24 f"and contains a row dated {asm.date.max().strftime('%Y-%m')} (outside the vintage coverage); "25 "from 2016-08 onward every value is shifted one month later than its true date "26 "(e.g. the 2016-07 value 240.101 is labelled 2016-08).",27 "handling": "cpi_assembled.csv not used for estimation; series rebuilt directly from the two authoritative vintage files."})28ov = asm[(asm.date >= "2014-01-01") & (asm.date <= "2016-06-01")]29if (ov.vintage == "2016-07-15").any():30 issues.append({"issue": f"cpi_assembled.csv uses the older 2016-07-15 vintage for {int((ov.vintage=='2016-07-15').sum())} overlap months "31 "2014-01..2016-06, although the 2026-10-28 vintage covers them (values differ by seasonal-factor revisions).",32 "handling": "Per the availability rule, the 2026-10-28 vintage value is used for every month it covers."})33nnew_na = new[new.cpi.isna()]34for d in nnew_na.date:35 issues.append({"issue": f"2026-10-28 vintage reports no value for {d.strftime('%Y-%m')} (empty cpi cell); "36 "the 2016-07-15 vintage does not cover that month either.",37 "handling": f"No interpolation; g_t for {d.strftime('%Y-%m')} (and for {(d+pd.DateOffset(months=12)).strftime('%Y-%m')}, outside the window) is not computable and is counted as missing."})38asm_na = asm[asm.cpi.isna()]39if len(asm_na):40 issues.append({"issue": f"cpi_assembled.csv has an empty cpi at row labelled {asm_na.date.iloc[0].strftime('%Y-%m')} (this is the 2025-10 gap, mislabelled by the one-month shift).",41 "handling": "Not used; see shift issue."})42issues.append({"issue": "Growth rates for 2014-01..2014-12 combine a 2026-10-28-vintage numerator with a 2016-07-15-vintage lag (2013 months not in the recent extract); seasonal factors for 2013 were later revised.",43 "handling": "Applied the stated rule (most recent vintage available per month) as specified; flagged as a limitation in the report."})4445# --- build the series per availability rule46new_av = new.dropna(subset=["cpi"])47ser = pd.concat([new_av.assign(vintage="2026-10-28"),48 old[~old.date.isin(new.date)].assign(vintage="2016-07-15")])49ser = full.merge(ser, on="date", how="left")50s = ser.set_index("date")51lag = s.cpi.shift(12) # index is complete monthly, so shift(12) = 12 months52lag_v = s.vintage.shift(12)53s["cpi_lag12"] = lag54s["growth_12m"] = 100 * (s.cpi / s.cpi_lag12 - 1)55s["vintage_used"] = s.vintage.where(lag_v == s.vintage, s.vintage + "/lag:" + lag_v.fillna(""))56win = s.loc["2001-01-01":"2026-09-01"].reset_index()57n_win = len(win)58ok = win.dropna(subset=["growth_12m"])59an = ok[["date", "cpi", "cpi_lag12", "vintage_used", "growth_12m"]].copy()60an["date"] = an.date.dt.strftime("%Y-%m-01")61an.to_csv(O / "analytic.csv", index=False, float_format="%.10g")6263excl = pd.DataFrame([64 dict(step=1, rule="All calendar months in window 2001-01..2026-09", n_before=n_win, n_after=n_win, n_removed=0,65 rationale="Starting set: one row per month of the estimation window"),66 dict(step=2, rule="Drop months whose g_t cannot be computed (P_t or P_{t-12} absent in all vintages)",67 n_before=n_win, n_after=len(ok), n_removed=n_win - len(ok),68 rationale="No interpolation/imputation allowed; missing: " + ", ".join(win[win.growth_12m.isna()].date.dt.strftime("%Y-%m"))),69])70excl.to_csv(O / "exclusions.csv", index=False)7172g = ok.set_index("date").growth_12m73def gm(m):74 v = s.growth_12m.get(pd.Timestamp(m + "-01"))75 return None if pd.isna(v) else round(float(v), 6)76last_level = ser.dropna(subset=["cpi"]).date.max().strftime("%Y-%m")77nmiss = n_win - len(ok)78res = {79 "task_id": "F1-conf-01",80 "first_growth_month": g.index.min().strftime("%Y-%m"),81 "last_month": last_level,82 "n_months_growth": int(len(ok)),83 "n_months_missing_growth": int(nmiss),84 "n_analysed": int(len(ok)),85 "mean_growth_12m": round(float(g.mean()), 6),86 "mean_growth_12m_2001_2013": round(float(g["2001-01-01":"2013-12-01"].mean()), 6),87 "mean_growth_12m_2014_2026": round(float(g["2014-01-01":"2026-09-01"].mean()), 6),88 "growth_2014_06": gm("2014-06"), "growth_2016_07": gm("2016-07"), "growth_2020_05": gm("2020-05"),89 "growth_2022_06": gm("2022-06"), "growth_2026_09": gm("2026-09"),90 "max_growth_12m": round(float(g.max()), 6), "max_growth_month": g.idxmax().strftime("%Y-%m"),91 "min_growth_12m": round(float(g.min()), 6), "min_growth_month": g.idxmin().strftime("%Y-%m"),92 "feasibility": "complete" if nmiss == 0 else "partial",93 "feasibility_reason": None if nmiss == 0 else94 f"{nmiss} month(s) of the window have no computable g_t because the source publishes no index level for 2025-10; "95 "summary statistics are over the remaining months.",96 "data_issues": issues,97 "n_input_rows": {F_OLD: int(len(old)), F_NEW: int(len(new)), F_ASM: int(len(asm))},98 "random_seed": None,99}100(O / "results.json").write_text(json.dumps(res, indent=2))101print(json.dumps(res, indent=2))102103# --- sensitivity (not part of results.json): 2014 growth using the 2016-07-15 vintage for both P_t and P_{t-12}104o = old.set_index("date").cpi105g_old = (100 * (o / o.shift(12) - 1)).loc["2014-01-01":"2014-12-01"]106g_mix = s.growth_12m.loc["2014-01-01":"2014-12-01"]107sens = {"mean_2014_mixed_vintage": round(float(g_mix.mean()), 6),108 "mean_2014_old_vintage_only": round(float(g_old.mean()), 6),109 "max_abs_diff_2014_months": round(float((g_mix - g_old).abs().max()), 6),110 "growth_2014_06_old_vintage_only": round(float(g_old.loc["2014-06-01"]), 6),111 "mean_growth_12m_2001_2013_n": int(g["2001-01-01":"2013-12-01"].size),112 "mean_growth_12m_2014_2026_n": int(g["2014-01-01":"2026-09-01"].size),113 "sd_growth_12m": round(float(g.std()), 6)}114(O / "sensitivity.json").write_text(json.dumps(sens, indent=2))115