"""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/figure1.png, out/table1.csv, 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.formula.api as smf
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

OUT = "out"; os.makedirs(OUT, exist_ok=True)
DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]
SPLIT = pd.Timestamp("2009-12-31")

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"] = pd.to_numeric(d["DGS10"], errors="raise")
d = d.sort_values("date").reset_index(drop=True)
assert not d["date"].duplicated().any()
d["dy"] = 100.0 * d["y"].diff()
d["gap"] = d["date"].diff().dt.days
ch = d.iloc[1:].copy()  # first available observation has no change
ch["wd"] = ch["date"].dt.dayofweek
assert ch["wd"].between(0, 4).all()
n_changes = len(ch)

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

train = ch[ch["date"] <= SPLIT]; test = ch[ch["date"] > SPLIT]
assert train["date"].max() < test["date"].min()

# (i) classical one-way ANOVA, all changes
groups = [ch.loc[ch.wd == k, "dy"].to_numpy() for k in range(5)]
F, p = stats.f_oneway(*groups)
# robustness (not a required key): Welch-type check via Kruskal-Wallis and Levene, reported in stdout only
kw = stats.kruskal(*groups); lev = stats.levene(*groups)

# (ii) Monday contrast, Welch two-sample t interval
mon = ch.loc[ch.wd == 0, "dy"].to_numpy(); oth = ch.loc[ch.wd != 0, "dy"].to_numpy()
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

# robustness (extra keys): Monday-indicator regression with HC1 and Newey-West(5) SEs; ANOVA excluding 2008-09 to 2009-06 crisis window
ch["mon"] = (ch["wd"] == 0).astype(int)
hc1 = smf.ols("dy ~ mon", ch).fit(cov_type="HC1").conf_int().loc["mon"].tolist()
nw5 = smf.ols("dy ~ mon", ch).fit(cov_type="HAC", cov_kwds={"maxlags": 5}).conf_int().loc["mon"].tolist()
xc = ch[(ch["date"] < "2008-09-01") | (ch["date"] > "2009-06-30")]
p_xc = stats.f_oneway(*[xc.loc[xc.wd == k, "dy"] for k in range(5)]).pvalue

# (iii) out-of-sample RMSE
wmeans = train.groupby("wd")["dy"].mean()
fc = test["wd"].map(wmeans).to_numpy()
rmse0 = float(np.sqrt(np.mean(test["dy"].to_numpy() ** 2)))
rmse1 = float(np.sqrt(np.mean((test["dy"].to_numpy() - fc) ** 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 and figure
rows = []
for k, g in enumerate(groups):
    rows.append(dict(weekday=DAYS[k], n=len(g), mean_bp=g.mean(), se_bp=g.std(ddof=1) / np.sqrt(len(g))))
tab = pd.DataFrame(rows)
tab.to_csv(f"{OUT}/table1.csv", index=False)
fig, ax = plt.subplots(figsize=(6, 4))
ci = [stats.t.ppf(0.975, r["n"] - 1) * r["se_bp"] for r in rows]
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([s[:3] for s in DAYS])
ax.set_ylabel("Mean one-day change in DGS10 (bp)")
ax.set_title("10y UST yield: mean daily change by weekday, 2000-2019\n(95%% CI, all %d changes)" % n_changes)
fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)

gap_counts = ch["gap"].value_counts().sort_index()
res = {
    "task_id": "F6-conf-01",
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": [
        {"issue": f"{n_missing} rows with empty DGS10 (US market holidays incl. 2001-09-11/12)",
         "handling": "dropped before differencing, per TASK rules"},
        {"issue": f"{int((ch.gap > 3).sum() + ((ch.gap == 2)).sum())} changes span a holiday (calendar gap 2 or 4 days), i.e. cover more than one trading session",
         "handling": "kept and attributed to weekday of t as the TASK prescribes; noted as a limitation"},
        {"issue": "extreme changes (e.g. -51 bp on 2009-03-18, -27 bp on 2008-09-15) during the 2008-09 crisis",
         "handling": "kept (genuine market events); no outlier trimming; discussed as limitation"},
        {"issue": "yields quoted to 0.01 percentage point, so changes are discretised at 1 bp",
         "handling": "none needed; noted"},
    ],
    "n_input_rows": {"DGS10.csv": int(n_raw)},
    "n_analysed": int(n_changes),
    "random_seed": None,
    "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,
    "effect_se_welch": float(se),
    "effect_df_welch": float(df),
    "train_weekday_means_bp": {DAYS[k]: float(v) for k, v in wmeans.items()},
    "rmse_ratio_model_over_baseline": rmse1 / rmse0,
    "robust_kruskal_wallis_p": float(kw.pvalue),
    "robust_levene_p": float(lev.pvalue),
    "robust_effect_ci_hc1": [float(v) for v in hc1],
    "robust_effect_ci_nw5": [float(v) for v in nw5],
    "robust_anova_p_excl_crisis_2008_09_to_2009_06": float(p_xc),
    "n_excl_crisis": int(len(xc)),
}
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_train", "n_test", "anova_f_stat", "anova_p_value", "effect_estimate",
                                      "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))