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"13OUT = "out"14PRE_N = 12015os.makedirs(OUT, exist_ok=True)1617raw = pd.read_csv(DATA, dtype={"date": str})18n_raw = len(raw)19sha = hashlib.sha256(open(DATA, "rb").read()).hexdigest()20issues = []2122# ---------------- audit ----------------23df = raw.copy()24df["date_parsed"] = pd.to_datetime(df["date"], format="%Y-%m-%d", errors="coerce")25n_bad_date = int(df["date_parsed"].isna().sum())26df["yield_5y"] = pd.to_numeric(df["yield_5y"], errors="coerce")27n_na_y = int(df["yield_5y"].isna().sum())28n_na_t = int(df["post_cut"].isna().sum())29n_dup = int(df["date"].duplicated().sum())30monotone = bool(df["date_parsed"].is_monotonic_increasing)31if n_bad_date or n_na_y or n_na_t:32 issues.append({"issue": f"{n_bad_date} unparseable dates, {n_na_y} missing yields, {n_na_t} missing post_cut",33 "handling": "reported; changes computed between consecutive rows as given"})34if n_dup:35 issues.append({"issue": f"{n_dup} duplicate dates", "handling": "reported; not removed"})36if not monotone:37 issues.append({"issue": "dates not sorted ascending", "handling": "file order kept as required by the change definition; reported"})3839sample_start = str(raw["date"].iloc[0]); sample_end = str(raw["date"].iloc[-1])40treated = df["post_cut"] == 141n_treated = int(treated.sum())42bad_flag = int(((df["date_parsed"] >= "2020-03-04") != treated).sum())43if bad_flag:44 issues.append({"issue": f"{bad_flag} rows where post_cut disagrees with the documented rule (date >= 2020-03-04)",45 "handling": "post_cut used exactly as provided, per task rules"})4647# no rows are removed for the descriptive baseline48excl = [{"step": "S0", "rule": "read file; no filtering applied (no missing/duplicate/unparseable rows found)",49 "n_before": n_raw, "n_after": n_raw, "n_removed": 0,50 "rationale": "audit found no rows requiring removal"}]5152# ---------------- changes and baseline ----------------53df["d_bp"] = 100.0 * df["yield_5y"].diff()54baseline_metric = float(df["d_bp"].mean()) # mean over n_raw-1 one-day changes5556# ---------------- estimand ----------------57model_metric = effect = lo = hi = None58n_analysed = None59if n_treated == 0:60 feasibility = "infeasible"; conclusion = "premise_invalid"61 feasibility_reason = (f"The file covers {sample_start} to {sample_end} and contains no rows with post_cut = 1 "62 "(the cut was announced 2020-03-04, after the last date). The post-cut mean change, "63 "and hence the effect, its interval and the two-window model RMSE, cannot be computed. "64 "No proxy estimate is substituted.")65 issues.append({"issue": f"Sample ends {sample_end}, before the 2020-03-04 cut; n_treated = 0 (all {n_raw} rows have post_cut = 0). "66 "The README retrieval URL requests end_date=2019-12-31, consistent with this truncation.",67 "handling": "Estimand declared not estimable: effect_estimate, CI, model_metric, n_analysed set to null; conclusion premise_invalid."})68 excl.append({"step": "S1", "rule": "keep rows in the two estimation windows (120 pre-window days + post_cut=1 days)",69 "n_before": n_raw, "n_after": 0, "n_removed": n_raw,70 "rationale": "no post_cut=1 rows exist, so the windows are undefined and no row enters the effect estimate"})71else:72 first = int(np.flatnonzero(treated.values)[0])73 post = df.loc[treated, "d_bp"].dropna()74 pre = df["d_bp"].iloc[max(first - PRE_N, 0):first].dropna()75 if len(pre) < PRE_N:76 issues.append({"issue": f"only {len(pre)} pre-window changes available (<{PRE_N})", "handling": "used what is available; reported"})77 if len(pre) >= 2 and len(post) >= 2:78 effect = float(post.mean() - pre.mean())79 v1, v2 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)80 se = np.sqrt(v1 + v2)81 dof = (v1 + v2) ** 2 / (v1 ** 2 / (len(post) - 1) + v2 ** 2 / (len(pre) - 1))82 q = stats.t.ppf(0.975, dof)83 lo, hi = float(effect - q * se), float(effect + q * se)84 resid = np.concatenate([post - post.mean(), pre - pre.mean()])85 model_metric = float(np.sqrt(np.mean(resid ** 2)))86 n_analysed = int(len(pre) + len(post))87 conclusion = "positive_effect" if lo > 0 else "negative_effect" if hi < 0 else "no_detectable_effect"88 feasibility, feasibility_reason = "complete", None89 else:90 feasibility, conclusion = "infeasible", "premise_invalid"91 feasibility_reason = "too few observations in one of the windows"92 excl.append({"step": "S1", "rule": "keep rows in the two estimation windows",93 "n_before": n_raw, "n_after": n_analysed or 0, "n_removed": n_raw - (n_analysed or 0),94 "rationale": "estimand is defined on these windows only"})9596pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False)9798# ---------------- table and figure ----------------99df["year"] = df["date_parsed"].dt.year100tab = df.groupby("year").agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),101 mean_change_bp=("d_bp", "mean")).reset_index()102tab["mean_yield"] = tab["mean_yield"].round(6); tab["mean_change_bp"] = tab["mean_change_bp"].round(6)103tab.to_csv(f"{OUT}/table1.csv", index=False)104105fig, ax = plt.subplots(figsize=(9, 4))106ax.plot(df["date_parsed"], df["yield_5y"], lw=0.8)107ax.set_title(f"GoC 5-year benchmark yield, {sample_start} to {sample_end}\n(BD.CDN.5YR.DQ.YLD; 2020-03-04 cut lies after sample end)")108ax.set_ylabel("percent"); ax.set_xlabel("date")109fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)110111res = {112 "task_id": "F6-conf-02",113 "sample_start": sample_start, "sample_end": sample_end,114 "n_treated": n_treated,115 "baseline_metric": baseline_metric,116 "model_metric": model_metric,117 "effect_estimate": effect, "effect_ci_low": lo, "effect_ci_high": hi,118 "conclusion": conclusion,119 "feasibility": feasibility, "feasibility_reason": feasibility_reason,120 "data_issues": issues,121 "n_input_rows": {"gc5y_daily.csv": n_raw},122 "n_analysed": n_analysed,123 "random_seed": None,124}125with open(f"{OUT}/results.json", "w") as fh:126 json.dump(res, fh, indent=2)127print(f"sha256={sha} rows={n_raw} {sample_start}..{sample_end} n_treated={n_treated}")128print(f"baseline_metric={baseline_metric:.6f} bp; conclusion={conclusion}; feasibility={feasibility}")129