Repository Public Research preview
1"""Replication against rugarch on a shared synthetic fixture: the log-likelihood maximised by our implementations2(arch-backed or custom) must be >= rugarch's minus a small tolerance, and the fitted GARCH(1,1) parameters must agree.3Skipped when Rscript/rugarch are unavailable. Results are written to results/metrics/rugarch_crosscheck.json."""4import json5import shutil6import subprocess7from pathlib import Path89import numpy as np10import pandas as pd11import pytest1213from wp46.models.arch_models import ArchModel14from wp46.models.custom import CGARCH11, IGARCH11, NAGARCH11, RealGARCH111516ROOT = Path(__file__).resolve().parents[2]17HERE = Path(__file__).resolve().parent181920def _fixture():21 rng = np.random.default_rng(2024); T = 300022 nu = 7.0; om, al, ga, be = 0.03, 0.04, 0.08, 0.8823 z = rng.standard_t(nu, T) * np.sqrt((nu - 2) / nu); s = np.empty(T); e = np.empty(T); s[0] = om / (1 - al - be - ga / 2)24 for t in range(T):25 if t > 0:26 s[t] = om + al * e[t - 1] ** 2 + ga * e[t - 1] ** 2 * (e[t - 1] < 0) + be * s[t - 1]27 e[t] = np.sqrt(s[t]) * z[t]28 r = 0.02 + e29 x = s * rng.lognormal(0, 0.35, T)30 return r, x313233@pytest.mark.r34def test_rugarch_crosscheck():35 if shutil.which("Rscript") is None:36 pytest.skip("Rscript not available")37 r, x = _fixture()38 fx = HERE / "fixture_returns.csv"39 pd.DataFrame({"r": r, "x": x}).to_csv(fx, index=False)40 out = ROOT / "results" / "metrics" / "rugarch_crosscheck.json"41 out.parent.mkdir(parents=True, exist_ok=True)42 p = subprocess.run(["Rscript", str(HERE / "rugarch_crosscheck.R"), str(fx), str(out)], capture_output=True, text=True, timeout=1200)43 if p.returncode != 0:44 pytest.skip("rugarch run failed: " + p.stderr[-400:])45 R = json.loads(out.read_text())46 ours = {}47 for mid, dist in [("garch11", "t"), ("gjr11", "t"), ("egarch11", "t"), ("aparch11", "t"), ("figarch11", "t"), ("garch11", "skewt")]:48 f = ArchModel(mid).fit(r, dist=dist)49 ours[f"{mid}_{dist}"] = dict(loglik=f.loglik, params=dict(zip(f.param_names, map(float, f.params))), converged=f.converged)50 for cls, key, dist, xx in [(IGARCH11, "igarch11_t", "t", None), (CGARCH11, "cgarch11_t", "t", None), (NAGARCH11, "nagarch11_t", "t", None), (RealGARCH11, "realgarch11_n", "normal", x)]:51 f = cls().fit(r, xx, dist=dist)52 ours[key] = dict(loglik=f.loglik, params=dict(zip(f.param_names, map(float, f.params))), converged=f.converged)53 comp = {}54 for key in ["garch11_t", "gjr11_t", "egarch11_t", "aparch11_t", "figarch11_t", "igarch11_t", "cgarch11_t", "nagarch11_t", "garch11_skewt", "realgarch11_n"]:55 rk = R.get(key, {})56 ok = ours.get(key, {})57 comp[key] = dict(rugarch_loglik=rk.get("loglik"), ours_loglik=ok.get("loglik"), rugarch_coef=rk.get("coef"), ours_params=ok.get("params"),58 diff=(ok.get("loglik", np.nan) - rk.get("loglik", np.nan)) if rk.get("converged") else None)59 out.write_text(json.dumps({"rugarch": R, "ours": ours, "comparison": comp}, indent=1))60 # GARCH(1,1)-t: parameters and log-likelihood agree61 g = R["garch11_t"]62 assert g["converged"]63 assert abs(ours["garch11_t"]["loglik"] - g["loglik"]) < 0.564 assert abs(ours["garch11_t"]["params"]["alpha[1]"] - g["coef"]["alpha1"]) < 5e-365 assert abs(ours["garch11_t"]["params"]["beta[1]"] - g["coef"]["beta1"]) < 5e-366 # the remaining models: our maximised log-likelihood is not materially below rugarch's (same model, same data)67 for key in ["gjr11_t", "igarch11_t", "cgarch11_t", "nagarch11_t"]:68 if R[key].get("converged"):69 assert ours[key]["loglik"] >= R[key]["loglik"] - 1.0, key707172@pytest.mark.r73def test_extended_crosschecks():74 """TGARCH (arch) vs rugarch fGARCH/TGARCH; GARCH-X return equation vs rugarch external regressor; GAS-t vs R GAS;75 MS-GARCH likelihood and 22-step forecasts vs R MSGARCH. Each block is skipped if the R package did not run."""76 out = ROOT / "results" / "metrics" / "rugarch_crosscheck.json"77 if not out.exists():78 pytest.skip("run test_rugarch_crosscheck first")79 R = json.loads(out.read_text())80 r, x = _fixture()81 ours = R.setdefault("ours", {})82 # TGARCH83 f = ArchModel("tarch11").fit(r, dist="t")84 ours["tgarch11_t"] = dict(loglik=f.loglik, params=dict(zip(f.param_names, map(float, f.params))))85 if R["rugarch"].get("tgarch11_t", {}).get("converged"):86 assert abs(f.loglik - R["rugarch"]["tgarch11_t"]["loglik"]) < 1.087 # GARCH-X return equation only (our joint likelihood includes the RM equation; recompute the return part)88 from wp46.models.custom import GARCHXRV11, GAST11, MSGARCH289 from wp46.models.distributions import loglik_z90 fx = GARCHXRV11().fit(r, x, dist="normal")91 s2 = GARCHXRV11().filter(fx, r, x)92 ll_r = float(np.sum(loglik_z("normal", (r - fx.params[0]) / np.sqrt(s2), np.empty(0)) - 0.5 * np.log(s2)))93 ours["garchx_return_eq_n"] = dict(loglik_return_eq=ll_r, params=dict(zip(fx.param_names, map(float, fx.params))))94 rg = R["rugarch"].get("garchx_return_eq_n", {})95 if rg.get("converged"):96 # rugarch drops the first observation (lagged regressor) and may stop at gamma = 0 (nested in our model):97 # our return-equation likelihood per observation must not be below rugarch's98 assert ll_r / len(r) >= rg["loglik"] / (len(r) - 1) - 0.005, (ll_r / len(r), rg["loglik"] / (len(r) - 1))99 # GAS-t100 fg = GAST11().fit(r, dist="t")101 ours["gas_t"] = dict(loglik=fg.loglik, params=dict(zip(fg.param_names, map(float, fg.params))))102 rg = R["rugarch"].get("gas_t", {})103 if rg.get("converged") and np.size(rg.get("loglik")) > 0:104 ours["gas_t"]["r_gas_loglik_reference"] = float(np.ravel(rg["loglik"])[0]) # different parameterization: reference only105 # MS-GARCH106 fm = MSGARCH2().fit(r, dist="normal")107 fc = MSGARCH2().forecast(fm, r, None, 22)108 ours["msgarch2_n"] = dict(loglik=fm.loglik, params=dict(zip(fm.param_names, map(float, fm.params))), var_forecast=[float(v) for v in fc.step_var])109 rg = R["rugarch"].get("msgarch2_n", {})110 if rg.get("converged"):111 rl = float(np.ravel(rg["loglik"])[0])112 assert fm.loglik >= rl - 2.0, (fm.loglik, rl)113 vf = np.asarray(rg["vol_forecast"]) ** 2114 rel = np.abs(fc.step_var - vf) / vf115 assert np.max(rel) < 0.10, rel # same multi-step recursion up to the parameter estimates116 out.write_text(json.dumps(R, indent=1))117