#!/usr/bin/env python3
"""Generates the synthetic 1-minute OHLCV fixtures used by CI (no vendor data) and the reference values
(numpy/statsmodels/exchange_calendars/arch) consumed by tests/reference. Deterministic (seed 20261102).

Fixture: 2 symbols × 5 sessions around Thanksgiving 2024 (2024-11-25 … 2024-11-29, incl. the early close on 11-29
and the holiday 11-28), extended hours 04:00–19:59 ET with regular bars 09:30–16:00 inclusive; ~3 % of regular
minutes randomly absent (zero-volume minutes are absent at the source); the 16:00 bar is present on every session.
"""
import json, math, os, sys
import numpy as np
import pandas as pd

root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
out_dir = os.path.join(root, "datasets", "fixtures"); os.makedirs(out_dir, exist_ok=True)
rng = np.random.default_rng(20261102)
sessions = ["2024-11-25", "2024-11-26", "2024-11-27", "2024-11-29", "2024-12-02"]
early = {"2024-11-29"}
rows = []
truth = {}
for sym, p0, vol in [("SYN1", 100.0, 0.0008), ("SYN2", 25.0, 0.0015)]:
    price = p0
    truth[sym] = {}
    for s in sessions:
        close_min = 13 * 60 if s in early else 16 * 60
        ext_end = 17 * 60 if s in early else 20 * 60
        minutes = list(range(4 * 60, ext_end))
        day = pd.Timestamp(s, tz="America/New_York")
        sess_rv = 0.0; n_reg = 0; closes = []
        for m in minutes:
            regular = 9 * 60 + 30 <= m <= close_min
            # drop ~3 % of regular minutes and ~40 % of extended minutes, but never the open/close prints
            if m not in (9 * 60 + 30, close_min):
                if regular and rng.random() < 0.03: continue
                if not regular and rng.random() < 0.40: continue
            r = rng.normal(0, vol * (1.0 if regular else 0.5))
            o = price; c = price * math.exp(r); h = max(o, c) * (1 + abs(rng.normal(0, 1e-4))); l = min(o, c) * (1 - abs(rng.normal(0, 1e-4)))
            v = float(rng.integers(100, 5000)) if regular else float(rng.integers(0, 300) + 1)
            ts = (day + pd.Timedelta(minutes=m)).tz_convert("UTC")
            rows.append((sym, ts.strftime("%Y-%m-%dT%H:%M:%SZ"), round(o, 4), round(h, 4), round(l, 4), round(c, 4), v))
            price = c
            if regular:
                closes.append(round(c, 4)); n_reg += 1   # truth from the rounded prices written to the CSV
        lr = np.diff(np.log(np.array(closes)))
        truth[sym][s] = {"n_regular_bars": n_reg, "rv": float(np.sum(lr ** 2)), "bpv": float((math.pi / 2) * (len(lr) / (len(lr) - 1)) * np.sum(np.abs(lr[1:]) * np.abs(lr[:-1]))),
                         "rsv_neg": float(np.sum(lr[lr < 0] ** 2)), "rsv_pos": float(np.sum(lr[lr > 0] ** 2)), "rq": float(len(lr) / 3 * np.sum(lr ** 4)),
                         "expected_bars": 211 if s in early else 391}
df = pd.DataFrame(rows, columns=["symbol", "timestamp_utc", "open", "high", "low", "close", "volume"])
df.to_csv(os.path.join(out_dir, "synthetic_bars_1m.csv"), index=False)
json.dump(truth, open(os.path.join(out_dir, "synthetic_bars_1m.expected.json"), "w"), indent=1)
print("wrote", len(df), "bars for", df.symbol.nunique(), "symbols;", df.groupby("symbol").size().to_dict())

# --- reference values from established software -------------------------------------------------------------
ref = {}
try:
    import exchange_calendars as xc
    cal = xc.get_calendar("XNYS", start="2000-01-03", end="2026-12-31")
    sess = cal.sessions_in_range("2000-01-03", "2026-12-31")
    ref["xnys"] = {"n_sessions_2000_2026": int(len(sess)), "sessions_2024": [str(d.date()) for d in cal.sessions_in_range("2024-01-01", "2024-12-31")],
                   "early_closes_2000_2026": [str(d.date()) for d in cal.early_closes if "2000-01-03" <= str(d.date()) <= "2026-12-31"],
                   "holidays_2000_2026": [str(d.date()) for d in pd.bdate_range("2000-01-03", "2026-12-31") if d not in sess],
                   "exchange_calendars_version": xc.__version__}
