"""Replication entry point for F6-conf-01. Run from workspace root:
   ~/wp29/venv/bin/python out/analysis.py"""
import json, os
import numpy as np, pandas as pd
from scipy import stats
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

SEED = None  # no randomness used
os.makedirs("out", exist_ok=True)
raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)
n_input = len(raw)
assert not raw.observation_date.duplicated().any()
n_missing = int((raw.DGS10.str.strip() == "").sum())
d = raw[raw.DGS10.str.strip() != ""].copy()
d["date"] = pd.to_datetime(d.observation_date)
d["y"] = d.DGS10.astype(float)
d = d.sort_values("date").reset_index(drop=True)
d["dy"] = 100 * d.y.diff()
d["gap_days"] = d.date.diff().dt.days
d = d.dropna(subset=["dy"]).copy()
d["dy"] = d.dy.round(10)
d["wd"] = d.date.dt.dayofweek
names = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]
assert d.wd.isin(range(5)).all()

train = d[d.date <= "2009-12-31"]; test = d[d.date >= "2010-01-01"]

# (i) one-way ANOVA, all changes
groups = [d.dy[d.wd == k].values for k in range(5)]
F, p = stats.f_oneway(*groups)

# (ii) Monday contrast, Welch two-sample interval
mon = d.dy[d.wd == 0].values; oth = d.dy[d.wd != 0].values
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 = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))
tq = stats.t.ppf(0.975, df)
lo, hi = eff - tq * se, eff + tq * se

# (iii) out-of-sample
wd_means = train.groupby("wd").dy.mean()
pred = test.wd.map(wd_means).values
rmse0 = float(np.sqrt(np.mean(test.dy.values ** 2)))
rmse1 = float(np.sqrt(np.mean((test.dy.values - pred) ** 2)))
# Diebold-Mariano (squared loss, h=1) as supplementary
dl = (test.dy.values - pred) ** 2 - test.dy.values ** 2
dm = dl.mean() / (dl.std(ddof=1) / np.sqrt(len(dl)))
dm_p = 2 * stats.norm.sf(abs(dm))
# supplementary: train-only mean (unconditional) forecast
rmse_mean = float(np.sqrt(np.mean((test.dy.values - train.dy.mean()) ** 2)))

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

# table + figure
tab = pd.DataFrame({"weekday": names,
                    "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("out/table1.csv", index=False, float_format="%.6f")
fig, ax = plt.subplots(figsize=(6, 4))
ci = [stats.t.ppf(0.975, n - 1) * s for n, s in zip(tab.n, tab.se_bp)]
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(names)
ax.set_ylabel("Mean one-day change in 10y yield (bp)")
ax.set_title(f"DGS10 daily changes by weekday, 2000-2019 (n={len(d)})\n95% CIs; ANOVA F={F:.2f}, p={p:.3f}")
fig.tight_layout(); fig.savefig("out/figure1.png", dpi=150)

post_hol = int(((d.gap_days > 1) & ~((d.wd == 0) & (d.gap_days == 3))).sum())
res = {
    "task_id": "F6-conf-01",
    "n_missing_rows": n_missing,
    "n_train": int(len(train)), "n_test": int(len(test)),
    "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,
    "feasibility": "complete", "feasibility_reason": None,
    "data_issues": [
        {"issue": f"{n_missing} rows with empty DGS10 field (US market holidays / no publication)",
         "handling": "Dropped before differencing, as prescribed; changes computed between consecutive available observations."},
        {"issue": f"{post_hol} changes (excluding ordinary Fri->Mon weekends) span a holiday gap and thus cover more than one business day",
         "handling": "Kept and attributed to the weekday of t, as prescribed by TASK.md; noted as a limitation (inflates variance, mostly for Tuesdays after Monday holidays)."},
    ],
    "n_input_rows": {"DGS10.csv": n_input},
    "n_analysed": int(len(d)),
    "random_seed": SEED,
    "supplementary": {
        "monday_contrast_method": "Welch two-sample t interval (unequal variances)",
        "welch_se_bp": float(se), "welch_df": float(df),
        "welch_p_value": float(2 * stats.t.sf(abs(eff / se), df)),
        "train_weekday_means_bp": {names[k]: float(v) for k, v in wd_means.items()},
        "rmse_train_unconditional_mean_bp": rmse_mean,
        "dm_stat_model_minus_zero_sqloss": float(dm), "dm_p_value": float(dm_p),
        "levene_p_value": float(stats.levene(*groups).pvalue),
        "kruskal_p_value": float(stats.kruskal(*groups).pvalue),
    },
}
with open("out/results.json", "w") as f:
    json.dump(res, f, indent=2)
print(json.dumps(res, indent=2)); print(tab)