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

OUT = "out"
os.makedirs(OUT, exist_ok=True)
SEED = None  # no randomness used
DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]

raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)
n_raw = len(raw)
empty = raw["DGS10"].str.strip() == ""
n_missing = int(empty.sum())
d = raw.loc[~empty].copy()
d["date"] = pd.to_datetime(d["observation_date"])
d["y"] = d["DGS10"].astype(float)
assert d["date"].is_monotonic_increasing and not d["date"].duplicated().any()
d = d.reset_index(drop=True)
d["dy"] = 100.0 * d["y"].diff()
d["gap_days"] = d["date"].diff().dt.days
ch = d.iloc[1:].copy()  # first available obs has no change
ch["wd"] = ch["date"].dt.dayofweek
assert ch["wd"].between(0, 4).all()
ch["weekday"] = ch["wd"].map(dict(enumerate(DAYS)))
n_ch = len(ch)

excl = pd.DataFrame([
    {"step": 1, "rule": "drop rows with empty DGS10 (market holidays)", "n_before": n_raw,
     "n_after": n_raw - n_missing, "n_removed": n_missing, "rationale": "no yield published; required by TASK info rules"},
    {"step": 2, "rule": "first available observation has no change", "n_before": n_raw - n_missing,
     "n_after": n_ch, "n_removed": n_raw - n_missing - n_ch, "rationale": "difference undefined"},
])
excl.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)

# (i) ANOVA, all changes
groups = [ch.loc[ch["wd"] == k, "dy"].values for k in range(5)]
F, p = stats.f_oneway(*groups)

# table 1
tab = pd.DataFrame({"weekday": DAYS,
                    "n": [len(g) for g in groups],
                    "mean_bp": [g.mean() for g in groups],
                    "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups]})
tab.to_csv(os.path.join(OUT, "table1.csv"), index=False)

# (ii) Monday contrast, Welch two-sample interval
mon, oth = groups[0], np.concatenate(groups[1:])
eff = mon.mean() - oth.mean()
v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)
se = np.sqrt(v1 + v2)
df_w = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))
tcrit = stats.t.ppf(0.975, df_w)
lo, hi = eff - tcrit * se, eff + tcrit * se
welch_p = 2 * stats.t.sf(abs(eff / se), df_w)

# (iii) out-of-sample
cut = pd.Timestamp("2009-12-31")
tr, te = ch[ch["date"] <= cut], ch[ch["date"] > cut]
wd_means = tr.groupby("wd")["dy"].mean()
assert tr["date"].max() < te["date"].min()
pred = te["wd"].map(wd_means).values
rmse0 = float(np.sqrt(np.mean(te["dy"].values ** 2)))
rmse1 = float(np.sqrt(np.mean((te["dy"].values - pred) ** 2)))
# supplementary: unconditional train-mean forecast
rmse_mean = float(np.sqrt(np.mean((te["dy"].values - tr["dy"].mean()) ** 2)))

if p >= 0.05:
    concl = "no_detectable_effect"
elif lo > 0:
    concl = "positive_effect"
elif hi < 0:
    concl = "negative_effect"
else:
    concl = "no_detectable_effect"

n_long_gap = int((ch["gap_days"] > 4).sum())
n_holiday_spanning = int(((ch["gap_days"] > 1) & ~((ch["wd"] == 0) & (ch["gap_days"] == 3))).sum())
issues = [
    {"issue": f"{n_missing} rows with empty DGS10 (US market holidays)",
     "handling": "dropped before differencing, as required; changes computed between consecutive available observations"},
    {"issue": f"{n_holiday_spanning} changes span a holiday (gap longer than a normal weekday/weekend gap); {n_long_gap} span more than 4 calendar days",
     "handling": "kept and attributed to the weekday of t per the task rules; noted as limitation (multi-day changes have larger variance)"},
]

res = {
    "task_id": "F6-conf-01",
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": issues,
    "n_input_rows": {"DGS10.csv": int(n_raw)},
    "n_analysed": int(n_ch),
    "random_seed": SEED,
    "n_missing_rows": n_missing,
    "n_train": int(len(tr)),
    "n_test": int(len(te)),
    "anova_f_stat": float(F),
    "anova_p_value": float(p),
    "effect_estimate": float(eff),
    "effect_ci_low": float(lo),
    "effect_ci_high": float(hi),
    "baseline_metric": rmse0,
    "model_metric": rmse1,
    "conclusion": concl,
    # supplementary diagnostics
    "supp_effect_se": float(se),
    "supp_effect_welch_df": float(df_w),
    "supp_effect_welch_p": float(welch_p),
    "supp_rmse_ratio_model_over_baseline": rmse1 / rmse0,
    "supp_rmse_train_mean_forecast": rmse_mean,
    "supp_levene_p": float(stats.levene(*groups, center="median").pvalue),
    "supp_kruskal_p": float(stats.kruskal(*groups).pvalue),
    "supp_train_weekday_means_bp": {DAYS[k]: float(v) for k, v in wd_means.items()},
    "supp_n_holiday_spanning_changes": n_holiday_spanning,
}
with open(os.path.join(OUT, "results.json"), "w") as f:
    json.dump(res, f, indent=2)

# figure 1
ci = stats.t.ppf(0.975, tab["n"] - 1) * tab["se_bp"]
fig, ax = plt.subplots(figsize=(6, 4))
ax.errorbar(range(5), tab["mean_bp"], yerr=ci, fmt="o", capsize=5)
ax.axhline(0, color="grey", lw=0.8)
ax.set_xticks(range(5)); ax.set_xticklabels(DAYS)
ax.set_ylabel("Mean one-day change in DGS10 (bp)")
ax.set_title("10y Treasury yield: mean daily change by weekday, 2000-2019\n(95% CI)")
fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)

print(json.dumps({k: res[k] for k in ["n_analysed", "n_train", "n_test", "anova_f_stat", "anova_p_value",
      "effect_estimate", "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))
print(tab.to_string(index=False))