1"""F5-conf-02: finite-sample bias of OLS AR(1) coefficient. Regenerates out/results.json from data/ only."""2import json, os, hashlib3import numpy as np4import pandas as pd56ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))7DATA = os.path.join(ROOT, "data", "design.csv")8OUT = os.path.join(ROOT, "out")910df = pd.read_csv(DATA)11d = dict(zip(df["parameter"].astype(str).str.strip(), df["value"].astype(str).str.strip()))12issues = []13required = ["n", "rho", "sigma", "intercept", "initial_condition", "n_replications", "seed"]14missing = [k for k in required if k not in d]15if missing:16 raise SystemExit(f"design missing parameters: {missing}")17if df["parameter"].duplicated().any():18 issues.append({"issue": "duplicated parameter rows", "handling": "first occurrence used"})1920n = int(d["n"]); rho = float(d["rho"]); sigma = float(d["sigma"])21R = int(d["n_replications"]); seed = int(d["seed"])22assert d["intercept"].lower() == "yes" and d["initial_condition"].lower() == "stationary"23assert abs(rho) < 1 and sigma > 0 and n >= 3 and R >= 22425# ---- Simulation (c = 0). Draw order: y_1 for all reps, then innovations (R x (n-1)).26rng = np.random.default_rng(seed)27y = np.empty((R, n))28y[:, 0] = rng.normal(0.0, sigma / np.sqrt(1.0 - rho**2), size=R)29eps = rng.normal(0.0, sigma, size=(R, n - 1))30for t in range(1, n):31 y[:, t] = rho * y[:, t - 1] + eps[:, t - 1]3233# ---- OLS of y_t on (1, y_{t-1}), t = 2..n (slope = demeaned cross-product ratio)34x = y[:, :-1]; z = y[:, 1:]35xc = x - x.mean(axis=1, keepdims=True); zc = z - z.mean(axis=1, keepdims=True)36rho_hat = (xc * zc).sum(axis=1) / (xc * xc).sum(axis=1)3738# cross-check one replication with lstsq39X0 = np.column_stack([np.ones(n - 1), x[0]])40b0 = np.linalg.lstsq(X0, z[0], rcond=None)[0]41assert abs(b0[1] - rho_hat[0]) < 1e-104243finite = np.isfinite(rho_hat)44n_excl = int((~finite).sum())45rho_hat_used = rho_hat[finite]46Rused = rho_hat_used.size4748mean_ = float(rho_hat_used.mean())49sd_ = float(rho_hat_used.std(ddof=1))50bias_mc = mean_ - rho51bias_mc_se = sd_ / np.sqrt(Rused)52bias_an = -(1 + 3 * rho) / n53se_asy = float(np.sqrt((1 - rho**2) / n))5455# supplementary (not required keys): formula with effective sample size n-1, median bias56bias_an_nm1 = -(1 + 3 * rho) / (n - 1)5758results = {59 "task_id": "F5-conf-02",60 "bias_analytical": bias_an,61 "bias_mc": bias_mc,62 "bias_mc_se": float(bias_mc_se),63 "rho_hat_mean": mean_,64 "rho_hat_sd_mc": sd_,65 "discrepancy_se_units": float((bias_mc - bias_an) / bias_mc_se),66 "se_asymptotic": se_asy,67 "n_replications": int(R),68 "feasibility": "complete",69 "feasibility_reason": None,70 "data_issues": issues,71 "n_input_rows": {"design.csv": int(len(df))},72 "n_analysed": int(Rused),73 "random_seed": seed,74 # supplementary diagnostics75 "supp_bias_analytical_n_minus_1": bias_an_nm1,76 "supp_discrepancy_se_units_n_minus_1": float((bias_mc - bias_an_nm1) / bias_mc_se),77 "supp_bias_mc_ci95_lower": float(bias_mc - 1.96 * bias_mc_se),78 "supp_bias_mc_ci95_upper": float(bias_mc + 1.96 * bias_mc_se),79 "supp_relative_error_analytical": float((bias_an - bias_mc) / bias_mc),80 "supp_rho_hat_median": float(np.median(rho_hat_used)),81 "supp_frac_rho_hat_ge_1": float((rho_hat_used >= 1).mean()),82 "supp_sd_ratio_mc_to_asymptotic": float(sd_ / se_asy),83 "supp_n_nonfinite_excluded": n_excl,84 "supp_design_sha256": hashlib.sha256(open(DATA, "rb").read()).hexdigest(),85}8687os.makedirs(OUT, exist_ok=True)88with open(os.path.join(OUT, "results.json"), "w") as f:89 json.dump(results, f, indent=2)90pd.DataFrame([91 {"step": "S0", "rule": "Monte Carlo replications simulated", "n_before": R, "n_after": R, "n_removed": 0,92 "rationale": "all R replications from design"},93 {"step": "S1", "rule": "drop replications with non-finite rho_hat", "n_before": R, "n_after": Rused,94 "n_removed": n_excl, "rationale": "degenerate regressor would make OLS undefined"},95]).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)96print(json.dumps({k: results[k] for k in ["bias_analytical", "bias_mc", "bias_mc_se", "discrepancy_se_units"]}, indent=2))97