"""F5-conf-02: finite-sample bias of OLS rho_hat in a Gaussian AR(1) with intercept.
Run from workspace root: ~/wp29/venv/bin/python out/analysis.py
Reads only data/design.csv; writes out/results.json and out/exclusions.csv deterministically."""
import json, hashlib, os
import numpy as np
import pandas as pd

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

d = pd.read_csv(DESIGN)
p = dict(zip(d["parameter"].astype(str).str.strip(), d["value"].astype(str).str.strip()))
issues = []
n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])
R = int(p["n_replications"]); seed = int(p["seed"])
intercept = p["intercept"].lower() == "yes"
init = p["initial_condition"].lower()
assert intercept, "design requires intercept"
assert init == "stationary" and abs(rho) < 1, "stationary init requires |rho|<1"
assert n >= 3 and R >= 2 and sigma > 0

# ---- P2 baseline quantities (closed form) ----
bias_analytical = -(1 + 3 * rho) / n
se_asymptotic = float(np.sqrt((1 - rho ** 2) / n))
# sensitivity: same formula with effective number of regression observations n-1
bias_analytical_neff = -(1 + 3 * rho) / (n - 1)

# ---- P3 Monte Carlo ----
# Draw order: one (R, n) standard-normal matrix; column 0 scaled to the stationary
# sd for y_1, columns 1..n-1 scaled by sigma as innovations e_2..e_n. c = 0.
rng = np.random.default_rng(seed)
z = rng.standard_normal((R, n))
y = np.empty((R, n))
y[:, 0] = z[:, 0] * sigma / np.sqrt(1 - rho ** 2)
eps = sigma * z[:, 1:]
for t in range(1, n):
    y[:, t] = rho * y[:, t - 1] + eps[:, t - 1]
x = y[:, :-1]; yy = y[:, 1:]          # t = 2..n  (n-1 pairs)
xc = x - x.mean(axis=1, keepdims=True)
yc = yy - yy.mean(axis=1, keepdims=True)
rho_hat = (xc * yc).sum(axis=1) / (xc ** 2).sum(axis=1)

finite = np.isfinite(rho_hat)
n_bad = int((~finite).sum())
if n_bad:
    issues.append({"issue": f"{n_bad} non-finite rho_hat", "handling": "excluded; see exclusions.csv"})
rh = rho_hat[finite]
Rused = 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(Rused)
disc = (bias_mc - bias_analytical) / bias_mc_se
disc_neff = (bias_mc - bias_analytical_neff) / bias_mc_se
q = np.quantile(rh, [0.025, 0.5, 0.975])

pd.DataFrame([{"step": 1, "rule": "drop replications with non-finite rho_hat",
               "n_before": R, "n_after": Rused, "n_removed": R - Rused,
               "rationale": "OLS undefined if regressor has zero variance"}]
             ).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)

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": R,
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": issues,
    "n_input_rows": {"design.csv": int(len(d))},
    "n_analysed": Rused,
    "random_seed": seed,
    # supplementary keys
    "design_n": n, "design_rho": rho, "design_sigma": sigma,
    "bias_mc_ci95": [float(bias_mc - 1.96 * bias_mc_se), float(bias_mc + 1.96 * bias_mc_se)],
    "bias_analytical_neff": float(bias_analytical_neff),
    "discrepancy_se_units_neff": float(disc_neff),
    "bias_ratio_mc_to_analytical": float(bias_mc / bias_analytical),
    "rho_hat_median": float(q[1]),
    "rho_hat_q025": float(q[0]), "rho_hat_q975": float(q[2]),
    "sd_ratio_mc_to_asymptotic": float(rho_hat_sd / se_asymptotic),
    "share_rho_hat_below_rho": float((rh < rho).mean()),
    "design_sha256": hashlib.sha256(open(DESIGN, "rb").read()).hexdigest(),
}
with open(os.path.join(OUT, "results.json"), "w") as f:
    json.dump(res, f, indent=2)
print(json.dumps(res, indent=2))