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

DATA = "data/DGS10.csv"; OUT = "out"
os.makedirs(OUT, exist_ok=True)
DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]

sha = hashlib.sha256(open(DATA, "rb").read()).hexdigest()
raw = pd.read_csv(DATA, dtype={"observation_date": str, "DGS10": 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)
n_avail = len(d)
d["dy"] = 100.0 * d["y"].diff()
d["dy"] = d["dy"].round(6)  # remove float noise from 2-decimal percent inputs
ch = d.iloc[1:].copy()  # first available observation has no change
ch["wd"] = ch["date"].dt.weekday
assert ch["wd"].between(0, 4).all()
ch["weekday"] = ch["wd"].map(dict(enumerate(DAYS)))

pd.DataFrame([
    ["1", "drop rows with empty DGS10 (US market holidays)", n_raw, n_avail, n_raw - n_avail,
     "no yield published; required by TASK information rules"],
    ["2", "first available observation has no change", n_avail, len(ch), n_avail - len(ch),
     "change undefined for first observation"],
], columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).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)

# (ii) Monday contrast, Welch two-sample t interval
mon = groups[0]; 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
# robustness (reported in stdout only, not a required key): pooled t and HC1 regression
import statsmodels.api as sm
X = sm.add_constant((ch["wd"] == 0).astype(float).values)
hc = sm.OLS(ch["dy"].values, X).fit(cov_type="HC1")
print(f"[robustness] HC1 regression Monday coef {hc.params[1]:.4f}, CI {hc.conf_int()[1]}")

# (iii) out-of-sample
wd_means_train = train.groupby("wd")["dy"].mean()
fc = test["wd"].map(wd_means_train)
assert fc.notna().all()
rmse_base = float(np.sqrt(np.mean(test["dy"].values ** 2)))
rmse_model = float(np.sqrt(np.mean((test["dy"].values - fc.values) ** 2)))
rmse_trainmean = float(np.sqrt(np.mean((test["dy"].values - train["dy"].mean()) ** 2)))
print(f"[context] RMSE of constant estimation-period mean forecast: {rmse_trainmean:.4f}")
print(f"[context] estimation-period weekday means (bp): {wd_means_train.round(4).to_dict()}")

# table + figure
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(f"{OUT}/table1.csv", index=False)
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(DAYS)
ax.set_ylabel("Mean one-day change in DGS10 (bp)")
ax.set_title(f"10y UST daily change by weekday, 2000-2019 (n={len(ch)})\n95% CIs; ANOVA F={F:.2f}, p={p:.3f}")
fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)

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

gaps = d["date"].diff().dt.days.iloc[1:].values
normal = np.where(ch["wd"].values == 0, 3, 1)
n_span = int((gaps > normal).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(est), "effect_ci_low": float(lo), "effect_ci_high": float(hi),
    "baseline_metric": rmse_base, "model_metric": rmse_model,
    "conclusion": conclusion,
    "feasibility": "complete", "feasibility_reason": None,
    "data_issues": [
        {"issue": f"{n_missing} rows with empty DGS10 field (US market holidays, incl. 2001-09-11/12)",
         "handling": "dropped before differencing, per TASK rules"},
        {"issue": f"{n_span} changes span a holiday (gap longer than the normal 1 day, or 3 days for Mondays); Monday changes always span a weekend",
         "handling": "kept and attributed to the weekday of t, as specified; no rescaling by elapsed days"},
        {"issue": "large outliers, e.g. -51 bp on 2009-03-18 (FOMC QE1 announcement)",
         "handling": "kept; verified as genuine market moves, not data errors"},
    ],
    "n_input_rows": {"DGS10.csv": int(n_raw)},
    "n_analysed": int(len(ch)),
    "random_seed": None,
    "monday_contrast_method": "Welch two-sample t interval (unequal variances), Monday vs Tue-Fri, all changes",
    "monday_contrast_se": float(se), "monday_contrast_welch_df": float(df_w),
    "robust_hc1_ci_low": float(hc.conf_int()[1][0]), "robust_hc1_ci_high": float(hc.conf_int()[1][1]),
    "rmse_ratio_model_over_baseline": rmse_model / rmse_base,
    "rmse_constant_train_mean": rmse_trainmean,
    "data_sha256": sha,
}
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))