"""F6-conf-01: weekday effects in daily changes of the 10y Treasury yield (DGS10), 2000-2019.
Run from workspace root: ~/wp29/venv/bin/python out/analysis.py
Reads only data/DGS10.csv. No randomness is used (random_seed = null)."""
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)
DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]

raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)
n_raw = len(raw)
raw["date"] = pd.to_datetime(raw["observation_date"], format="%Y-%m-%d")
val = raw["DGS10"].str.strip()
missing = val == ""
n_missing = int(missing.sum())
obs = raw.loc[~missing, ["date"]].copy()
obs["y"] = val[~missing].astype(float).values
obs = obs.sort_values("date").reset_index(drop=True)
obs["dy"] = 100.0 * obs["y"].diff()
ch = obs.dropna(subset=["dy"]).copy()
ch["dy"] = ch["dy"].round(10)  # remove float noise from percent->bp conversion
ch["weekday"] = ch["date"].dt.dayofweek.map(dict(enumerate(DAYS)))
assert ch["weekday"].notna().all()

# exclusion ledger
excl = pd.DataFrame([
    dict(step=1, rule="drop rows with empty DGS10 field (no published yield)", n_before=n_raw,
         n_after=n_raw - n_missing, n_removed=n_missing, rationale="TASK rule: holidays/no publication"),
    dict(step=2, rule="first available observation has no prior value -> no change", n_before=n_raw - n_missing,
         n_after=len(ch), n_removed=n_raw - n_missing - len(ch), rationale="TASK rule: change undefined"),
])
excl.to_csv(f"{OUT}/exclusions.csv", index=False)

cut = pd.Timestamp("2009-12-31")
train = ch[ch["date"] <= cut]; test = ch[ch["date"] >= pd.Timestamp("2010-01-01")]
assert len(train) + len(test) == len(ch)

# (i) classical one-way ANOVA, all changes
groups = [ch.loc[ch.weekday == d, "dy"].values for d in DAYS]
F, p = stats.f_oneway(*groups)

# (ii) Monday contrast, Welch two-sample t interval
mon = ch.loc[ch.weekday == "Monday", "dy"].values
oth = ch.loc[ch.weekday != "Monday", "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))
tq = stats.t.ppf(0.975, dfw)
lo, hi = eff - tq * se, eff + tq * se
welch_p = stats.ttest_ind(mon, oth, equal_var=False).pvalue

# (iii) out-of-sample: baseline zero change vs train weekday means
wk_means_train = train.groupby("weekday")["dy"].mean()
pred = test["weekday"].map(wk_means_train)
assert pred.notna().all()
rmse_base = float(np.sqrt(np.mean(test["dy"].values ** 2)))
rmse_model = float(np.sqrt(np.mean((test["dy"].values - pred.values) ** 2)))
# supplementary: unconditional train mean forecast and DM-type test of model vs zero
rmse_trainmean = float(np.sqrt(np.mean((test["dy"].values - train["dy"].mean()) ** 2)))
dloss = (test["dy"].values - pred.values) ** 2 - test["dy"].values ** 2
dm_stat = dloss.mean() / (dloss.std(ddof=1) / np.sqrt(len(dloss)))
dm_p = 2 * stats.norm.sf(abs(dm_stat))

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

# table + figure
tab = ch.groupby("weekday")["dy"].agg(n="size", mean_bp="mean", sd="std").reindex(DAYS)
tab["se_bp"] = tab["sd"] / np.sqrt(tab["n"])
tab = tab.reset_index()[["weekday", "n", "mean_bp", "se_bp"]]
tab.to_csv(f"{OUT}/table1.csv", index=False, float_format="%.6f")
crit = stats.t.ppf(0.975, tab["n"] - 1)
fig, ax = plt.subplots(figsize=(6.5, 4))
ax.errorbar(range(5), tab["mean_bp"], yerr=crit * tab["se_bp"], fmt="o", capsize=5, color="C0")
ax.axhline(0, color="grey", lw=0.8, ls="--")
ax.set_xticks(range(5)); ax.set_xticklabels([d[:3] for d in DAYS])
ax.set_ylabel("Mean one-day change in DGS10 (bp)")
ax.set_title(f"10y Treasury yield change by weekday, 2000-2019 (n={len(ch)})\n"
             f"95% CIs; ANOVA F={F:.2f}, p={p:.3f}", fontsize=9)
fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=150); plt.close(fig)

issues = [
    {"issue": f"{n_missing} rows with empty DGS10 field (US market holidays / closures)",
     "handling": "dropped before computing changes, per TASK rule (see out/exclusions.csv)"},
    {"issue": "changes following a holiday span more than one business day (e.g. Tuesday after a Monday holiday)",
     "handling": "kept and attributed to the weekday of t, per TASK definition; noted as a limitation"},
    {"issue": "large outlier change of -51 bp on 2009-03-18 (FOMC QE announcement)",
     "handling": "genuine market move; kept (no outlier removal)"},
]
res = dict(
    task_id="F6-conf-01", feasibility="complete", feasibility_reason=None, data_issues=issues,
    n_input_rows={"DGS10.csv": int(n_raw)}, n_analysed=int(len(ch)), 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=rmse_base, model_metric=rmse_model, conclusion=conclusion,
    # supplementary keys
    effect_se=float(se), effect_welch_df=float(dfw), effect_welch_p_value=float(welch_p),
    effect_ci_method="Welch two-sample t interval (Monday vs Tue-Fri, all changes)",
    oos_rmse_ratio_model_to_baseline=rmse_model / rmse_base,
    oos_rmse_train_unconditional_mean=rmse_trainmean,
    oos_dm_stat_model_vs_zero=float(dm_stat), oos_dm_p_value=float(dm_p),
    train_weekday_means_bp={d: float(wk_means_train[d]) for d in DAYS},
    n_changes_total=int(len(ch)),
)
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))