spb/wp46_uqo

1"""Standardized innovation distributions (zero mean, unit variance) used by the custom likelihoods.23Parameterizations follow the `arch` package (version 8.0.0) so that log-likelihoods are directly comparable:4  normal          no parameter5  t               nu > 2, Bollerslev (1987) standardized Student-t6  ged             nu > 0, Nelson (1991) generalized error distribution, standardized7  skewt           eta > 2, -1 < lambda < 1, Hansen (1994) skewed Student-t (arch's SkewStudent)8All log densities are those of the STANDARDIZED variable z; the density of eps = sigma z adds -log(sigma).9"""10from __future__ import annotations1112import math1314import numpy as np15from numba import njit16from scipy import stats17from scipy.special import gammaln1819DIST_NPARAMS = {"normal": 0, "t": 1, "ged": 1, "skewt": 2}20DIST_BOUNDS = {"normal": [], "t": [(2.05, 200.0)], "ged": [(0.5, 20.0)], "skewt": [(2.05, 200.0), (-0.99, 0.99)]}21DIST_START = {"normal": [], "t": [6.0], "ged": [1.5], "skewt": [6.0, -0.05]}222324@njit(cache=True)25def _ll_normal(z: np.ndarray) -> np.ndarray:26    return -0.5 * math.log(2.0 * math.pi) - 0.5 * z * z272829@njit(cache=True)30def _ll_t(z: np.ndarray, nu: float, c: float) -> np.ndarray:31    # c = lgamma((nu+1)/2) - lgamma(nu/2) - 0.5*log(pi*(nu-2))32    return c - 0.5 * (nu + 1.0) * np.log1p(z * z / (nu - 2.0))333435@njit(cache=True)36def _ll_ged(z: np.ndarray, nu: float, lam: float, c: float) -> np.ndarray:37    # c = log(nu/lam) - (1+1/nu)*log(2) - lgamma(1/nu)38    return c - 0.5 * np.abs(z / lam) ** nu394041@njit(cache=True)42def _ll_skewt(z: np.ndarray, eta: float, lam: float, a: float, b: float, c: float) -> np.ndarray:43    out = np.empty_like(z)44    for i in range(z.shape[0]):45        if z[i] < -a / b:46            y = (b * z[i] + a) / (1.0 - lam)47        else:48            y = (b * z[i] + a) / (1.0 + lam)49        out[i] = c + math.log(b) - 0.5 * (eta + 1.0) * math.log1p(y * y / (eta - 2.0))50    return out515253def dist_constants(dist: str, dp: np.ndarray) -> tuple:54    if dist == "normal":55        return ()56    if dist == "t":57        nu = float(dp[0])58        return (nu, float(gammaln((nu + 1) / 2) - gammaln(nu / 2) - 0.5 * math.log(math.pi * (nu - 2))))59    if dist == "ged":60        nu = float(dp[0])61        lam = math.sqrt(2.0 ** (-2.0 / nu) * math.exp(gammaln(1.0 / nu) - gammaln(3.0 / nu)))62        c = math.log(nu / lam) - (1.0 + 1.0 / nu) * math.log(2.0) - float(gammaln(1.0 / nu))63        return (nu, lam, c)64    if dist == "skewt":65        eta, lam = float(dp[0]), float(dp[1])66        c = float(gammaln((eta + 1) / 2) - gammaln(eta / 2) - 0.5 * math.log(math.pi * (eta - 2)))67        cc = math.exp(c)68        a = 4.0 * lam * cc * (eta - 2.0) / (eta - 1.0)69        b = math.sqrt(1.0 + 3.0 * lam * lam - a * a)70        return (eta, lam, a, b, c)71    raise ValueError(dist)727374def loglik_z(dist: str, z: np.ndarray, dp: np.ndarray) -> np.ndarray:75    k = dist_constants(dist, dp)76    if dist == "normal":77        return _ll_normal(z)78    if dist == "t":79        return _ll_t(z, k[0], k[1])80    if dist == "ged":81        return _ll_ged(z, k[0], k[1], k[2])82    return _ll_skewt(z, k[0], k[1], k[2], k[3], k[4])838485# ---------------------------------------------------------------------------- quantiles / moments (via arch, verified)86def _arch_dist(dist: str):87    from arch.univariate import distribution as D88    return {"normal": D.Normal, "t": D.StudentsT, "ged": D.GeneralizedError, "skewt": D.SkewStudent}[dist]()899091def ppf(dist: str, p: np.ndarray | float, dp: np.ndarray) -> np.ndarray:92    """Quantiles of the standardized innovation (arch 8.0.0 implementation, parameter order as above)."""93    return np.asarray(_arch_dist(dist).ppf(np.atleast_1d(p), list(dp)))949596def expected_shortfall(dist: str, p: float, dp: np.ndarray, n: int = 2000) -> float:97    """ES_p = E[z | z <= q_p] of the standardized innovation, by midpoint quadrature over the lower tail."""98    u = (np.arange(n) + 0.5) / n * p99    return float(np.mean(ppf(dist, u, dp)))100101102def neg_half_moment(dist: str, dp: np.ndarray, n: int = 20000) -> float:103    """kappa = E[z^2 1{z<0}] (0.5 for symmetric distributions); needed for GJR/TARCH multi-step forecasts."""104    if dist in ("normal", "t", "ged"):105        return 0.5106    u = (np.arange(n) + 0.5) / n107    q = ppf(dist, u, dp)108    return float(np.mean(q * q * (q < 0)))109110111def var_z2(dist: str, dp: np.ndarray, n: int = 20000) -> float:112    """Var(z^2) (= 2 for the normal); used by Realized GARCH diagnostics."""113    u = (np.arange(n) + 0.5) / n114    q = ppf(dist, u, dp)115    return float(np.var(q * q))116117118def std_t_rvs(nu: float, size: int, rng: np.random.Generator) -> np.ndarray:119    return rng.standard_t(nu, size) * math.sqrt((nu - 2.0) / nu)120121122def rvs(dist: str, dp: np.ndarray, size: int, rng: np.random.Generator) -> np.ndarray:123    """Random standardized innovations (inverse-cdf for ged/skewt via arch ppf)."""124    if dist == "normal":125        return rng.standard_normal(size)126    if dist == "t":127        return std_t_rvs(float(dp[0]), size, rng)128    u = rng.uniform(1e-7, 1 - 1e-7, size)129    return ppf(dist, u, dp)130131132def cdf(dist: str, z: np.ndarray, dp: np.ndarray) -> np.ndarray:133    return np.asarray(_arch_dist(dist).cdf(np.atleast_1d(z), list(dp)))134