"""F6-conf-01: weekday effects in daily changes of the 10-year Treasury yield (DGS10).
Run from workspace root: ~/wp29/venv/bin/python out/analysis.py
Reads only data/DGS10.csv; writes out/results.json, out/table1.csv, out/figure1.png, 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 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")
assert raw["date"].is_monotonic_increasing and not raw["date"].duplicated().any()
assert (raw["date"].dt.dayofweek < 5).all()
empty = raw["DGS10"].str.strip() == ""
n_missing = int(empty.sum())
obs = raw.loc[~empty].copy()
obs["y"] = obs["DGS10"].astype(float)
obs = obs.reset_index(drop=True)

# One-day change between consecutive available observations, dated t
obs["dy"] = 100.0 * obs["y"].diff()
obs["gap_days"] = obs["date"].diff().dt.days
ch = obs.iloc[1:].copy()  # first available observation has no change
ch["dy"] = ch["dy"].round(10)  # remove float noise from 2-decimal percent data
ch["wd"] = ch["date"].dt.dayofweek
n_changes = len(ch)

excl = pd.DataFrame([
    ["1", "drop rows with empty DGS10 (market holidays)", n_raw, n_raw - n_missing, n_missing,
     "TASK rule: no yield published; dropped before differencing"],
    ["2", "first available observation has no change", n_raw - n_missing, n_changes, 1,
     "Δy undefined for first observation"],
], columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"])
excl.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) == n_changes

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

# Table 1
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, float_format="%.6f")

# (ii) Monday contrast with Welch 95% CI
mon, oth = groups[0], np.concatenate(groups[1:])
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
welch_p = 2 * stats.t.sf(abs(est / se), df_w)

# (iii) out-of-sample RMSE
wmeans = train.groupby("wd")["dy"].mean()
pred = 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() - pred) ** 2)))
# Diebold-Mariano-type check on squared-error loss differential (HAC not needed for h=1; iid SE)
d = (test["dy"].to_numpy() - pred) ** 2 - test["dy"].to_numpy() ** 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:
    conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")
else:
    conclusion = "no_detectable_effect"

gap_long = ch[ch["gap_days"] > 3]
res = {
    "task_id": "F6-conf-01",
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": [
        {"issue": f"{n_missing} rows with empty DGS10 field (US market holidays)",
         "handling": "dropped before differencing, per TASK rule"},
        {"issue": f"{len(gap_long)} changes span more than 3 calendar days (holiday gaps, e.g. after Monday holidays or 2001-09-11 closure)",
         "handling": "kept as one-day changes between consecutive available observations, attributed to weekday of t, per TASK rule"},
        {"issue": "heavy-tailed changes (max |dy| = %.0f bp on %s)" % (ch.dy.abs().max(), ch.loc[ch.dy.abs().idxmax(), "observation_date"]),
         "handling": "kept: large moves are plausible market events (the largest coincides with the 2009-03-18 FOMC announcement); no outlier removal; Welch CI used for the contrast"},
    ],
    "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": rmse0,
    "model_metric": rmse1,
    "conclusion": conclusion,
    # supplementary
    "effect_se": float(se),
    "effect_welch_df": float(df_w),
    "effect_welch_p_value": float(welch_p),
    "oos_mse_diff_mean": float(d.mean()),
    "oos_dm_t_stat": float(dm_t),
    "oos_dm_p_value": float(dm_p),
    "train_weekday_means_bp": {DAYS[k]: float(wmeans[k]) for k in range(5)},
    "rmse_ratio_model_to_baseline": rmse1 / rmse0,
}
with open(f"{OUT}/results.json", "w") as f:
    json.dump(res, f, indent=2)

# Figure 1
ci = stats.t.ppf(0.975, tab["n"] - 1) * tab["se_bp"]
fig, ax = plt.subplots(figsize=(6, 4))
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 Treasury yield change by weekday, 2000-2019 (n={n_changes})\nANOVA F={F:.2f}, p={p:.3f}; bars = 95% CI")
fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120)
print(json.dumps({k: res[k] for k in ["n_analysed", "anova_f_stat", "anova_p_value", "effect_estimate",
      "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))