spb/wp46_uqo

1"""Seeded Monte Carlo recovery: parameters estimated on simulated data are within sampling error of the truth."""2import numpy as np3import pytest45from wp46.models.custom import GJR11, NAGARCH11, HEAVY11, MSGARCH2, RealGARCH116from wp46.models.distributions import rvs789def sim_gjr(T, om, al, ga, be, dist, dp, seed):10    rng = np.random.default_rng(seed)11    z = rvs(dist, np.array(dp), T, rng)12    s = np.empty(T); e = np.empty(T); s[0] = om / (1 - al - be - ga / 2)13    for t in range(T):14        if t > 0:15            s[t] = om + al * e[t - 1] ** 2 + ga * e[t - 1] ** 2 * (e[t - 1] < 0) + be * s[t - 1]16        e[t] = np.sqrt(s[t]) * z[t]17    return e181920@pytest.mark.slow21def test_gjr_recovery_median_over_replications():22    truth = dict(om=0.04, al=0.03, ga=0.10, be=0.88)23    est = []24    for s in range(12):25        r = sim_gjr(4000, dist="t", dp=[6.0], seed=100 + s, **truth)26        f = GJR11().fit(r, dist="t")27        est.append(f.params[1:5])28    med = np.median(np.array(est), axis=0)29    assert np.allclose(med, [truth["om"], truth["al"], truth["ga"], truth["be"]], atol=[0.02, 0.02, 0.03, 0.03])303132@pytest.mark.slow33def test_nagarch_recovery():34    rng = np.random.default_rng(7); T = 500035    om, al, th, be = 0.03, 0.05, 0.8, 0.9036    z = rng.standard_normal(T); s = np.empty(T); e = np.empty(T); s[0] = om / (1 - al * (1 + th ** 2) - be)37    for t in range(T):38        if t > 0:39            s[t] = om + al * (e[t - 1] - th * np.sqrt(s[t - 1])) ** 2 + be * s[t - 1]40        e[t] = np.sqrt(s[t]) * z[t]41    f = NAGARCH11().fit(e, dist="normal")42    assert abs(f.params[3] - th) < 0.3 and abs(f.params[4] - be) < 0.05   # params = [mu, omega, alpha, theta, beta]434445@pytest.mark.slow46def test_realgarch_recovery():47    rng = np.random.default_rng(11); T = 400048    om, be, ga, xi, phi, t1, t2, su = 0.06, 0.55, 0.40, -0.2, 1.0, -0.07, 0.07, 0.449    ls = np.empty(T); r = np.empty(T); lx = np.empty(T); ls[0] = (om + ga * xi) / (1 - be - ga * phi)50    for t in range(T):51        if t > 0:52            ls[t] = om + be * ls[t - 1] + ga * lx[t - 1]53        z = rng.standard_normal(); r[t] = np.exp(0.5 * ls[t]) * z54        lx[t] = xi + phi * ls[t] + t1 * z + t2 * (z * z - 1) + su * rng.standard_normal()55    f = RealGARCH11().fit(r, np.exp(lx), dist="normal")56    assert f.converged57    assert abs(f.params[2] - be) < 0.1 and abs(f.params[3] - ga) < 0.1 and abs(f.params[5] - phi) < 0.158    assert abs(np.exp(f.params[8]) - su) < 0.05596061@pytest.mark.slow62def test_heavy_recovery():63    rng = np.random.default_rng(5); T = 400064    om, al, be, omr, alr, ber = 0.05, 0.4, 0.55, 0.05, 0.45, 0.565    h = np.empty(T); m = np.empty(T); rm = np.empty(T); r = np.empty(T); h[0] = 1.0; m[0] = 1.066    for t in range(T):67        if t > 0:68            h[t] = om + al * rm[t - 1] + be * h[t - 1]; m[t] = omr + alr * rm[t - 1] + ber * m[t - 1]69        r[t] = np.sqrt(h[t]) * rng.standard_normal(); rm[t] = m[t] * rng.gamma(4.0, 0.25)70    f = HEAVY11().fit(r, rm, dist="normal")71    assert f.converged and abs(f.params[2] - al) < 0.1 and abs(f.params[3] - be) < 0.1 and abs(f.params[5] - alr) < 0.1727374@pytest.mark.slow75def test_msgarch_recovery_regime_separation():76    rng = np.random.default_rng(9); T = 600077    P = np.array([[0.98, 0.02], [0.03, 0.97]]); S = 0; s1 = s2 = 1.0; e = np.empty(T)78    par = [(0.02, 0.03, 0.90), (0.3, 0.10, 0.80)]79    for t in range(T):80        if t > 0:81            s1 = par[0][0] + par[0][1] * e[t - 1] ** 2 + par[0][2] * s182            s2 = par[1][0] + par[1][1] * e[t - 1] ** 2 + par[1][2] * s283            S = S if rng.uniform() < P[S, S] else 1 - S84        e[t] = np.sqrt(s1 if S == 0 else s2) * rng.standard_normal()85    f = MSGARCH2().fit(e, dist="normal")86    v = f.params87    assert v[1] > 0.9 and v[2] > 0.9             # persistent regimes recovered88    u1 = v[3] / (1 - v[4] - v[5]); u2 = v[6] / (1 - v[7] - v[8])89    assert u2 / u1 > 3                             # high-variance regime clearly separated90