"""Replication entry point for F6-conf-01: weekday effects in daily DGS10 changes."""
import json, os
import numpy as np, pandas as pd
from scipy import stats
import statsmodels.api as sm
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
OUT = os.path.join(ROOT, "out")
raw = pd.read_csv(os.path.join(ROOT, "data", "DGS10.csv"), dtype=str, keep_default_na=False)
n_input = len(raw)
issues = []
raw["date"] = pd.to_datetime(raw["observation_date"])
miss = raw["DGS10"].str.strip() == ""
n_missing = int(miss.sum())
issues.append({"issue": f"{n_missing} rows with empty DGS10 field (US market holidays)",
               "handling": "dropped before differencing; change computed between consecutive available observations (so post-holiday changes span >1 business day)"})
df = raw.loc[~miss].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"])
dup = int(df["date"].duplicated().sum())
if dup:
    issues.append({"issue": f"{dup} duplicated dates", "handling": "kept first"})
    df = df.drop_duplicates("date")
if not df["date"].is_monotonic_increasing:
    issues.append({"issue": "dates not sorted", "handling": "sorted"})
df = df.sort_values("date").reset_index(drop=True)
wkend = int((df["date"].dt.dayofweek >= 5).sum())
assert wkend == 0, "weekend observations present"

df["dy"] = 100 * df["y"].diff()
df["dy"] = df["dy"].round(6)  # remove float noise from decimal differencing
ch = df.dropna(subset=["dy"]).copy()
ch["wd"] = ch["date"].dt.dayofweek
names = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]
gap = (ch["date"] - df["date"].shift(1).loc[ch.index]).dt.days
n_long = int(((gap > 1) & ~((ch["wd"] == 0) & (gap == 3))).sum())
issues.append({"issue": f"{n_long} changes span a holiday gap (longer than the normal 1-day or Monday 3-day gap)",
               "handling": "kept, per task definition (attributed to weekday of t); noted as limitation"})

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

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

# Monday contrast: OLS on Monday indicator with HC3 robust SE
X = sm.add_constant((ch["wd"] == 0).astype(float).values)
fit = sm.OLS(ch["dy"].values, X).fit(cov_type="HC3")
eff = float(fit.params[1]); lo, hi = map(float, fit.conf_int(0.05)[1])
# Welch cross-check
w = stats.ttest_ind(groups[0], np.concatenate(groups[1:]), equal_var=False)

tm = train.groupby("wd")["dy"].mean()
pred = test["wd"].map(tm)
rmse_base = float(np.sqrt(np.mean(test["dy"] ** 2)))
rmse_model = float(np.sqrt(np.mean((test["dy"] - pred) ** 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(os.path.join(OUT, "table1.csv"), index=False, float_format="%.4f")

crit = stats.t.ppf(0.975, tab["n"] - 1)
fig, ax = plt.subplots(figsize=(6, 4))
ax.errorbar(range(5), tab["mean_bp"], yerr=crit * tab["se_bp"], 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(f"DGS10 daily change by weekday, 2000-2019 (95% CI)\nANOVA F={F:.2f}, p={p:.3f}")
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": eff, "effect_ci_low": lo, "effect_ci_high": hi,
       "baseline_metric": rmse_base, "model_metric": rmse_model,
       "conclusion": concl, "feasibility": "complete", "feasibility_reason": None,
       "data_issues": issues, "n_input_rows": {"DGS10.csv": n_input},
       "n_analysed": int(len(ch)), "random_seed": None,
       "supplementary": {"welch_t": float(w.statistic), "welch_p": float(w.pvalue),
                         "monday_contrast_hc3_p": float(fit.pvalues[1]),
                         "train_weekday_means_bp": {names[k]: float(v) for k, v in tm.items()},
                         "rmse_ratio_model_over_baseline": rmse_model / rmse_base,
                         "test_mean_change_bp": float(test["dy"].mean())}}
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)