except Exception as e:
    ref["xnys"] = {"error": str(e)}
import statsmodels.api as sm
from scipy import stats as sps
X = np.column_stack([np.ones(500), rng.normal(size=500), rng.normal(size=500) * 2 + 1]); y = 0.5 + 1.5 * X[:, 1] - 0.7 * X[:, 2] + rng.normal(size=500) * (1 + 0.5 * np.abs(X[:, 1]))
ols = sm.OLS(y, X).fit(); hc0 = sm.OLS(y, X).fit(cov_type="HC0"); hc3 = sm.OLS(y, X).fit(cov_type="HC3"); hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": 5, "use_correction": False})
groups = np.arange(500) // 25; cl = sm.OLS(y, X).fit(cov_type="cluster", cov_kwds={"groups": groups})
np.savetxt(os.path.join(out_dir, "ols_design.csv"), np.column_stack([y, X[:, 1], X[:, 2]]), delimiter=",", header="y,x1,x2", comments="", fmt="%.17g")
ref["ols"] = {"params": ols.params.tolist(), "bse": ols.bse.tolist(), "pvalues": ols.pvalues.tolist(), "rsquared": float(ols.rsquared), "ssr": float(ols.ssr),
              "hc0_bse": hc0.bse.tolist(), "hc3_bse": hc3.bse.tolist(), "hac5_bse": hac.bse.tolist(), "cluster_bse": cl.bse.tolist(), "statsmodels_version": sm.__version__,
              "t_pvalues": {"1.96_inf": float(2 * sps.norm.sf(1.96)), "2.0_10": float(2 * sps.t.sf(2.0, 10)), "-3.0_30": float(2 * sps.t.sf(3.0, 30))}}
try:
    from arch import arch_model
    import arch
    s2 = 0.04 / (1 - 0.05 - 0.9); e = 0.0; r = np.empty(3000); g = np.random.default_rng(7)
    for t in range(3000):
        s2 = 0.04 + 0.05 * e * e + 0.9 * s2; e = math.sqrt(s2) * g.standard_normal(); r[t] = 0.01 + e
    np.savetxt(os.path.join(out_dir, "garch_returns.csv"), r, header="r", comments="", fmt="%.17g")
    out = {"arch_version": arch.__version__}
    for name, kw in [("garch", dict(vol="GARCH", p=1, q=1)), ("gjr", dict(vol="GARCH", p=1, o=1, q=1)), ("egarch", dict(vol="EGARCH", p=1, o=1, q=1))]:
        for dist in ["normal", "t"]:
            am = arch_model(r, mean="Constant", dist=dist, rescale=False, **kw); res = am.fit(disp="off")
            out[f"{name}.{dist}"] = {"params": {k: float(v) for k, v in res.params.items()}, "std_err": {k: float(v) for k, v in res.std_err.items()}, "loglik": float(res.loglikelihood),
                                     "converged": bool(res.convergence_flag == 0), "forecast_h5": res.forecast(horizon=5, reindex=False).variance.values[-1].tolist() if name != "egarch" else res.forecast(horizon=5, method="simulation", simulations=20000, reindex=False).variance.values[-1].tolist()}
    ref["garch"] = out
except Exception as e:
    ref["garch"] = {"error": repr(e)}
json.dump(ref, open(os.path.join(out_dir, "reference_values.json"), "w"), indent=1)
print("reference values written:", list(ref.keys()))

