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. Deterministic (no randomness used)."""5import json, os6import numpy as np7import pandas as pd8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt11from scipy import stats1213OUT = "out"14FN = "gc5y_daily.csv"15PRE_WINDOW = 12016os.makedirs(OUT, exist_ok=True)1718raw = pd.read_csv(os.path.join("data", FN), dtype=str)19n_raw = len(raw)20df = raw.copy()21df["date"] = pd.to_datetime(df["date"], format="%Y-%m-%d")22df["yield_5y"] = pd.to_numeric(df["yield_5y"], errors="coerce")23df["post_cut"] = pd.to_numeric(df["post_cut"], errors="coerce").astype("Int64")2425issues = []26# --- audit ---27n_missing = int(df[["date", "yield_5y", "post_cut"]].isna().any(axis=1).sum())28n_dup = int(df["date"].duplicated().sum())29monotonic = bool(df["date"].is_monotonic_increasing)30n_treated = int((df["post_cut"] == 1).sum())31gaps = df["date"].diff().dt.days32long_gaps = df.loc[gaps > 4, "date"].dt.strftime("%Y-%m-%d").tolist()33if n_missing:34 issues.append({"issue": f"{n_missing} rows with missing/non-numeric values",35 "handling": "reported; differences across them left NaN"})36if n_dup:37 issues.append({"issue": f"{n_dup} duplicated dates", "handling": "reported, not removed"})38if not monotonic:39 issues.append({"issue": "dates not sorted ascending", "handling": "file order kept per availability rule"})40sample_start = df["date"].min().strftime("%Y-%m-%d")41sample_end = df["date"].max().strftime("%Y-%m-%d")42if n_treated == 0:43 issues.append({44 "issue": (f"Coverage ends {sample_end} (README: Valet query end_date=2019-12-31), "45 "before the 2020-03-04 cut; post_cut is 0 for all "46 f"{n_raw} rows, so there are zero treated days"),47 "handling": ("Estimand not computable; effect, CI and model RMSE set to null, "48 "conclusion premise_invalid. post_cut not redefined and no proxy "49 "event or extrapolation substituted.")})50if long_gaps:51 issues.append({52 "issue": (f"{len(long_gaps)} calendar gaps of 5 days (year-end holidays, e.g. "53 f"{long_gaps[0]}); one-day changes span non-trading days"),54 "handling": "kept; changes defined between consecutive rows per the task rule"})5556# --- one-day changes (bp) between consecutive rows, dated t ---57df["chg_bp"] = 100.0 * df["yield_5y"].diff()58n_chg = int(df["chg_bp"].notna().sum())59baseline_metric = float(df["chg_bp"].mean())6061excl = [{"step": 1, "rule": "first row has no previous row for a one-day change",62 "n_before": n_raw, "n_after": n_chg, "n_removed": n_raw - n_chg,63 "rationale": "change defined as 100*(y_t - y_{t-1}) between consecutive rows"}]6465# --- effect: post-cut mean change minus mean over 120 days before first treated day ---66effect = lo = hi = rmse = None67n_analysed = None68treated_idx = np.flatnonzero((df["post_cut"] == 1).to_numpy())69if len(treated_idx) > 0:70 first = treated_idx[0]71 post = df.loc[df["post_cut"] == 1, "chg_bp"].dropna()72 pre = df["chg_bp"].iloc[max(0, first - PRE_WINDOW):first].dropna()73 if len(post) >= 2 and len(pre) >= 2:74 effect = float(post.mean() - pre.mean())75 v1, v0 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)76 se = np.sqrt(v1 + v0)77 dof = (v1 + v0) ** 2 / (v1 ** 2 / (len(post) - 1) + v0 ** 2 / (len(pre) - 1))78 q = stats.t.ppf(0.975, dof)79 lo, hi = float(effect - q * se), float(effect + q * se)80 resid = np.concatenate([post - post.mean(), pre - pre.mean()])81 rmse = float(np.sqrt(np.mean(resid ** 2)))82 n_analysed = int(len(post) + len(pre))83 n_win = 0 if n_analysed is None else n_analysed84else:85 n_win = 086excl.append({"step": 2, "rule": "keep post_cut=1 days and the 120 rows preceding the first post_cut=1 day",87 "n_before": n_chg, "n_after": n_win, "n_removed": n_chg - n_win,88 "rationale": ("estimand windows; no post_cut=1 rows exist so no window can be formed"89 if n_win == 0 else "estimand windows")})90pd.DataFrame(excl).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)9192if effect is None:93 conclusion, feas = "premise_invalid", "infeasible"94 reason = (f"The file covers {sample_start} to {sample_end} and contains no rows with post_cut=1 "95 "(the cut was announced 2020-03-04), so the post-cut mean change and the "96 "120-day pre-window anchored on the first treated day do not exist.")97else:98 feas, reason = "complete", None99 conclusion = ("positive_effect" if lo > 0 else "negative_effect" if hi < 0100 else "no_detectable_effect")101102# --- table1: per calendar year ---103df["year"] = df["date"].dt.year104t1 = (df.groupby("year")105 .agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),106 mean_change_bp=("chg_bp", "mean")).reset_index())107t1["mean_yield"] = t1["mean_yield"].round(6)108t1["mean_change_bp"] = t1["mean_change_bp"].round(6)109t1.to_csv(os.path.join(OUT, "table1.csv"), index=False)110111# --- figure1 ---112fig, ax = plt.subplots(figsize=(9, 4))113ax.plot(df["date"], df["yield_5y"], lw=0.8, color="C0")114ax.axvline(pd.Timestamp("2020-03-04"), color="red", ls="--", lw=1,115 label="2020-03-04 cut (outside file coverage)")116ax.set_xlabel("date"); ax.set_ylabel("5-year GoC benchmark yield (%)")117ax.set_title(f"BD.CDN.5YR.DQ.YLD, {sample_start} to {sample_end} (post_cut=1 rows: {n_treated})")118ax.legend(loc="upper right"); fig.tight_layout()119fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)120plt.close(fig)121122res = {123 "task_id": "F6-conf-02",124 "sample_start": sample_start, "sample_end": sample_end,125 "n_treated": n_treated,126 "baseline_metric": baseline_metric,127 "model_metric": rmse,128 "effect_estimate": effect, "effect_ci_low": lo, "effect_ci_high": hi,129 "conclusion": conclusion,130 "feasibility": feas, "feasibility_reason": reason,131 "data_issues": issues,132 "n_input_rows": {FN: n_raw},133 "n_analysed": n_analysed,134 "random_seed": None,135 # auxiliary audit values136 "n_changes": n_chg,137 "n_missing_rows": n_missing, "n_duplicate_dates": n_dup,138 "baseline_change_sd_bp": float(df["chg_bp"].std(ddof=1)),139}140with open(os.path.join(OUT, "results.json"), "w") as f:141 json.dump(res, f, indent=2)142print(json.dumps({k: res[k] for k in ["sample_start", "sample_end", "n_treated",143 "baseline_metric", "effect_estimate", "conclusion", "feasibility"]}, indent=2))144