"""F5-conf-02: Monte Carlo bias of OLS AR(1) coefficient vs Marriott-Pope/Kendall approximation.
Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""
import json, os
import numpy as np
import pandas as pd

d = pd.read_csv("data/design.csv")
p = dict(zip(d["parameter"], d["value"]))
n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])
R = int(p["n_replications"]); seed = int(p["seed"])
assert str(p["intercept"]).strip().lower() == "yes"
assert str(p["initial_condition"]).strip().lower() == "stationary"
c = 0.0

rng = np.random.default_rng(seed)
# One (R, n) standard-normal draw: column 0 -> stationary y_1, columns 1..n-1 -> innovations.
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] = c + rho * y[:, t - 1] + eps[:, t - 1]

# OLS of y_t on (1, y_{t-1}), t = 2..n  -> slope = cov/var (demeaned within each sample)
x = y[:, :-1]; yy = y[:, 1:]
xd = x - x.mean(axis=1, keepdims=True); yd = yy - yy.mean(axis=1, keepdims=True)
rho_hat = (xd * yd).sum(axis=1) / (xd**2).sum(axis=1)

ok = np.isfinite(rho_hat)
rho_hat = rho_hat[ok]
m = float(rho_hat.mean()); sd = float(rho_hat.std(ddof=1))
bias_mc = m - rho; se = sd / np.sqrt(len(rho_hat))
bias_an = -(1 + 3 * rho) / n
res = {
    "task_id": "F5-conf-02",
    "bias_analytical": bias_an,
    "bias_mc": bias_mc,
    "bias_mc_se": se,
    "rho_hat_mean": m,
    "rho_hat_sd_mc": sd,
    "discrepancy_se_units": (bias_mc - bias_an) / se,
    "se_asymptotic": float(np.sqrt((1 - rho**2) / n)),
    "n_replications": R,
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": [],
    "n_input_rows": {"design.csv": int(len(d))},
    "n_analysed": int(ok.sum()),
    "random_seed": seed,
    # supplementary
    "bias_analytical_n_minus_1": -(1 + 3 * rho) / (n - 1),
    "bias_mc_ci95": [bias_mc - 1.96 * se, bias_mc + 1.96 * se],
    "rho_hat_median": float(np.median(rho_hat)),
    "share_rho_hat_below_rho": float((rho_hat < rho).mean()),
}
os.makedirs("out", exist_ok=True)
with open("out/results.json", "w") as f:
    json.dump(res, f, indent=2)
print(json.dumps(res, indent=2))