# --- synthetic Arrow IPC page + Parquet in the FinData bars schema v1 (one regular session, 391 bars) --------------
import pyarrow as pa, pyarrow.ipc as ipc, pyarrow.parquet as pq
day = pd.Timestamp("2024-11-27", tz="America/New_York"); n = 391
ts = [(day + pd.Timedelta(minutes=9 * 60 + 30 + i)).tz_convert("UTC").to_pydatetime() for i in range(n)]
g = np.random.default_rng(3); c = 100 * np.exp(np.cumsum(g.normal(0, 5e-4, n))); o = np.r_[100.0, c[:-1]]
schema = pa.schema([("instrument_id", pa.string()), ("symbol_as_reported", pa.string()), ("venue", pa.string()), ("asset_class", pa.string()), ("interval", pa.string()),
                    ("timestamp_utc", pa.timestamp("us", tz="UTC")), ("session_date", pa.date32()), ("open", pa.float64()), ("high", pa.float64()), ("low", pa.float64()), ("close", pa.float64()),
                    ("volume", pa.float64()), ("currency", pa.string()), ("trade_count", pa.int32()), ("vwap", pa.float64()), ("adjustment_type", pa.string()), ("source_id", pa.string()),
                    ("ingested_at", pa.timestamp("us", tz="UTC")), ("quality_flags", pa.list_(pa.string())), ("revision_id", pa.int32())],
                   metadata={"findata.schema": "bars", "findata.schema_version": "1", "findata.timestamp_convention": "bar_open"})
tbl = pa.table({"instrument_id": ["synthetic:stock:SYN1"] * n, "symbol_as_reported": ["SYN1"] * n, "venue": [None] * n, "asset_class": ["stock"] * n, "interval": ["1m"] * n,
                "timestamp_utc": pa.array(ts, pa.timestamp("us", tz="UTC")), "session_date": [day.date()] * n, "open": o, "high": np.maximum(o, c) * 1.0001, "low": np.minimum(o, c) * 0.9999, "close": c,
                "volume": g.integers(100, 1000, n).astype(float), "currency": ["USD"] * n, "trade_count": pa.array([None] * n, pa.int32()), "vwap": pa.array([None] * n, pa.float64()),
                "adjustment_type": ["raw"] * n, "source_id": ["synthetic"] * n, "ingested_at": pa.array([ts[0]] * n, pa.timestamp("us", tz="UTC")), "quality_flags": [[]] * n, "revision_id": pa.array([1] * n, pa.int32())}, schema=schema)
with open(os.path.join(out_dir, "synthetic_bars_page.arrow"), "wb") as fh:
    with ipc.new_stream(fh, schema) as w:
        w.write_table(tbl)
pq.write_table(tbl, os.path.join(out_dir, "synthetic_bars_page.parquet"))
print("wrote synthetic_bars_page.arrow/.parquet", tbl.num_rows, "rows")

# --- HAR reference: simulated daily RV (log-AR(1)), HAR-RV h=1 and h=5 (average target) with HAC, log-HAR, HARQ -----
g = np.random.default_rng(40); T = 1500; lv = np.empty(T); lv[0] = np.log(1e-4)
for t in range(1, T): lv[t] = 0.3 * np.log(1e-4) + 0.7 * lv[t - 1] + 0.4 * g.standard_normal()
rv = np.exp(lv); rq = rv ** 2 * (1 + 0.5 * g.random(T))
np.savetxt(os.path.join(out_dir, "har_rv.csv"), np.column_stack([rv, rq]), delimiter=",", header="rv,rq", comments="", fmt="%.17g")
def har_design(x, lags=(1, 5, 22), h=1, log=False, harq=None, target="average"):
    z = np.log(x) if log else x; L = max(lags); rows = []; ys = []
    sqrt_rq_mean = np.sqrt(np.mean(np.maximum(harq, 0.0))) if harq is not None else 0.0
    for t in range(L - 1, len(x) - h):
        row = [1.0] + [z[t - l + 1:t + 1].mean() for l in lags]
        if harq is not None: row.append((np.sqrt(max(harq[t], 0.0)) - sqrt_rq_mean) * x[t])
        if target == "average": yv = x[t + 1:t + h + 1].mean(); yv = np.log(yv) if log else yv
        else: yv = z[t + h]
        rows.append(row); ys.append(yv)
    return np.array(rows), np.array(ys)
