"""Day-of-week effect in daily 10y Treasury yield changes (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
OUT = "out"
raw = pd.read_csv("data/DGS10.csv", dtype={"DGS10": str}, keep_default_na=False)
n_input = len(raw)
issues = []
raw["observation_date"] = pd.to_datetime(raw["observation_date"])
val = raw["DGS10"].str.strip()
missing = (val == "") | (val == ".")
n_missing = int(missing.sum())
issues.append({"issue": f"{n_missing} rows with empty DGS10 field (US market holidays / no publication)",
               "handling": "Dropped before computing changes, per TASK rules; the change after a holiday spans the gap and is attributed to the weekday of its end date."})
df = raw.loc[~missing].copy()
df["y"] = pd.to_numeric(df["DGS10"], errors="coerce")
bad = int(df["y"].isna().sum())
if bad:
    issues.append({"issue": f"{bad} non-numeric DGS10 values", "handling": "Dropped."})
df = df.dropna(subset=["y"])
ndup = int(df["observation_date"].duplicated().sum())
nwkend = int((df["observation_date"].dt.dayofweek >= 5).sum())
if ndup or nwkend:
    issues.append({"issue": f"{ndup} duplicate dates, {nwkend} weekend dates", "handling": "Reported only; none expected."})
df = df.sort_values("observation_date").reset_index(drop=True)
df["dy"] = 100 * df["y"].diff()
gap = df["observation_date"].diff().dt.days
n_multi = int(((gap > 1) & (df["observation_date"].dt.dayofweek != 0)).sum() + ((gap > 3) & (df["observation_date"].dt.dayofweek == 0)).sum())
issues.append({"issue": f"{n_multi} changes span a holiday (more calendar/business days than usual), e.g. Tuesday after a Monday holiday",
               "handling": "Kept as defined by the TASK (consecutive available observations); noted as a limitation."})
issues.append({"issue": "Yields are quoted to 2 decimals, so changes are discrete in 1 bp steps",
               "handling": "No action; immaterial for means/ANOVA."})
ch = df.dropna(subset=["dy"]).copy()
ch["wd"] = ch["observation_date"].dt.dayofweek
names = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]

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

# Monday contrast: Welch two-sample
mon, oth = groups[0], ch.loc[ch.wd != 0, "dy"].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)
dfw = (v1 + v2) ** 2 / (v1**2 / (len(mon) - 1) + v2**2 / (len(oth) - 1))
t = stats.t.ppf(0.975, dfw)
lo, hi = eff - t * se, eff + t * se
welch_p = 2 * stats.t.sf(abs(eff / se), dfw)

# Out-of-sample
train = ch[ch.observation_date <= "2009-12-31"]
test = ch[ch.observation_date >= "2010-01-01"]
wmeans = train.groupby("wd")["dy"].mean()
pred = test["wd"].map(wmeans).values
rmse_base = float(np.sqrt(np.mean(test["dy"].values ** 2)))
rmse_model = float(np.sqrt(np.mean((test["dy"].values - pred) ** 2)))
# Diebold-Mariano (squared loss, simple HAC-free since 1-step) for descriptive purposes
d = test["dy"].values ** 2 - (test["dy"].values - pred) ** 2
dm_t = d.mean() / (d.std(ddof=1) / np.sqrt(len(d)))
dm_p = 2 * stats.norm.sf(abs(dm_t))

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
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(os.path.join(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(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(os.path.join(OUT, "figure1.png"), dpi=150)

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": rmse_base, "model_metric": rmse_model,
    "conclusion": concl, "feasibility": "complete", "feasibility_reason": None,
    "data_issues": issues, "n_input_rows": {"DGS10.csv": int(n_input)},
    "n_analysed": int(len(ch)), "random_seed": SEED,
    "supplementary": {"monday_contrast_method": "Welch two-sample t interval",
                      "welch_se_bp": float(se), "welch_df": float(dfw), "welch_p_value": float(welch_p),
                      "dm_stat_squared_loss": float(dm_t), "dm_p_value": float(dm_p),
                      "train_weekday_means_bp": {names[k]: float(v) for k, v in wmeans.items()},
                      "sd_all_changes_bp": float(ch.dy.std(ddof=1))},
}
with open(os.path.join(OUT, "results.json"), "w") as f:
    json.dump(res, f, indent=2)
print(json.dumps(res, indent=2)); print(tab)