"""Replication script for F6-conf-01: day-of-week effects in daily 10y Treasury yield changes."""
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_input = len(raw)
issues = []
raw["date"] = pd.to_datetime(raw["observation_date"], format="%Y-%m-%d")
empty = raw["DGS10"].str.strip() == ""
n_missing = int(empty.sum())
issues.append({"issue": f"{n_missing} rows with empty DGS10 field (US market holidays / no publication)",
               "handling": "rows dropped before computing changes, as required; the change after a gap spans >1 business day and is attributed to the weekday of its end date"})
vals = pd.to_numeric(raw.loc[~empty, "DGS10"], errors="coerce")
n_nonnum = int(vals.isna().sum())
if n_nonnum:
    issues.append({"issue": f"{n_nonnum} non-numeric DGS10 values", "handling": "treated as missing and dropped"})
df = raw.loc[~empty].assign(y=vals).dropna(subset=["y"]).sort_values("date")
n_dup = int(df["date"].duplicated().sum())
if n_dup:
    issues.append({"issue": f"{n_dup} duplicated dates", "handling": "kept first occurrence"})
    df = df.drop_duplicates("date")
n_weekend = int((df["date"].dt.dayofweek >= 5).sum())
if n_weekend:
    issues.append({"issue": f"{n_weekend} weekend observations", "handling": "dropped"})
    df = df[df["date"].dt.dayofweek < 5]
df["dy"] = 100 * df["y"].diff()
df = df.dropna(subset=["dy"]).copy()
gap = df["date"].diff().dt.days
df["wd"] = df["date"].dt.dayofweek
# count changes spanning a holiday (calendar gap larger than normal)
normal = np.where(df["wd"] == 0, 3, 1)
prev_dates = df["date"].shift(1)
span = (df["date"] - prev_dates).dt.days.values
long_span = int(np.nansum(span[1:] > normal[1:]))
issues.append({"issue": f"{long_span} changes span a holiday gap (longer than the usual 1 or 3 calendar days); this mixes horizons across weekdays (e.g. Tuesdays after Monday holidays)",
               "handling": "kept as specified by TASK.md; noted as a limitation"})
names = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]
df["weekday"] = df["wd"].map(dict(enumerate(names)))

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

groups = [df.loc[df["wd"] == k, "dy"].values for k in range(5)]
F, p = stats.f_oneway(*groups)

mon = df.loc[df["wd"] == 0, "dy"].values; oth = df.loc[df["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)
dof = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))
tq = stats.t.ppf(0.975, dof)
lo, hi = est - tq * se, est + tq * se
welch_p = float(2 * stats.t.sf(abs(est / se), dof))

wmeans = train.groupby("wd")["dy"].mean()
pred = test["wd"].map(wmeans).values
rmse0 = float(np.sqrt(np.mean(test["dy"].values ** 2)))
rmse1 = float(np.sqrt(np.mean((test["dy"].values - pred) ** 2)))
# supplementary: training-mean (no weekday) forecast
rmse_mean = float(np.sqrt(np.mean((test["dy"].values - train["dy"].mean()) ** 2)))

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"

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

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, ls="--")
ax.set_xticks(range(5)); ax.set_xticklabels(names)
ax.set_ylabel("Mean one-day change in DGS10 (bp)")
ax.set_title(f"10y Treasury yield change by weekday, 2000-2019\n(95% CI; ANOVA F={F:.2f}, p={p:.3f})")
fig.tight_layout(); fig.savefig(f"{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(est), "effect_ci_low": float(lo), "effect_ci_high": float(hi),
       "baseline_metric": rmse0, "model_metric": rmse1, "conclusion": concl,
       "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,
       "n_input_rows": {"DGS10.csv": int(n_input)}, "n_analysed": int(len(df)), "random_seed": None,
       "supplementary": {"monday_contrast_method": "Welch two-sample t (Monday vs other weekdays)",
                         "welch_se_bp": float(se), "welch_df": float(dof), "welch_p_value": welch_p,
                         "train_weekday_means_bp": {names[k]: float(v) for k, v in wmeans.items()},
                         "rmse_train_overall_mean_bp": rmse_mean,
                         "rmse_ratio_model_to_baseline": rmse1 / rmse0,
                         "n_valid_obs": int(len(df) + 1)}}
with open(f"{OUT}/results.json", "w") as f:
    json.dump(res, f, indent=2)
print(json.dumps(res, indent=2)); print(tab)