"""Replication entry point for F6-conf-01. Run from workspace root."""
import json, os
import numpy as np, 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)
raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)
n_raw = len(raw)
missing = raw["DGS10"].str.strip() == ""
n_missing = int(missing.sum())
d = raw[~missing].copy()
d["date"] = pd.to_datetime(d["observation_date"])
d["y"] = d["DGS10"].astype(float)
d = d.sort_values("date").reset_index(drop=True)
gapdays = d["date"].diff().dt.days
d["dy"] = 100 * d["y"].diff()
ch = d.iloc[1:].copy()
ch["dy"] = ch["dy"].round(10)  # remove float noise from 2-decimal inputs
ch["wd"] = ch["date"].dt.dayofweek
names = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]
n_changes = len(ch)

# exclusion ledger
ex = pd.DataFrame([
    ["1", "drop rows with empty DGS10", n_raw, n_raw - n_missing, n_missing, "no yield published (market holiday)"],
    ["2", "first available observation has no change", n_raw - n_missing, n_changes, 1, "change undefined"],
], columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"])
ex.to_csv(f"{OUT}/exclusions.csv", index=False)

train = ch[ch["date"] <= "2009-12-31"]
test = ch[ch["date"] >= "2010-01-01"]

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

# Monday contrast, Welch
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 = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))
tq = stats.t.ppf(0.975, df)
lo, hi = est - tq * se, est + tq * se

# Out-of-sample
wd_means = train.groupby("wd")["dy"].mean()
pred = test["wd"].map(wd_means).values
rmse_base = float(np.sqrt(np.mean(test["dy"].values ** 2)))
rmse_model = float(np.sqrt(np.mean((test["dy"].values - pred) ** 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"

# Table and 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(f"{OUT}/table1.csv", index=False)
fig, ax = plt.subplots(figsize=(6, 4))
crit = [stats.t.ppf(0.975, n - 1) for n in tab["n"]]
ax.errorbar(range(5), tab["mean_bp"], yerr=np.array(crit) * tab["se_bp"], fmt="o", capsize=4)
ax.axhline(0, color="grey", lw=0.8)
ax.set_xticks(range(5)); ax.set_xticklabels([n[:3] for n in names])
ax.set_ylabel("Mean daily change in 10y yield (bp)")
ax.set_title("DGS10 daily change by weekday, 2000-2019 (95% CI)")
fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)

n_multiday = int(((gapdays.iloc[1:] == 2) | ((gapdays.iloc[1:] >= 3) & (d["date"].iloc[1:].dt.dayofweek != 0))
                  | (gapdays.iloc[1:] >= 4)).sum())
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, as required by TASK.md"},
        {"issue": f"{n_multiday} changes span more than one trading day because the previous weekday was a holiday",
         "handling": "kept and attributed to weekday of t, per TASK.md rules; noted as limitation"},
    ],
    "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(est),
    "effect_ci_low": float(lo),
    "effect_ci_high": float(hi),
    "baseline_metric": rmse_base,
    "model_metric": rmse_model,
    "conclusion": concl,
    "effect_se_welch": float(se),
    "effect_df_welch": float(df),
    "train_weekday_means_bp": {names[k]: float(v) for k, v in wd_means.items()},
    "rmse_ratio_model_over_baseline": rmse_model / rmse_base,
    "test_mean_bp": float(test["dy"].mean()),
}
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))