"""F6-conf-01: weekday effects in daily changes of the 10y Treasury yield (DGS10).
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.
No randomness is used (random_seed = null)."""
import json, os
import numpy as np, pandas as pd
from scipy import stats
import statsmodels.api as sm
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

OUT = "out"; os.makedirs(OUT, exist_ok=True)
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[~empty].copy()
d["date"] = pd.to_datetime(d["observation_date"], format="%Y-%m-%d")
d["y"] = pd.to_numeric(d["DGS10"], errors="raise")
assert d["date"].is_monotonic_increasing and not d["date"].duplicated().any()
d["dy"] = 100.0 * d["y"].diff()
d["dy"] = d["dy"].round(8)  # remove float noise; data have 1 bp resolution
ch = d.dropna(subset=["dy"]).copy()
ch["wd"] = ch["date"].dt.dayofweek
assert ch["wd"].between(0, 4).all()
ch["weekday"] = ch["wd"].map(dict(enumerate(DAYS)))

pd.DataFrame([
    dict(step=1, rule="drop rows with empty DGS10 field (holidays)", n_before=n_raw,
         n_after=n_raw - n_missing, n_removed=n_missing, rationale="TASK rule: no published yield"),
    dict(step=2, rule="first available observation has no change", n_before=n_raw - n_missing,
         n_after=len(ch), n_removed=n_raw - n_missing - len(ch), rationale="TASK rule: change needs y_{t-1}"),
]).to_csv(f"{OUT}/exclusions.csv", index=False)

train = ch[ch["date"] <= "2009-12-31"]; test = ch[ch["date"] >= "2010-01-01"]
assert len(train) + len(test) == len(ch)

# (i) classical one-way ANOVA, all changes
groups = [ch.loc[ch["wd"] == k, "dy"].values for k in range(5)]
F, p = stats.f_oneway(*groups)
kw_H, kw_p = stats.kruskal(*groups)                       # sensitivity (rank-based)
lev_W, lev_p = stats.levene(*groups, center="median")     # equal-variance check
from statsmodels.stats.oneway import anova_oneway
wa = anova_oneway(groups, use_var="unequal")              # Welch ANOVA sensitivity

# (ii) Monday contrast, Welch two-sample interval (primary)
mon = ch.loc[ch["wd"] == 0, "dy"].values; oth = ch.loc[ch["wd"] != 0, "dy"].values
est = 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))
tq = stats.t.ppf(0.975, df_w)
lo, hi = est - tq * se, est + tq * se
welch_p = 2 * stats.t.sf(abs(est / se), df_w)
# sensitivity: OLS on Monday dummy with HAC (Newey-West, 5 lags) SE
X = sm.add_constant((ch["wd"] == 0).astype(float).values)
hac = sm.OLS(ch["dy"].values, X).fit(cov_type="HAC", cov_kwds={"maxlags": 5})

# (iii) out-of-sample
wmeans = train.groupby("wd")["dy"].mean()
f_model = test["wd"].map(wmeans).values
e0 = test["dy"].values; e1 = e0 - f_model
rmse0 = float(np.sqrt(np.mean(e0 ** 2))); rmse1 = float(np.sqrt(np.mean(e1 ** 2)))
e_unc = e0 - train["dy"].mean()
rmse_unc = float(np.sqrt(np.mean(e_unc ** 2)))
# Diebold-Mariano on squared-error loss differential (model - baseline), HAC 5 lags
dl = e1 ** 2 - e0 ** 2
dm = sm.OLS(dl, np.ones_like(dl)).fit(cov_type="HAC", cov_kwds={"maxlags": 5})

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"

# table + figure
tab = ch.groupby("wd")["dy"].agg(n="count", mean_bp="mean", sd="std").reindex(range(5))
tab["se_bp"] = tab["sd"] / np.sqrt(tab["n"])
tab.insert(0, "weekday", DAYS)
tab[["weekday", "n", "mean_bp", "se_bp"]].to_csv(f"{OUT}/table1.csv", index=False, float_format="%.6f")
crit = stats.t.ppf(0.975, tab["n"] - 1)
fig, ax = plt.subplots(figsize=(6, 4))
ax.errorbar(range(5), tab["mean_bp"], yerr=crit * tab["se_bp"], fmt="o", capsize=5, color="C0")
ax.axhline(0, color="grey", lw=0.8, ls="--")
ax.set_xticks(range(5)); ax.set_xticklabels([s[:3] for s in DAYS])
ax.set_ylabel("Mean one-day change in DGS10 (bp)")
ax.set_title(f"10y Treasury yield changes by weekday, 2000-2019 (n={len(ch)})\n"
             f"95% CI; one-way ANOVA F={F:.2f}, p={p:.3f}", fontsize=9)
fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=150); plt.close(fig)

r = lambda x: float(round(float(x), 10))
res = {
    "task_id": "F6-conf-01", "feasibility": "complete", "feasibility_reason": None,
    "data_issues": [
        {"issue": f"{n_missing} rows with empty DGS10 field (US market holidays)",
         "handling": "dropped before differencing, per TASK rule"},
        {"issue": "changes following a holiday span more than one trading day (e.g. Fri->Tue after Monday holidays)",
         "handling": "kept and attributed to weekday of t, per TASK rule; noted as limitation"},
        {"issue": "extreme changes (min -51 bp on 2009-03-18, QE1 announcement; max +25 bp); heavy tails",
         "handling": "kept (genuine market events); rank-based Kruskal-Wallis reported as sensitivity"},
        {"issue": "yields recorded to 0.01 pct, so changes are discretised at 1 bp",
         "handling": "none needed; noted"},
    ],
    "n_input_rows": {"DGS10.csv": int(n_raw)},
    "n_analysed": int(len(ch)), "random_seed": None,
    "n_missing_rows": n_missing, "n_train": int(len(train)), "n_test": int(len(test)),
    "anova_f_stat": r(F), "anova_p_value": r(p),
    "effect_estimate": r(est), "effect_ci_low": r(lo), "effect_ci_high": r(hi),
    "baseline_metric": r(rmse0), "model_metric": r(rmse1),
    "conclusion": concl,
    # supplementary
    "effect_se_welch": r(se), "effect_df_welch": r(df_w), "effect_p_welch": r(welch_p),
    "effect_ci_method": "Welch two-sample t interval (Monday vs Tue-Fri), all changes",
    "sens_hac_effect_se": r(hac.bse[1]), "sens_hac_effect_p": r(hac.pvalues[1]),
    "sens_kruskal_H": r(kw_H), "sens_kruskal_p": r(kw_p),
    "sens_welch_anova_F": r(wa.statistic), "sens_welch_anova_p": r(wa.pvalue),
    "sens_levene_W": r(lev_W), "sens_levene_p": r(lev_p),
    "oos_rmse_ratio_model_over_baseline": r(rmse1 / rmse0),
    "oos_rmse_unconditional_train_mean": r(rmse_unc),
    "oos_dm_mean_loss_diff": r(dm.params[0]), "oos_dm_t": r(dm.tvalues[0]), "oos_dm_p": r(dm.pvalues[0]),
    "train_weekday_means_bp": {DAYS[k]: r(v) for k, v in wmeans.items()},
    "n_changes_total": int(len(ch)),
}
with open(f"{OUT}/results.json", "w") as f:
    json.dump(res, f, indent=2)
print(json.dumps({k: res[k] for k in ["n_analysed", "anova_f_stat", "anova_p_value", "effect_estimate",
      "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))