1"""F1-conf-01: 12-month CPIAUCSL growth from two ALFRED vintages.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/; writes out/results.json, out/analytic.csv, out/exclusions.csv.4No randomness is used (random_seed = null)."""5import json, os6import numpy as np7import pandas as pd89D, O = "data", "out"10F_OLD, F_NEW, F_ASM = "cpi_vintage_2016-07-15.csv", "cpi_vintage_2026-10-28.csv", "cpi_assembled.csv"11W0, W1 = pd.Timestamp("2001-01-01"), pd.Timestamp("2026-09-01")1213old = pd.read_csv(os.path.join(D, F_OLD), dtype={"date": str})14new = pd.read_csv(os.path.join(D, F_NEW), dtype={"date": str})15asm = pd.read_csv(os.path.join(D, F_ASM), dtype={"date": str, "vintage": str}) # read for audit only16for d in (old, new, asm):17 d["date"] = pd.to_datetime(d["date"])18n_input_rows = {F_OLD: int(len(old)), F_NEW: int(len(new)), F_ASM: int(len(asm))}1920issues = []21# --- integrity checks on the authoritative files22for name, d in ((F_OLD, old), (F_NEW, new)):23 assert not d["date"].duplicated().any(), name24 full = pd.date_range(d["date"].min(), d["date"].max(), freq="MS")25 assert len(full) == len(d), f"{name} has gaps in date sequence"26miss_new = new.loc[new["cpi"].isna(), "date"].dt.strftime("%Y-%m").tolist()27miss_old = old.loc[old["cpi"].isna(), "date"].dt.strftime("%Y-%m").tolist()28if miss_new:29 issues.append({"issue": f"{F_NEW}: empty cpi value (source reports no value) for month(s) {miss_new}",30 "handling": "Not imputed or interpolated and not replaced by an older-vintage value "31 "(the 2026-10-28 vintage covers the month, and the older vintage has no value for it anyway); "32 "g_t for that month (and for t+12, which lies outside the window) is reported missing and counted."})33if miss_old:34 issues.append({"issue": f"{F_OLD}: empty cpi for {miss_old}", "handling": "treated as missing; no imputation"})3536# --- overlap / revisions between vintages37ov = old.merge(new, on="date", suffixes=("_old", "_new")).dropna()38diff = ov["cpi_new"] - ov["cpi_old"]39issues.append({"issue": f"Vintages overlap for {len(ov)} months ({ov.date.min():%Y-%m}..{ov.date.max():%Y-%m}) with revised values "40 f"(new minus old: min {diff.min():.3f}, max {diff.max():.3f}, {int((diff.abs() > 0).sum())} months differ; "41 "revisions between vintages, e.g. re-estimated seasonal factors)",42 "handling": "Per the availability rule, the 2026-10-28 value is used for every month it covers; the 2016-07-15 value only before 2014-01."})4344# --- audit of the hand-assembled file (not used for computation)45a = asm.merge(old.rename(columns={"cpi": "c_old"}), on="date", how="left").merge(46 new.rename(columns={"cpi": "c_new"}), on="date", how="left")47gaps = sorted(set(pd.date_range(asm.date.min(), asm.date.max(), freq="MS")) - set(asm.date))48if gaps:49 issues.append({"issue": f"{F_ASM}: month(s) absent from the sequence: {[g.strftime('%Y-%m') for g in gaps]}",50 "handling": "Assembled file not used for computation; series rebuilt from the two vintage files."})51wrong_v = a[(a.vintage == "2016-07-15") & a.c_new.notna()]52if len(wrong_v):53 issues.append({"issue": f"{F_ASM}: {len(wrong_v)} months ({wrong_v.date.min():%Y-%m}..{wrong_v.date.max():%Y-%m}) take values from the 2016-07-15 vintage although the 2026-10-28 vintage covers them (violates the most-recent-vintage rule)",54 "handling": "Assembled file not used; 2026-10-28 values used for those months."})55nv = a[a.vintage == "2026-10-28"]56mism = nv[~np.isclose(nv.cpi, nv.c_new, equal_nan=True)]57if len(mism):58 # test for a one-month shift: assembled value at t equals vintage value at t-159 sh = new.assign(date=new.date + pd.DateOffset(months=1)).rename(columns={"cpi": "c_prev"})60 m2 = mism.merge(sh, on="date", how="left")61 n_shift = int(np.isclose(m2.cpi, m2.c_prev, equal_nan=True).sum())62 issues.append({"issue": f"{F_ASM}: {len(mism)} rows labelled vintage 2026-10-28 ({mism.date.min():%Y-%m}..{mism.date.max():%Y-%m}) do not match that vintage; "63 f"{n_shift} of them equal the vintage value of the previous month (dates shifted forward by one month, creating a spurious {asm.date.max():%Y-%m} row and misplacing the missing value)",64 "handling": "Assembled file not used for computation; only the two authoritative vintage files are used."})6566# --- combined series under the most-recent-vintage rule67new_start, new_end = new.date.min(), new.date.max()68o2 = old[(old.date < new_start)].assign(vintage_used="2016-07-15")69n2 = new.assign(vintage_used="2026-10-28") # covers 2014-01..2026-09 incl. months with empty value70s = pd.concat([o2, n2]).sort_values("date").set_index("date")71assert s.index.is_unique72assert len(s) == len(pd.date_range(s.index.min(), s.index.max(), freq="MS"))73s["cpi_lag12"] = s["cpi"].shift(12)74assert (s.index[12:] - pd.DateOffset(months=12) == s.index[:-12]).all()75s["growth_12m"] = 100.0 * (s["cpi"] / s["cpi_lag12"] - 1.0)7677n_mixed = int(((s.vintage_used == "2026-10-28") & (s.vintage_used.shift(12) == "2016-07-15") & s.index.to_series().between(W0, W1)).sum())78issues.append({"issue": f"{n_mixed} window months (2014-01..2014-12) have P_t from the 2026-10-28 vintage but P_(t-12) from the 2016-07-15 vintage (the newer extract starts 2014-01), so their g_t mixes vintages",79 "handling": "Kept, as the TASK rule prescribes the most recent vintage per month; sensitivity value growth using the 2016 vintage only is reported as sens_growth_2014_06_old_vintage_only."})80_o = old.set_index("date")["cpi"]81sens_2014_06 = float(round(100 * (_o[pd.Timestamp("2014-06-01")] / _o[pd.Timestamp("2013-06-01")] - 1), 10))8283# --- exclusion ledger84led = []85def step(name, rule, before, after, why):86 led.append(dict(step=name, rule=rule, n_before=len(before), n_after=len(after),87 n_removed=len(before) - len(after), rationale=why))88s0 = s89s1 = s0[(s0.index >= W0) & (s0.index <= W1)]90step("1_window", "keep months 2001-01..2026-09 of combined monthly series", s0, s1,91 "estimation window defined in TASK.md (combined series spans 1947-01..2026-09)")92s2 = s1[s1.cpi.notna()]93step("2_missing_level", "drop months whose own index level P_t is missing in the source", s1, s2,94 "no imputation allowed; g_t not computable")95s3 = s2[s2.cpi_lag12.notna()]96step("3_missing_lag", "drop months whose lagged level P_{t-12} is missing in the source", s2, s3,97 "no imputation allowed; g_t not computable")98pd.DataFrame(led).to_csv(os.path.join(O, "exclusions.csv"), index=False)99100an = s3.reset_index()101out = an[["date", "cpi", "cpi_lag12", "vintage_used", "growth_12m"]].copy()102out["date"] = out.date.dt.strftime("%Y-%m-01")103out.to_csv(os.path.join(O, "analytic.csv"), index=False, float_format="%.10g")104105g = s3["growth_12m"]106def gm(m):107 t = pd.Timestamp(m + "-01")108 return float(round(g[t], 10)) if t in g.index else None109n_win = len(s1)110missing_months = s1.index[s1.growth_12m.isna()].strftime("%Y-%m").tolist()111last_level = s.index[s.cpi.notna()].max()112res = {113 "task_id": "F1-conf-01",114 "first_growth_month": g.index.min().strftime("%Y-%m"),115 "last_month": last_level.strftime("%Y-%m"),116 "n_months_growth": int(len(g)),117 "n_months_missing_growth": int(n_win - len(g)),118 "n_analysed": int(len(g)),119 "mean_growth_12m": float(round(g.mean(), 10)),120 "mean_growth_12m_2001_2013": float(round(g[g.index <= "2013-12-01"].mean(), 10)),121 "mean_growth_12m_2014_2026": float(round(g[g.index >= "2014-01-01"].mean(), 10)),122 "growth_2014_06": gm("2014-06"), "growth_2016_07": gm("2016-07"), "growth_2020_05": gm("2020-05"),123 "growth_2022_06": gm("2022-06"), "growth_2026_09": gm("2026-09"),124 "max_growth_12m": float(round(g.max(), 10)), "max_growth_month": g.idxmax().strftime("%Y-%m"),125 "min_growth_12m": float(round(g.min(), 10)), "min_growth_month": g.idxmin().strftime("%Y-%m"),126 "feasibility": "complete" if not missing_months else "partial",127 "feasibility_reason": None if not missing_months else128 f"g_t not computable for {missing_months} because the source (2026-10-28 vintage) reports no index level for 2025-10; "129 "no imputation allowed. All other window months computed.",130 "data_issues": issues,131 "n_input_rows": n_input_rows,132 "random_seed": None,133 # supplementary134 "missing_growth_months": missing_months,135 "n_window_months": int(n_win),136 "n_months_2001_2013": int((g.index <= "2013-12-01").sum()),137 "n_months_2014_2026": int((g.index >= "2014-01-01").sum()),138 "n_months_mixed_vintage_growth": n_mixed,139 "sens_growth_2014_06_old_vintage_only": sens_2014_06,140 "sd_growth_12m": float(round(g.std(ddof=1), 10)),141}142with open(os.path.join(O, "results.json"), "w") as f:143 json.dump(res, f, indent=2)144print(json.dumps({k: v for k, v in res.items() if k != "data_issues"}, indent=1))145