1"""Replication entry point for F6-conf-02. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads data/ only; writes out/results.json, out/figure1.png, out/table1.csv, out/exclusions.csv.4No randomness is used (random_seed = null)."""5import json, hashlib, os6import 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)14EXPECTED_SHA = "2cdda6ccb2147c8861997166cb0d32cf85ff2744e1e7c0f47eded0f354c2db51"15sha = hashlib.sha256(open(DATA, "rb").read()).hexdigest()1617df = pd.read_csv(DATA)18n_raw = len(df)19df["dt"] = pd.to_datetime(df["date"], format="%Y-%m-%d", errors="coerce")20issues = []21if sha != EXPECTED_SHA:22 issues.append({"issue": f"SHA-256 {sha} differs from README", "handling": "reported; analysis continued"})23n_bad_date = int(df["dt"].isna().sum()); n_missing_y = int(df["yield_5y"].isna().sum())24n_dup = int(df["date"].duplicated().sum()); monotone = bool(df["dt"].is_monotonic_increasing)25for cond, txt in [(n_bad_date, f"{n_bad_date} unparseable dates"), (n_missing_y, f"{n_missing_y} missing yields"),26 (n_dup, f"{n_dup} duplicate dates"), (not monotone, "dates not sorted ascending")]:27 if cond:28 issues.append({"issue": txt, "handling": "reported; rows kept as provided"})2930# One-day changes between consecutive rows, in bp, dated t (first row has none)31df["chg_bp"] = 100 * df["yield_5y"].diff()32excl = [dict(step=1, rule="drop first row (no previous row for one-day change)", n_before=n_raw,33 n_after=n_raw - 1, n_removed=1, rationale="change defined between consecutive rows")]34chg = df["chg_bp"].dropna()35baseline = float(chg.mean())3637n_treated = int((df["post_cut"] == 1).sum())38inconsistent = int(((df["dt"] >= "2020-03-04") != (df["post_cut"] == 1)).sum())39if inconsistent:40 issues.append({"issue": f"{inconsistent} rows where post_cut disagrees with date>=2020-03-04",41 "handling": "post_cut used as provided (task rule); reported"})4243effect = lo = hi = rmse = None; n_analysed = None44if n_treated == 0:45 issues.append({"issue": f"File covers only {df['date'].iloc[0]} to {df['date'].iloc[-1]} (README retrieval "46 "end_date=2019-12-31); no row has post_cut=1 and no observation exists on/after the "47 "2020-03-04 cut, so the post-cut window is empty",48 "handling": "estimand not computable; effect, CI, model_metric and n_analysed set to null; "49 "conclusion premise_invalid; no proxy window substituted"})50 feas, reason, concl = ("infeasible",51 "The data contain zero rows with post_cut=1 (sample ends 2019-12-31, before the 2020-03-04 cut), so the "52 "post-cut mean change, the 120-day pre-window anchored at the first flagged day, the effect and its "53 "interval cannot be computed.", "premise_invalid")54 excl.append(dict(step=2, rule="keep post_cut=1 rows plus the 120 rows preceding the first flagged day",55 n_before=n_raw - 1, n_after=0, n_removed=n_raw - 1,56 rationale="estimation windows; no flagged day exists so both windows are empty"))57else:58 first = int(np.flatnonzero(df["post_cut"].values == 1)[0])59 post = df.loc[df["post_cut"] == 1, "chg_bp"].dropna()60 pre = df["chg_bp"].iloc[max(first - 120, 0):first].dropna()61 if len(pre) < 120:62 issues.append({"issue": f"pre-window has only {len(pre)} changes", "handling": "used available"})63 effect = float(post.mean() - pre.mean())64 v1, v0 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)65 se = np.sqrt(v1 + v0)66 dof = (v1 + v0) ** 2 / (v1 ** 2 / (len(post) - 1) + v0 ** 2 / (len(pre) - 1))67 t = stats.t.ppf(0.975, dof); lo, hi = float(effect - t * se), float(effect + t * se)68 resid = np.concatenate([post - post.mean(), pre - pre.mean()])69 rmse = float(np.sqrt(np.mean(resid ** 2))); n_analysed = int(len(post) + len(pre))70 feas, reason = "complete", None71 concl = "positive_effect" if lo > 0 else "negative_effect" if hi < 0 else "no_detectable_effect"72 excl.append(dict(step=2, rule="keep estimation windows", n_before=n_raw - 1, n_after=n_analysed,73 n_removed=n_raw - 1 - n_analysed, rationale="pre/post windows"))74pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False)7576# Table 177df["year"] = df["dt"].dt.year78t1 = df.groupby("year").agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),79 mean_change_bp=("chg_bp", "mean")).reset_index().round(4)80t1.to_csv(f"{OUT}/table1.csv", index=False)8182# Figure 183fig, ax = plt.subplots(figsize=(9, 4))84ax.plot(df["dt"], df["yield_5y"], lw=0.8)85ax.axvline(pd.Timestamp("2020-03-04"), color="red", ls="--", lw=1, label="2020-03-04 cut (outside sample)")86ax.set_xlabel("date"); ax.set_ylabel("GoC 5-year yield (%)")87ax.set_title(f"BD.CDN.5YR.DQ.YLD, {df['date'].iloc[0]} to {df['date'].iloc[-1]}")88ax.legend(); fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)8990res = {"task_id": "F6-conf-02", "sample_start": df["date"].iloc[0], "sample_end": df["date"].iloc[-1],91 "n_treated": n_treated, "baseline_metric": round(baseline, 6), "model_metric": rmse,92 "effect_estimate": effect, "effect_ci_low": lo, "effect_ci_high": hi, "conclusion": concl,93 "feasibility": feas, "feasibility_reason": reason, "data_issues": issues,94 "n_input_rows": {"gc5y_daily.csv": n_raw}, "n_analysed": n_analysed, "random_seed": None,95 "baseline_n_changes": int(len(chg)), "baseline_sd_change_bp": round(float(chg.std(ddof=1)), 6),96 "input_sha256": sha}97json.dump(res, open(f"{OUT}/results.json", "w"), indent=2)98print(f"conclusion={concl}; n_treated={n_treated}; baseline mean change={baseline:.4f} bp over {len(chg)} changes")99