1import numpy as np2import pytest3from arch.univariate import GARCH, ConstantMean, StudentsT45from wp46.models.custom import (CGARCH11, GARCH11, GAST11, GJR11, IGARCH11, MODELS, MSGARCH2, NAGARCH11, QGARCH11,6 HEAVY11, GARCHXRV11, RealGARCH11)7from wp46.models.distributions import loglik_z8910def simulate_garch(T=2500, seed=1, nu=6.0, om=0.05, al=0.08, be=0.90, mu=0.03):11 rng = np.random.default_rng(seed)12 z = rng.standard_t(nu, T) * np.sqrt((nu - 2) / nu)13 s = np.empty(T); e = np.empty(T); s[0] = om / (1 - al - be)14 for t in range(T):15 if t > 0:16 s[t] = om + al * e[t - 1] ** 2 + be * s[t - 1]17 e[t] = np.sqrt(s[t]) * z[t]18 return mu + e192021R = simulate_garch()222324def test_custom_garch_equals_arch_loglik_and_forecast():25 am = ConstantMean(R); am.volatility = GARCH(1, 0, 1); am.distribution = StudentsT()26 res = am.fit(disp="off")27 m = GARCH11(); f = m.fit(R, dist="t")28 assert abs(f.loglik - res.loglikelihood) < 0.0529 assert np.max(np.abs(f.params - res.params.values)) < 2e-330 f.params = res.params.values.copy(); f.aux["backcast"] = float(am.volatility.backcast(R - res.params["mu"]))31 s2 = m.filter(f, R)32 ref = np.empty(len(R)); resid = R - res.params["mu"]33 am.volatility.compute_variance(res.params.values[1:4], resid, ref, am.volatility.backcast(resid), am.volatility.variance_bounds(resid))34 assert np.max(np.abs(s2 - ref)) < 1e-9 # identical recursion and initialization (arch 8.0.0 compute_variance)35 fa = res.forecast(horizon=22, reindex=False).variance.values[0]36 fc = m.forecast(f, R, None, 22).step_var37 assert np.max(np.abs(fa - fc) / fa) < 1e-10383940@pytest.mark.parametrize("cls", [GJR11, NAGARCH11, QGARCH11, CGARCH11, IGARCH11, GAST11])41def test_nesting_and_admissibility(cls):42 f = cls().fit(R, dist="t")43 assert f.converged and f.admissible44 # all nest GARCH(1,1): their log-likelihood is at least that of GARCH(1,1) minus a tolerance (IGARCH/GAS are not nests)45 g = GARCH11().fit(R, dist="t")46 if cls not in (IGARCH11, GAST11):47 assert f.loglik >= g.loglik - 0.548 fc = cls().forecast(f, R, None, 22)49 assert fc.step_var.shape == (22,) and np.all(fc.step_var > 0) and np.all(np.isfinite(fc.step_var))505152def test_msgarch_two_regimes_ordered_and_forecast_positive():53 f = MSGARCH2().fit(R, dist="t")54 assert f.converged55 v = f.params56 u1 = v[3] / (1 - v[4] - v[5]); u2 = v[6] / (1 - v[7] - v[8])57 assert u1 <= u2 * 1.0158 fc = MSGARCH2().forecast(f, R, None, 22)59 assert np.all(fc.step_var > 0)606162def test_realized_models_fit_on_synthetic_measure():63 rng = np.random.default_rng(3)64 x = (R - R.mean()) ** 2 * rng.lognormal(0, 0.3, len(R)) + 0.0565 for cls in (RealGARCH11, HEAVY11, GARCHXRV11):66 f = cls().fit(R, x, dist="normal")67 assert f.converged, cls.model_id68 fc = cls().forecast(f, R, x, 22)69 assert np.all(fc.step_var > 0) and np.all(np.isfinite(fc.step_var))707172def test_no_lookahead_in_forecasts():73 """Perturbing returns AFTER the origin must not change the forecast made at the origin."""74 for mid in ("garch11_custom", "gjr11_custom", "nagarch11", "cgarch11", "msgarch2"):75 m = MODELS[mid]()76 f = m.fit(R[:2000], dist="t" if mid != "msgarch2" else "normal")77 a = m.forecast(f, R[:2000], None, 22).step_var78 R2 = R.copy(); R2[2000:] *= 5.079 b = m.forecast(f, R2[:2000], None, 22).step_var80 assert np.array_equal(a, b)818283def test_edge_cases_constant_and_extreme():84 const = np.full(600, 0.1)85 f = GARCH11().fit(const + np.random.default_rng(0).normal(0, 1e-6, 600), dist="normal")86 assert np.isfinite(f.loglik)87 spike = R[:1500].copy(); spike[700] = 60.088 f = GJR11().fit(spike, dist="t")89 assert f.converged and np.all(np.isfinite(GJR11().forecast(f, spike, None, 22).step_var))909192def test_garchm_nests_garch_and_recovers_delta():93 from wp46.models.custom import GARCHM1194 g = GARCH11().fit(R, dist="t"); m = GARCHM11().fit(R, dist="t")95 assert m.loglik >= g.loglik - 0.5 # GARCH-M nests GARCH(1,1) (delta = 0)96 rng = np.random.default_rng(21); T = 6000; om, al, be, delta = 0.05, 0.08, 0.90, 0.0897 s = np.empty(T); r = np.empty(T); e = np.empty(T); s[0] = om / (1 - al - be)98 for t in range(T):99 if t > 0:100 s[t] = om + al * e[t - 1] ** 2 + be * s[t - 1]101 e[t] = np.sqrt(s[t]) * rng.standard_normal(); r[t] = 0.0 + delta * s[t] + e[t]102 f = GARCHM11().fit(r, dist="normal")103 assert abs(f.param("delta") - delta) < 0.05 and abs(f.param("beta") - be) < 0.05104 fc = GARCHM11().forecast(f, r, None, 22)105 assert np.all(fc.step_var > 0) and abs(fc.mean - (f.param("mu") + f.param("delta") * fc.step_var[0])) < 1e-12106107108def test_qgarch_recovery_and_positivity():109 from wp46.models.custom import QGARCH11110 rng = np.random.default_rng(31); T = 6000; om, al, ga, be = 0.06, 0.08, -0.05, 0.88111 s = np.empty(T); e = np.empty(T); s[0] = om / (1 - al - be)112 for t in range(T):113 if t > 0:114 s[t] = om + al * e[t - 1] ** 2 + ga * e[t - 1] + be * s[t - 1]115 e[t] = np.sqrt(s[t]) * rng.standard_normal()116 f = QGARCH11().fit(e, dist="normal")117 assert f.converged and f.admissible118 assert abs(f.param("gamma") - ga) < 0.03 and abs(f.param("beta") - be) < 0.05119 assert f.param("gamma") ** 2 < 4 * f.param("alpha") * f.param("omega")120121122def test_gas_t_nests_garch_as_nu_to_infinity():123 """With nu -> inf the GAS-t score is eps^2 - f, so the recursion equals GARCH(1,1) with alpha_G = alpha, beta_G = beta - alpha."""124 from wp46.models.custom import rec_garch, rec_gas_t125 eps = R[:1500] - R[:1500].mean(); s0 = float(np.var(eps))126 om, al, be = 0.05, 0.07, 0.95127 f = rec_gas_t(eps, om, al, be, 1e9, s0)128 g = rec_garch(eps, om, al, be - al, s0)129 # initializations differ (GAS starts at s0, GARCH at omega + (alpha+beta) s0) and the gap decays at rate beta_G130 assert np.max(np.abs(f[-300:] - g[-300:]) / g[-300:]) < 1e-6131