"""F5-conf-02: finite-sample bias of OLS AR(1) coefficient (with intercept).
Reads only data/design.csv; writes out/results.json (and out/exclusions.csv). Deterministic."""
import json, os, hashlib
import numpy as np
import pandas as pd

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
DATA = os.path.join(ROOT, "data", "design.csv")
OUT = os.path.join(ROOT, "out")

raw = pd.read_csv(DATA)
n_rows = int(len(raw))
d = dict(zip(raw["parameter"].astype(str).str.strip(), raw["value"].astype(str).str.strip()))
data_issues = []
if raw["parameter"].duplicated().any():
    data_issues.append({"issue": "duplicated parameter rows", "handling": "first occurrence used"})

n = int(d["n"]); rho = float(d["rho"]); sigma = float(d["sigma"])
intercept = d["intercept"].lower() == "yes"
init = d["initial_condition"].lower()
R = int(d["n_replications"]); seed = int(d["seed"])
assert intercept and init == "stationary" and abs(rho) < 1 and n > 3 and sigma > 0

# --- baseline / analytical quantities ---
bias_analytical = -(1.0 + 3.0 * rho) / n
se_asymptotic = float(np.sqrt((1.0 - rho ** 2) / n))

# --- Monte Carlo ---
# One (R, n) standard-normal draw: column 0 -> stationary y_1, columns 1..n-1 -> innovations.
rng = np.random.default_rng(seed)
z = rng.standard_normal((R, n))
y = np.empty((R, n))
y[:, 0] = sigma / np.sqrt(1.0 - rho ** 2) * z[:, 0]
eps = sigma * z[:, 1:]
for t in range(1, n):
    y[:, t] = rho * y[:, t - 1] + eps[:, t - 1]   # c = 0
x = y[:, :-1]; yy = y[:, 1:]                       # t = 2..n  (n-1 obs)
xc = x - x.mean(axis=1, keepdims=True)
yc = yy - yy.mean(axis=1, keepdims=True)
sxx = (xc * xc).sum(axis=1)
rho_hat = (xc * yc).sum(axis=1) / sxx              # OLS slope with intercept (FWL)

ok = np.isfinite(rho_hat) & (sxx > 0)
n_bad = int((~ok).sum())
rows = [{"step": 1, "rule": "drop replications with non-finite rho_hat or zero regressor variance",
         "n_before": R, "n_after": R - n_bad, "n_removed": n_bad,
         "rationale": "OLS undefined; none expected under continuous Gaussian DGP"}]
pd.DataFrame(rows).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)
if n_bad:
    data_issues.append({"issue": f"{n_bad} degenerate replications", "handling": "excluded"})
rh = rho_hat[ok]
R_used = int(rh.size)

rho_hat_mean = float(rh.mean())
rho_hat_sd = float(rh.std(ddof=1))
bias_mc = rho_hat_mean - rho
bias_mc_se = rho_hat_sd / np.sqrt(R_used)
disc = (bias_mc - bias_analytical) / bias_mc_se

res = {
    "task_id": "F5-conf-02",
    "bias_analytical": float(bias_analytical),
    "bias_mc": float(bias_mc),
    "bias_mc_se": float(bias_mc_se),
    "rho_hat_mean": rho_hat_mean,
    "rho_hat_sd_mc": rho_hat_sd,
    "discrepancy_se_units": float(disc),
    "se_asymptotic": se_asymptotic,
    "n_replications": int(R),
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": data_issues,
    "n_input_rows": {"design.csv": n_rows},
    "n_analysed": R_used,
    "random_seed": seed,
    # supplementary diagnostics
    "supp_bias_analytical_neff": float(-(1.0 + 3.0 * rho) / (n - 1)),
    "supp_rho_hat_median": float(np.median(rh)),
    "supp_bias_mc_ci95_low": float(bias_mc - 1.96 * bias_mc_se),
    "supp_bias_mc_ci95_high": float(bias_mc + 1.96 * bias_mc_se),
    "supp_ratio_bias_mc_to_analytical": float(bias_mc / bias_analytical),
    "supp_ratio_sd_mc_to_se_asymptotic": float(rho_hat_sd / se_asymptotic),
    "supp_bias_mc_over_sd": float(bias_mc / rho_hat_sd),
    "supp_rmse_mc": float(np.sqrt(np.mean((rh - rho) ** 2))),
    "supp_share_rho_hat_below_rho": float(np.mean(rh < rho)),
    "design_sha256": hashlib.sha256(open(DATA, "rb").read()).hexdigest(),
}
with open(os.path.join(OUT, "results.json"), "w") as f:
    json.dump(res, f, indent=2)
print(json.dumps({k: res[k] for k in ["bias_analytical", "bias_mc", "bias_mc_se", "discrepancy_se_units",
                                      "rho_hat_mean", "rho_hat_sd_mc", "se_asymptotic", "n_analysed"]}, indent=2))