"""Replication entry point for F6-conf-02. Run from workspace root:
   ~/wp29/venv/bin/python out/analysis.py
Reads only data/gc5y_daily.csv; writes out/results.json, out/figure1.png,
out/table1.csv, out/exclusions.csv. No randomness is used."""
import json, os, hashlib
import numpy as np, pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy import stats

DATA = "data/gc5y_daily.csv"; OUT = "out"
os.makedirs(OUT, exist_ok=True)
raw = pd.read_csv(DATA)
n_raw = len(raw)
sha = hashlib.sha256(open(DATA, "rb").read()).hexdigest()
d = raw.copy(); d["date"] = pd.to_datetime(d["date"])

# ---- audit ----
issues = []
n_treated = int((d["post_cut"] == 1).sum())
n_na = int(d["yield_5y"].isna().sum()); n_dup = int(d["date"].duplicated().sum())
sorted_ok = bool(d["date"].is_monotonic_increasing)
bad_flag = int(((d["date"] >= "2020-03-04").astype(int) != d["post_cut"]).sum())
if n_treated == 0:
    issues.append({"issue": f"No rows with post_cut = 1: the file covers {raw['date'].iloc[0]} to {raw['date'].iloc[-1]} "
                   "(README retrieval used end_date=2019-12-31), entirely before the 2020-03-04 cut, so the post-cut window is empty.",
                   "handling": "Estimand declared not computable (premise_invalid); effect, CI, model_metric and n_analysed set to null. No proxy or substitute data used."})
if n_na: issues.append({"issue": f"{n_na} missing yields", "handling": "Changes involving missing values dropped."})
if n_dup: issues.append({"issue": f"{n_dup} duplicate dates", "handling": "Reported; not removed."})
if not sorted_ok: issues.append({"issue": "Dates not sorted", "handling": "Sorted by date before differencing."})
if 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."})
gaps = d["date"].diff().dt.days
issues.append({"issue": f"Calendar gaps between consecutive rows up to {int(gaps.max())} days (weekends/holidays; max gaps around year-end).",
               "handling": "Per task rules, one-day change is defined between consecutive rows regardless of calendar gap; no adjustment."})
d = d.sort_values("date").reset_index(drop=True)

# ---- changes and baseline ----
d["chg_bp"] = 100 * d["yield_5y"].diff()
excl = [("1", "drop first row (no prior row for one-day change)", n_raw, n_raw - 1, 1,
         "Change undefined for first observation")]
chg = d.dropna(subset=["chg_bp"])
baseline = float(chg["chg_bp"].mean())

# ---- main estimate ----
effect = lo = hi = rmse = None; n_an = None
if n_treated > 0:
    first = d.index[d["post_cut"] == 1][0]
    post = d.loc[d["post_cut"] == 1, "chg_bp"].dropna()
    pre = d.loc[max(first - 120, 1):first - 1, "chg_bp"].dropna()
    if len(pre) >= 2 and len(post) >= 2:
        effect = float(post.mean() - pre.mean())
        v1, v0 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)
        se = np.sqrt(v1 + v0)
        df = (v1 + v0) ** 2 / (v1 ** 2 / (len(post) - 1) + v0 ** 2 / (len(pre) - 1))
        t = stats.t.ppf(0.975, df); lo, hi = effect - t * se, effect + t * se
        resid = np.concatenate([post - post.mean(), pre - pre.mean()])
        rmse = float(np.sqrt(np.mean(resid ** 2))); n_an = int(len(post) + len(pre))
if effect is None:
    concl, feas = "premise_invalid", "infeasible"
    reason = ("The file contains no rows with post_cut = 1 (sample ends 2019-12-31, before the 2020-03-04 cut), "
              "so neither the post-cut window nor the '120 trading days before the first flagged day' window exists; "
              "the effect, its interval and the two-window RMSE cannot be computed.")
else:
    feas, reason = "complete", None
    concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")
excl.append(("2", "main-estimate windows (post_cut=1 rows + 120 preceding rows)", n_raw - 1,
             0 if n_an is None else n_an, (n_raw - 1) - (0 if n_an is None else n_an),
             "No post_cut=1 rows exist; window empty" if n_an is None else "Restrict to estimand windows"))
pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(f"{OUT}/exclusions.csv", index=False)

# ---- table and figure ----
d["year"] = d["date"].dt.year
tab = d.groupby("year").agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),
                            mean_change_bp=("chg_bp", "mean")).reset_index().round(4)
tab.to_csv(f"{OUT}/table1.csv", index=False)
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(d["date"], d["yield_5y"], lw=0.8)
ax.set_title("Government of Canada 5-year benchmark yield (BD.CDN.5YR.DQ.YLD)\n"
             "No post-cut (2020-03-04+) observations in file")
ax.set_ylabel("Yield (%)"); ax.set_xlabel("Date"); fig.tight_layout()
fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)

r = lambda x: None if x is None else round(float(x), 6)
res = {"task_id": "F6-conf-02", "sample_start": raw["date"].iloc[0], "sample_end": raw["date"].iloc[-1],
       "n_treated": n_treated, "baseline_metric": r(baseline), "model_metric": r(rmse),
       "effect_estimate": r(effect), "effect_ci_low": r(lo), "effect_ci_high": r(hi),
       "conclusion": concl, "feasibility": feas, "feasibility_reason": reason,
       "data_issues": issues, "n_input_rows": {"gc5y_daily.csv": int(n_raw)},
       "n_analysed": n_an, "random_seed": None,
       "audit_n_changes": int(len(chg)), "audit_sd_change_bp": r(chg["chg_bp"].std(ddof=1)),
       "audit_input_sha256": sha}
json.dump(res, open(f"{OUT}/results.json", "w"), indent=2)
print(f"rows={n_raw} treated={n_treated} baseline={baseline:.4f} bp conclusion={concl}")