1"""Replication entry point for F6-conf-02. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads only data/gc5y_daily.csv; writes out/results.json, out/figure1.png,4out/table1.csv, out/exclusions.csv. No randomness is used."""5import json, os, hashlib6import numpy as np, pandas as pd7import matplotlib8matplotlib.use("Agg")9import matplotlib.pyplot as plt10from scipy import stats1112DATA = "data/gc5y_daily.csv"; OUT = "out"13os.makedirs(OUT, exist_ok=True)14raw = pd.read_csv(DATA)15n_raw = len(raw)16sha = hashlib.sha256(open(DATA, "rb").read()).hexdigest()17d = raw.copy(); d["date"] = pd.to_datetime(d["date"])1819# ---- audit ----20issues = []21n_treated = int((d["post_cut"] == 1).sum())22n_na = int(d["yield_5y"].isna().sum()); n_dup = int(d["date"].duplicated().sum())23sorted_ok = bool(d["date"].is_monotonic_increasing)24bad_flag = int(((d["date"] >= "2020-03-04").astype(int) != d["post_cut"]).sum())25if n_treated == 0:26 issues.append({"issue": f"No rows with post_cut = 1: the file covers {raw['date'].iloc[0]} to {raw['date'].iloc[-1]} "27 "(README retrieval used end_date=2019-12-31), entirely before the 2020-03-04 cut, so the post-cut window is empty.",28 "handling": "Estimand declared not computable (premise_invalid); effect, CI, model_metric and n_analysed set to null. No proxy or substitute data used."})29if n_na: issues.append({"issue": f"{n_na} missing yields", "handling": "Changes involving missing values dropped."})30if n_dup: issues.append({"issue": f"{n_dup} duplicate dates", "handling": "Reported; not removed."})31if not sorted_ok: issues.append({"issue": "Dates not sorted", "handling": "Sorted by date before differencing."})32if bad_flag: issues.append({"issue": f"{bad_flag} rows where post_cut disagrees with date>=2020-03-04", "handling": "post_cut used as provided per task rules."})33gaps = d["date"].diff().dt.days34issues.append({"issue": f"Calendar gaps between consecutive rows up to {int(gaps.max())} days (weekends/holidays; max gaps around year-end).",35 "handling": "Per task rules, one-day change is defined between consecutive rows regardless of calendar gap; no adjustment."})36d = d.sort_values("date").reset_index(drop=True)3738# ---- changes and baseline ----39d["chg_bp"] = 100 * d["yield_5y"].diff()40excl = [("1", "drop first row (no prior row for one-day change)", n_raw, n_raw - 1, 1,41 "Change undefined for first observation")]42chg = d.dropna(subset=["chg_bp"])43baseline = float(chg["chg_bp"].mean())4445# ---- main estimate ----46effect = lo = hi = rmse = None; n_an = None47if n_treated > 0:48 first = d.index[d["post_cut"] == 1][0]49 post = d.loc[d["post_cut"] == 1, "chg_bp"].dropna()50 pre = d.loc[max(first - 120, 1):first - 1, "chg_bp"].dropna()51 if len(pre) >= 2 and len(post) >= 2:52 effect = float(post.mean() - pre.mean())53 v1, v0 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)54 se = np.sqrt(v1 + v0)55 df = (v1 + v0) ** 2 / (v1 ** 2 / (len(post) - 1) + v0 ** 2 / (len(pre) - 1))56 t = stats.t.ppf(0.975, df); lo, hi = effect - t * se, effect + t * se57 resid = np.concatenate([post - post.mean(), pre - pre.mean()])58 rmse = float(np.sqrt(np.mean(resid ** 2))); n_an = int(len(post) + len(pre))59if effect is None:60 concl, feas = "premise_invalid", "infeasible"61 reason = ("The file contains no rows with post_cut = 1 (sample ends 2019-12-31, before the 2020-03-04 cut), "62 "so neither the post-cut window nor the '120 trading days before the first flagged day' window exists; "63 "the effect, its interval and the two-window RMSE cannot be computed.")64else:65 feas, reason = "complete", None66 concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")67excl.append(("2", "main-estimate windows (post_cut=1 rows + 120 preceding rows)", n_raw - 1,68 0 if n_an is None else n_an, (n_raw - 1) - (0 if n_an is None else n_an),69 "No post_cut=1 rows exist; window empty" if n_an is None else "Restrict to estimand windows"))70pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(f"{OUT}/exclusions.csv", index=False)7172# ---- table and figure ----73d["year"] = d["date"].dt.year74tab = d.groupby("year").agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),75 mean_change_bp=("chg_bp", "mean")).reset_index().round(4)76tab.to_csv(f"{OUT}/table1.csv", index=False)77fig, ax = plt.subplots(figsize=(9, 4))78ax.plot(d["date"], d["yield_5y"], lw=0.8)79ax.set_title("Government of Canada 5-year benchmark yield (BD.CDN.5YR.DQ.YLD)\n"80 "No post-cut (2020-03-04+) observations in file")81ax.set_ylabel("Yield (%)"); ax.set_xlabel("Date"); fig.tight_layout()82fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)8384r = lambda x: None if x is None else round(float(x), 6)85res = {"task_id": "F6-conf-02", "sample_start": raw["date"].iloc[0], "sample_end": raw["date"].iloc[-1],86 "n_treated": n_treated, "baseline_metric": r(baseline), "model_metric": r(rmse),87 "effect_estimate": r(effect), "effect_ci_low": r(lo), "effect_ci_high": r(hi),88 "conclusion": concl, "feasibility": feas, "feasibility_reason": reason,89 "data_issues": issues, "n_input_rows": {"gc5y_daily.csv": int(n_raw)},90 "n_analysed": n_an, "random_seed": None,91 "audit_n_changes": int(len(chg)), "audit_sd_change_bp": r(chg["chg_bp"].std(ddof=1)),92 "audit_input_sha256": sha}93json.dump(res, open(f"{OUT}/results.json", "w"), indent=2)94print(f"rows={n_raw} treated={n_treated} baseline={baseline:.4f} bp conclusion={concl}")95