"""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"
PRE_N = 120
os.makedirs(OUT, exist_ok=True)

raw = pd.read_csv(DATA, dtype={"date": str})
n_raw = len(raw)
sha = hashlib.sha256(open(DATA, "rb").read()).hexdigest()
issues = []

# ---------------- audit ----------------
df = raw.copy()
df["date_parsed"] = pd.to_datetime(df["date"], format="%Y-%m-%d", errors="coerce")
n_bad_date = int(df["date_parsed"].isna().sum())
df["yield_5y"] = pd.to_numeric(df["yield_5y"], errors="coerce")
n_na_y = int(df["yield_5y"].isna().sum())
n_na_t = int(df["post_cut"].isna().sum())
n_dup = int(df["date"].duplicated().sum())
monotone = bool(df["date_parsed"].is_monotonic_increasing)
if n_bad_date or n_na_y or n_na_t:
    issues.append({"issue": f"{n_bad_date} unparseable dates, {n_na_y} missing yields, {n_na_t} missing post_cut",
                   "handling": "reported; changes computed between consecutive rows as given"})
if n_dup:
    issues.append({"issue": f"{n_dup} duplicate dates", "handling": "reported; not removed"})
if not monotone:
    issues.append({"issue": "dates not sorted ascending", "handling": "file order kept as required by the change definition; reported"})

sample_start = str(raw["date"].iloc[0]); sample_end = str(raw["date"].iloc[-1])
treated = df["post_cut"] == 1
n_treated = int(treated.sum())
bad_flag = int(((df["date_parsed"] >= "2020-03-04") != treated).sum())
if bad_flag:
    issues.append({"issue": f"{bad_flag} rows where post_cut disagrees with the documented rule (date >= 2020-03-04)",
                   "handling": "post_cut used exactly as provided, per task rules"})

# no rows are removed for the descriptive baseline
excl = [{"step": "S0", "rule": "read file; no filtering applied (no missing/duplicate/unparseable rows found)",
         "n_before": n_raw, "n_after": n_raw, "n_removed": 0,
         "rationale": "audit found no rows requiring removal"}]

# ---------------- changes and baseline ----------------
df["d_bp"] = 100.0 * df["yield_5y"].diff()
baseline_metric = float(df["d_bp"].mean())  # mean over n_raw-1 one-day changes

# ---------------- estimand ----------------
model_metric = effect = lo = hi = None
n_analysed = None
if n_treated == 0:
    feasibility = "infeasible"; conclusion = "premise_invalid"
    feasibility_reason = (f"The file covers {sample_start} to {sample_end} and contains no rows with post_cut = 1 "
                          "(the cut was announced 2020-03-04, after the last date). The post-cut mean change, "
                          "and hence the effect, its interval and the two-window model RMSE, cannot be computed. "
                          "No proxy estimate is substituted.")
    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). "
                            "The README retrieval URL requests end_date=2019-12-31, consistent with this truncation.",
                   "handling": "Estimand declared not estimable: effect_estimate, CI, model_metric, n_analysed set to null; conclusion premise_invalid."})
    excl.append({"step": "S1", "rule": "keep rows in the two estimation windows (120 pre-window days + post_cut=1 days)",
                 "n_before": n_raw, "n_after": 0, "n_removed": n_raw,
                 "rationale": "no post_cut=1 rows exist, so the windows are undefined and no row enters the effect estimate"})
else:
    first = int(np.flatnonzero(treated.values)[0])
    post = df.loc[treated, "d_bp"].dropna()
    pre = df["d_bp"].iloc[max(first - PRE_N, 0):first].dropna()
    if len(pre) < PRE_N:
        issues.append({"issue": f"only {len(pre)} pre-window changes available (<{PRE_N})", "handling": "used what is available; reported"})
    if len(pre) >= 2 and len(post) >= 2:
        effect = float(post.mean() - pre.mean())
        v1, v2 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)
        se = np.sqrt(v1 + v2)
        dof = (v1 + v2) ** 2 / (v1 ** 2 / (len(post) - 1) + v2 ** 2 / (len(pre) - 1))
        q = stats.t.ppf(0.975, dof)
        lo, hi = float(effect - q * se), float(effect + q * se)
        resid = np.concatenate([post - post.mean(), pre - pre.mean()])
        model_metric = float(np.sqrt(np.mean(resid ** 2)))
        n_analysed = int(len(pre) + len(post))
        conclusion = "positive_effect" if lo > 0 else "negative_effect" if hi < 0 else "no_detectable_effect"
        feasibility, feasibility_reason = "complete", None
    else:
        feasibility, conclusion = "infeasible", "premise_invalid"
        feasibility_reason = "too few observations in one of the windows"
    excl.append({"step": "S1", "rule": "keep rows in the two estimation windows",
                 "n_before": n_raw, "n_after": n_analysed or 0, "n_removed": n_raw - (n_analysed or 0),
                 "rationale": "estimand is defined on these windows only"})

pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False)

# ---------------- table and figure ----------------
df["year"] = df["date_parsed"].dt.year
tab = df.groupby("year").agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),
                             mean_change_bp=("d_bp", "mean")).reset_index()
tab["mean_yield"] = tab["mean_yield"].round(6); tab["mean_change_bp"] = tab["mean_change_bp"].round(6)
tab.to_csv(f"{OUT}/table1.csv", index=False)

fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(df["date_parsed"], df["yield_5y"], lw=0.8)
ax.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)")
ax.set_ylabel("percent"); ax.set_xlabel("date")
fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)

res = {
    "task_id": "F6-conf-02",
    "sample_start": sample_start, "sample_end": sample_end,
    "n_treated": n_treated,
    "baseline_metric": baseline_metric,
    "model_metric": model_metric,
    "effect_estimate": effect, "effect_ci_low": lo, "effect_ci_high": hi,
    "conclusion": conclusion,
    "feasibility": feasibility, "feasibility_reason": feasibility_reason,
    "data_issues": issues,
    "n_input_rows": {"gc5y_daily.csv": n_raw},
    "n_analysed": n_analysed,
    "random_seed": None,
}
with open(f"{OUT}/results.json", "w") as fh:
    json.dump(res, fh, indent=2)
print(f"sha256={sha} rows={n_raw} {sample_start}..{sample_end} n_treated={n_treated}")
print(f"baseline_metric={baseline_metric:.6f} bp; conclusion={conclusion}; feasibility={feasibility}")