har = {}
for name, kw in [("h1", dict(h=1)), ("h5", dict(h=5)), ("h22", dict(h=22)), ("log_h1", dict(h=1, log=True)), ("harq_h1", dict(h=1, harq=rq))]:
    X, y = har_design(rv, **kw); n = len(y); h = kw["h"]; L = max(h - 1, int(np.floor(4 * (n / 100) ** (2 / 9))))
    f = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": False})
    har[name] = {"params": f.params.tolist(), "hac_bse": f.bse.tolist(), "hac_lags": L, "n": n, "ssr": float(f.ssr), "rsquared": float(f.rsquared)}
ref["har"] = har; json.dump(ref, open(os.path.join(out_dir, "reference_values.json"), "w"), indent=1)
print("HAR references:", {k: v["n"] for k, v in har.items()})

# --- VAR(2) reference: statsmodels VAR, orthogonalised FEVD, generalized FEVD (Pesaran–Shin) by hand, DY totals, LW ----
from statsmodels.tsa.api import VAR
g = np.random.default_rng(77); N = 4; T = 1200
A1 = np.array([[0.5, 0.1, 0.0, 0.05], [0.0, 0.4, 0.15, 0.0], [0.1, 0.0, 0.3, 0.1], [0.05, 0.05, 0.0, 0.45]]); A2 = 0.15 * np.eye(N)
L = np.linalg.cholesky(np.array([[1.0, 0.3, 0.1, 0.0], [0.3, 1.0, 0.2, 0.1], [0.1, 0.2, 1.0, 0.3], [0.0, 0.1, 0.3, 1.0]]))
Y = np.zeros((T, N))
for t in range(2, T): Y[t] = 0.1 + A1 @ Y[t - 1] + A2 @ Y[t - 2] + L @ g.standard_normal(N)
np.savetxt(os.path.join(out_dir, "var_panel.csv"), Y, delimiter=",", header=",".join(f"y{i}" for i in range(N)), comments="", fmt="%.17g")
res = VAR(Y).fit(2, trend="c"); H = 10
Phi = [np.eye(N)]
coefs = [res.coefs[i] for i in range(2)]
for h in range(1, H): Phi.append(sum(coefs[j - 1] @ Phi[h - j] for j in range(1, min(h, 2) + 1)))
S = res.sigma_u
num = np.zeros((N, N)); den = np.zeros(N)
for h in range(H):
    PS = Phi[h] @ S
    for i in range(N):
        den[i] += PS[i] @ Phi[h][i]
        for j in range(N): num[i, j] += PS[i, j] ** 2 / S[j, j]
th = num / den[:, None]; th = th / th.sum(axis=1, keepdims=True)
dy_total = 100 * (th.sum() - np.trace(th)) / N; dy_from = 100 * (th.sum(axis=1) - np.diag(th)); dy_to = 100 * (th.sum(axis=0) - np.diag(th))
sel = VAR(Y).select_order(4)
try:
    from sklearn.covariance import ledoit_wolf as sk_lw
    lw_cov, lw_s = sk_lw(Y)
except Exception:
    Xc = Y - Y.mean(axis=0); n = Xc.shape[0]; Sc = Xc.T @ Xc / n; mu = np.trace(Sc) / N
    d2 = np.sum((Sc - mu * np.eye(N)) ** 2) / N; b = sum(np.sum((np.outer(x, x) - Sc) ** 2) for x in Xc) / (n * n * N); lw_s = min(1.0, b / d2); lw_cov = (1 - lw_s) * Sc + lw_s * mu * np.eye(N)
ref = json.load(open(os.path.join(out_dir, "reference_values.json")))
ref["var"] = {"nobs": int(res.nobs), "params": res.params.tolist(), "sigma_u": res.sigma_u.tolist(), "llf": float(res.llf), "fevd10": res.fevd(H).decomp[:, -1, :].tolist(),
              "gfevd10": th.tolist(), "dy_total10": float(dy_total), "dy_from10": dy_from.tolist(), "dy_to10": dy_to.tolist(), "bic_order": int(sel.bic), "aic_order": int(sel.aic),
              "lw_shrinkage": float(lw_s), "lw_cov": np.asarray(lw_cov).tolist(), "statsmodels_version": sm.__version__}
json.dump(ref, open(os.path.join(out_dir, "reference_values.json"), "w"), indent=1)
print("VAR reference: nobs", res.nobs, "DY total", round(dy_total, 3), "bic order", int(sel.bic), "LW shrinkage", round(float(lw_s), 4))