1"""F5-conf-02: Monte Carlo finite-sample bias of OLS AR(1) estimator vs Marriott-Pope/Kendall approximation."""2import json, os3import numpy as np, pandas as pd45ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))6d = pd.read_csv(os.path.join(ROOT, "data", "design.csv"))7p = dict(zip(d["parameter"], d["value"]))8issues = []9n = int(float(p["n"])); rho = float(p["rho"]); sigma = float(p["sigma"])10R = int(float(p["n_replications"])); seed = int(float(p["seed"]))11assert str(p["intercept"]).strip().lower() == "yes"12assert str(p["initial_condition"]).strip().lower() == "stationary"13assert abs(rho) < 1 and sigma > 0 and n > 3 and R > 11415rng = np.random.default_rng(seed)16# Draw order: y_1 stationary draws (R), then innovations (R, n-1)17y = np.empty((R, n))18y[:, 0] = rng.normal(0.0, sigma / np.sqrt(1 - rho**2), size=R)19eps = rng.normal(0.0, sigma, size=(R, n - 1))20for t in range(1, n):21 y[:, t] = rho * y[:, t - 1] + eps[:, t - 1] # c = 02223x = y[:, :-1]; z = y[:, 1:] # t = 2..n24xc = x - x.mean(axis=1, keepdims=True)25zc = z - z.mean(axis=1, keepdims=True)26rho_hat = (xc * zc).sum(1) / (xc * xc).sum(1) # OLS slope with intercept27ok = np.isfinite(rho_hat)28if not ok.all():29 issues.append({"issue": f"{(~ok).sum()} non-finite estimates", "handling": "excluded"})30rh = rho_hat[ok]; Ru = rh.size3132mean = rh.mean(); sd = rh.std(ddof=1)33bias_mc = mean - rho; se = sd / np.sqrt(Ru)34bias_an = -(1 + 3 * rho) / n35res = {36 "task_id": "F5-conf-02",37 "bias_analytical": bias_an,38 "bias_mc": bias_mc,39 "bias_mc_se": se,40 "rho_hat_mean": mean,41 "rho_hat_sd_mc": sd,42 "discrepancy_se_units": (bias_mc - bias_an) / se,43 "se_asymptotic": float(np.sqrt((1 - rho**2) / n)),44 "n_replications": R,45 "feasibility": "complete",46 "feasibility_reason": None,47 "data_issues": issues,48 "n_input_rows": {"design.csv": int(len(d))},49 "n_analysed": int(Ru),50 "random_seed": seed,51 # supplementary diagnostics52 "supplementary": {53 "bias_analytical_using_n_minus_1": -(1 + 3 * rho) / (n - 1),54 "rho_hat_median": float(np.median(rh)),55 "share_rho_hat_below_rho": float((rh < rho).mean()),56 "ratio_bias_mc_to_analytical": bias_mc / bias_an,57 "ratio_sd_mc_to_se_asymptotic": sd / np.sqrt((1 - rho**2) / n),58 },59}60res = {k: (float(v) if isinstance(v, np.floating) else v) for k, v in res.items()}61res["supplementary"] = {k: float(v) for k, v in res["supplementary"].items()}62with open(os.path.join(ROOT, "out", "results.json"), "w") as f:63 json.dump(res, f, indent=2)64print(json.dumps(res, indent=2))65