Repository Public Research preview
1"""Custom maximum-likelihood volatility models (numba recursions + scipy optimizer), with exact admissibility domains2through smooth reparameterizations. Every model documents its equations, constraints and forecast recursion here and in3model_cards/. Returns r_t are in PERCENT; eps_t = r_t - mu; sigma2 in percent^2.45Models (model_id -> class): garch11 (cross-check only), igarch11, gjr11_custom (cross-check), nagarch11, qgarch11,6cgarch11, garchm11, gas_t11 (Beta-t-GARCH / GAS with inverse-information scaling), msgarch2 (Haas–Mittnik–Paolella7two-regime GARCH), realgarch11 (Hansen–Huang–Shek log-linear), heavy11 (Shephard–Sheppard), garchx_rv11 (GARCH(1,1) + lagged8realized measure, RM dynamics from the HEAVY-RM equation), har_rv (Corsi, direct projections per horizon, OLS).9"""10from __future__ import annotations1112import math13import time14import warnings1516import numpy as np17from numba import njit18from scipy import optimize1920from .base import FitResult, Forecast21from .distributions import DIST_BOUNDS, DIST_NPARAMS, DIST_START, loglik_z, neg_half_moment, rvs2223_BIG = 1e10242526# =============================================================================== transforms27def _sig(x): # R -> (0,1)28 return 1.0 / (1.0 + np.exp(-x))293031def _logit(p):32 p = min(max(p, 1e-9), 1 - 1e-9)33 return math.log(p / (1 - p))343536def _dist_theta_to_params(dist: str, th: np.ndarray) -> np.ndarray:37 if dist == "normal":38 return np.empty(0)39 if dist == "t":40 return np.array([2.05 + math.exp(th[0])])41 if dist == "ged":42 return np.array([0.3 + math.exp(th[0])])43 return np.array([2.05 + math.exp(th[0]), 0.99 * math.tanh(th[1])])444546def _dist_params_to_theta(dist: str, dp: np.ndarray) -> np.ndarray:47 if dist == "normal":48 return np.empty(0)49 if dist == "t":50 return np.array([math.log(max(dp[0] - 2.05, 1e-3))])51 if dist == "ged":52 return np.array([math.log(max(dp[0] - 0.3, 1e-3))])53 return np.array([math.log(max(dp[0] - 2.05, 1e-3)), math.atanh(min(max(dp[1] / 0.99, -0.999), 0.999))])545556# =============================================================================== numba recursions57@njit(cache=True)58def rec_garch(eps, omega, alpha, beta, s0):59 T = eps.shape[0]60 s = np.empty(T)61 s[0] = omega + (alpha + beta) * s062 for t in range(1, T):63 s[t] = omega + alpha * eps[t - 1] * eps[t - 1] + beta * s[t - 1]64 return s656667@njit(cache=True)68def rec_gjr(eps, omega, alpha, gamma, beta, s0, kappa):69 T = eps.shape[0]70 s = np.empty(T)71 s[0] = omega + (alpha + gamma * kappa + beta) * s072 for t in range(1, T):73 e2 = eps[t - 1] * eps[t - 1]74 s[t] = omega + alpha * e2 + gamma * e2 * (eps[t - 1] < 0.0) + beta * s[t - 1]75 return s767778@njit(cache=True)79def rec_nagarch(eps, omega, alpha, theta, beta, s0):80 T = eps.shape[0]81 s = np.empty(T)82 s[0] = omega + (alpha * (1.0 + theta * theta) + beta) * s083 for t in range(1, T):84 d = eps[t - 1] - theta * math.sqrt(s[t - 1])85 s[t] = omega + alpha * d * d + beta * s[t - 1]86 return s878889@njit(cache=True)90def rec_qgarch(eps, omega, alpha, gamma, beta, s0):91 T = eps.shape[0]92 s = np.empty(T)93 s[0] = omega + (alpha + beta) * s094 for t in range(1, T):95 s[t] = omega + alpha * eps[t - 1] * eps[t - 1] + gamma * eps[t - 1] + beta * s[t - 1]96 if s[t] < 1e-12:97 s[t] = 1e-1298 return s99100101@njit(cache=True)102def rec_cgarch(eps, omega, rho, phi, alpha, beta, s0):103 T = eps.shape[0]104 s = np.empty(T)105 q = np.empty(T)106 q[0] = s0107 s[0] = s0108 for t in range(1, T):109 e2 = eps[t - 1] * eps[t - 1]110 q[t] = omega + rho * q[t - 1] + phi * (e2 - s[t - 1])111 s[t] = q[t] + alpha * (e2 - q[t - 1]) + beta * (s[t - 1] - q[t - 1])112 if s[t] < 1e-12:113 s[t] = 1e-12114 if q[t] < 1e-12:115 q[t] = 1e-12116 return s, q117118119@njit(cache=True)120def rec_gas_t(eps, omega, alpha, beta, nu, s0):121 # Student-t GAS(1,1) for the variance with inverse-information scaling (Creal–Koopman–Lucas 2013; Harvey 2013)122 T = eps.shape[0]123 f = np.empty(T)124 f[0] = s0125 for t in range(1, T):126 e2 = eps[t - 1] * eps[t - 1]127 if nu > 1e6:128 sc = e2 - f[t - 1]129 else:130 sc = (nu + 3.0) / nu * ((nu + 1.0) * f[t - 1] * e2 / ((nu - 2.0) * f[t - 1] + e2) - f[t - 1])131 f[t] = omega + alpha * sc + beta * f[t - 1]132 if f[t] < 1e-12:133 f[t] = 1e-12134 return f135136137@njit(cache=True)138def rec_garchm(r, mu, delta, omega, alpha, beta, s0):139 T = r.shape[0]140 s = np.empty(T)141 eps = np.empty(T)142 s[0] = omega + (alpha + beta) * s0143 eps[0] = r[0] - mu - delta * s[0]144 for t in range(1, T):145 s[t] = omega + alpha * eps[t - 1] * eps[t - 1] + beta * s[t - 1]146 eps[t] = r[t] - mu - delta * s[t]147 return s, eps148149150@njit(cache=True)151def rec_realgarch(eps, logx, omega, beta, gamma, xi, phi, tau1, tau2, s0):152 # log sigma2_t = omega + beta log sigma2_{t-1} + gamma log x_{t-1}; log x_t = xi + phi log sigma2_t + tau1 z_t + tau2 (z_t^2-1) + u_t153 T = eps.shape[0]154 ls = np.empty(T)155 u = np.empty(T)156 ls[0] = math.log(s0)157 for t in range(T):158 if t > 0:159 ls[t] = omega + beta * ls[t - 1] + gamma * logx[t - 1]160 z = eps[t] / math.exp(0.5 * ls[t])161 u[t] = logx[t] - xi - phi * ls[t] - tau1 * z - tau2 * (z * z - 1.0)162 return ls, u163164165@njit(cache=True)166def rec_heavy_r(rm, omega, alpha, beta, s0):167 T = rm.shape[0]168 h = np.empty(T)169 h[0] = s0170 for t in range(1, T):171 h[t] = omega + alpha * rm[t - 1] + beta * h[t - 1]172 return h173174175@njit(cache=True)176def rec_garchx(eps, rm, omega, alpha, gamma, beta, s0):177 T = eps.shape[0]178 s = np.empty(T)179 s[0] = s0180 for t in range(1, T):181 s[t] = omega + alpha * eps[t - 1] * eps[t - 1] + gamma * rm[t - 1] + beta * s[t - 1]182 return s183184185@njit(cache=True)186def hamilton_msgarch(eps, p11, p22, om1, al1, be1, om2, al2, be2, s0, llz_fn_code, nu, c):187 """Haas–Mittnik–Paolella (2004) two-regime GARCH(1,1) with Markov switching; Hamilton filter.188 Returns: loglik contributions (T,), filtered probabilities of regime 1 (T,), regime variances (T,2).189 llz_fn_code: 0 = normal, 1 = Student-t (nu, c = log normalising constant)."""190 T = eps.shape[0]191 s1 = np.empty(T)192 s2 = np.empty(T)193 ll = np.empty(T)194 pf = np.empty(T) # filtered P(S_t = 1 | F_t)195 # stationary distribution as prior196 pi1 = (1.0 - p22) / (2.0 - p11 - p22)197 s1[0] = s0198 s2[0] = s0199 pred1 = pi1200 for t in range(T):201 if t > 0:202 e2 = eps[t - 1] * eps[t - 1]203 s1[t] = om1 + al1 * e2 + be1 * s1[t - 1]204 s2[t] = om2 + al2 * e2 + be2 * s2[t - 1]205 pred1 = p11 * pf[t - 1] + (1.0 - p22) * (1.0 - pf[t - 1])206 z1 = eps[t] / math.sqrt(s1[t])207 z2 = eps[t] / math.sqrt(s2[t])208 if llz_fn_code == 0:209 l1 = -0.5 * math.log(2.0 * math.pi) - 0.5 * z1 * z1 - 0.5 * math.log(s1[t])210 l2 = -0.5 * math.log(2.0 * math.pi) - 0.5 * z2 * z2 - 0.5 * math.log(s2[t])211 else:212 l1 = c - 0.5 * (nu + 1.0) * math.log1p(z1 * z1 / (nu - 2.0)) - 0.5 * math.log(s1[t])213 l2 = c - 0.5 * (nu + 1.0) * math.log1p(z2 * z2 / (nu - 2.0)) - 0.5 * math.log(s2[t])214 m = max(l1, l2)215 d1 = pred1 * math.exp(l1 - m)216 d2 = (1.0 - pred1) * math.exp(l2 - m)217 den = d1 + d2218 ll[t] = m + math.log(den)219 pf[t] = d1 / den220 out = np.empty((T, 2))221 out[:, 0] = s1222 out[:, 1] = s2223 return ll, pf, out224225226# =============================================================================== generic estimator227class CustomModel:228 model_id = "base"229 family = ""230 vol_names: list[str] = []231 supports = ("normal", "t", "ged", "skewt")232 needs_rm = False233 analytic_multistep = True234235 # --- to be implemented236 def theta_to_vol(self, th: np.ndarray) -> np.ndarray: ...237 def vol_to_theta(self, v: np.ndarray) -> np.ndarray: ...238 def start_vol(self, eps: np.ndarray, s0: float) -> list[np.ndarray]: ...239 def sigma2(self, v: np.ndarray, eps: np.ndarray, x: np.ndarray | None, s0: float, kappa: float) -> np.ndarray: ...240 def step_forecast(self, v: np.ndarray, eps: np.ndarray, s2: np.ndarray, x, h: int, kappa: float, dp, dist, rng) -> np.ndarray: ...241 def admissible(self, v: np.ndarray, kappa: float) -> bool:242 return True243244 # --- shared245 @property246 def param_names(self):247 return ["mu"] + list(self.vol_names)248249 @staticmethod250 def backcast(eps: np.ndarray) -> float:251 """arch-style exponentially weighted backcast of the initial variance (75 lags, decay 0.94)."""252 n = min(75, len(eps))253 w = 0.94 ** np.arange(n)254 w /= w.sum()255 return float(np.sum(w * eps[:n] ** 2))256257 def _negll(self, theta: np.ndarray, r: np.ndarray, x, dist: str, s0: float, kappa_cache: dict) -> float:258 kd = DIST_NPARAMS[dist]259 mu = theta[0]260 v = self.theta_to_vol(theta[1:len(theta) - kd])261 dp = _dist_theta_to_params(dist, theta[len(theta) - kd:]) if kd else np.empty(0)262 # kappa = E[z^2 1{z<0}] only enters the first-observation initialization of the GJR recursion (and the263 # admissibility/forecast rules); it is held at 0.5 inside the optimizer and computed exactly once after the fit264 # (A-04: computing it at every likelihood evaluation made skew-t fits ~100x slower).265 kappa = 0.5266 eps = r - mu267 try:268 s2 = self.sigma2(v, eps, x, s0, kappa)269 except (FloatingPointError, ValueError, ZeroDivisionError):270 return _BIG271 if not np.all(np.isfinite(s2)) or np.any(s2 <= 0):272 return _BIG273 z = eps / np.sqrt(s2)274 ll = loglik_z(dist, z, dp) - 0.5 * np.log(s2)275 ll = ll + self.extra_loglik(v, eps, s2, x)276 val = -float(np.sum(ll))277 return val if np.isfinite(val) else _BIG278279 def extra_loglik(self, v, eps, s2, x):280 return 0.0281282 def fit(self, r: np.ndarray, x: np.ndarray | None = None, dist: str = "normal", max_restarts: int = 3,283 maxiter: int = 2000) -> FitResult:284 t0 = time.perf_counter()285 r = np.asarray(r, float)286 if dist not in self.supports:287 raise ValueError(f"{self.model_id} does not support {dist}")288 mu0 = float(np.mean(r))289 eps0 = r - mu0290 s0 = self.backcast(eps0)291 kd = DIST_NPARAMS[dist]292 starts = []293 for v in self.start_vol(eps0, s0):294 starts.append(np.concatenate([[mu0], self.vol_to_theta(v), _dist_params_to_theta(dist, np.array(DIST_START[dist])) if kd else []]))295 best = None296 n_iter = 0297 n_rest = 0298 cache: dict = {}299 with warnings.catch_warnings():300 warnings.simplefilter("ignore")301 for i, th0 in enumerate(starts[: max_restarts + 1]):302 res = optimize.minimize(self._negll, th0, args=(r, x, dist, s0, cache), method="L-BFGS-B",303 options={"maxiter": maxiter, "ftol": 1e-10, "gtol": 1e-6})304 n_iter += int(res.nit)305 if best is None or res.fun < best.fun - 1e-8:306 best = res307 n_rest = i308 if res.success and i >= 1 and best.fun < _BIG:309 break310 # polish with Nelder–Mead if L-BFGS-B did not converge or stalled immediately (badly scaled gradients),311 # then a final L-BFGS-B pass from the Nelder–Mead optimum312 if best is not None and (not best.success or best.fun >= _BIG or n_iter <= 3 * len(starts[:max_restarts + 1])):313 res = optimize.minimize(self._negll, best.x, args=(r, x, dist, s0, cache), method="Nelder-Mead",314 options={"maxiter": 4000, "xatol": 1e-6, "fatol": 1e-8})315 n_iter += int(res.nit)316 if res.fun < best.fun:317 best = res318 res = optimize.minimize(self._negll, best.x, args=(r, x, dist, s0, cache), method="L-BFGS-B",319 options={"maxiter": maxiter, "ftol": 1e-10, "gtol": 1e-6})320 n_iter += int(res.nit)321 if res.fun < best.fun:322 best = res323 # A-05: likelihood floor against GARCH(1,1) with the same law (models at least as flexible as GARCH(1,1))324 degenerate = False325 if (self.model_id not in ("garch11_custom", "igarch11", "gas_t11") and not self.needs_rm # joint likelihoods not comparable326 and best is not None and best.fun < _BIG):327 try:328 ref = GARCH11().fit(r, None, dist=dist, max_restarts=1).loglik329 if np.isfinite(ref) and -best.fun < ref - 50.0:330 for i in range(2): # bounded retry from a perturbed optimum331 th0 = best.x * (0.8 + 0.4 * np.random.default_rng(i).uniform(size=len(best.x)))332 with warnings.catch_warnings():333 warnings.simplefilter("ignore")334 res = optimize.minimize(self._negll, th0, args=(r, x, dist, s0, cache), method="L-BFGS-B", options={"maxiter": maxiter, "ftol": 1e-10, "gtol": 1e-6})335 n_iter += int(res.nit); n_rest += 1336 if res.fun < best.fun:337 best = res338 if -best.fun >= ref - 50.0:339 break340 degenerate = -best.fun < ref - 50.0341 except Exception: # noqa: BLE001342 pass343 th = best.x344 mu = float(th[0])345 v = self.theta_to_vol(th[1:len(th) - kd])346 dp = _dist_theta_to_params(dist, th[len(th) - kd:]) if kd else np.empty(0)347 kappa = 0.5 if dist != "skewt" else neg_half_moment(dist, dp)348 params = np.concatenate([[mu], v, dp])349 names = self.param_names + (["nu"] if dist in ("t", "ged") else ["eta", "lambda"] if dist == "skewt" else [])350 adm = bool(self.admissible(v, kappa)) and best.fun < _BIG351 conv = bool(best.success) and best.fun < _BIG and not degenerate352 note = "" if adm else "inadmissible at optimum"353 if degenerate:354 note = (note + ";" if note else "") + "degenerate: loglik below GARCH(1,1) reference - 50"355 return FitResult(self.model_id, dist, names, params, -float(best.fun), conv,356 len(r), n_iter, n_rest, time.perf_counter() - t0, adm, note, {"backcast": s0, "kappa": kappa})357358 # --- filtering & forecasting with fixed parameters359 def _unpack(self, fit: FitResult):360 kd = DIST_NPARAMS[fit.dist]361 mu = float(fit.params[0])362 v = fit.params[1:len(fit.params) - kd] if kd else fit.params[1:]363 dp = fit.params[len(fit.params) - kd:] if kd else np.empty(0)364 return mu, np.asarray(v, float), np.asarray(dp, float)365366 def filter(self, fit: FitResult, r: np.ndarray, x: np.ndarray | None = None) -> np.ndarray:367 mu, v, dp = self._unpack(fit)368 eps = np.asarray(r, float) - mu369 return self.sigma2(v, eps, x, fit.aux["backcast"], fit.aux["kappa"])370371 def forecast(self, fit: FitResult, r: np.ndarray, x: np.ndarray | None, h: int, seed: int = 0, n_sim: int = 2000) -> Forecast:372 mu, v, dp = self._unpack(fit)373 eps = np.asarray(r, float) - mu374 s2 = self.sigma2(v, eps, x, fit.aux["backcast"], fit.aux["kappa"])375 rng = np.random.default_rng(seed)376 sv = self.step_forecast(v, eps, s2, x, h, fit.aux["kappa"], dp, fit.dist, rng)377 return Forecast(self.model_id, fit.dist, np.asarray(sv, float), mu, dp,378 "analytic" if self.analytic_multistep else "simulation", float(s2[-1]))379380381# =============================================================================== concrete models382class GARCH11(CustomModel):383 """sigma2_t = omega + alpha eps2_{t-1} + beta sigma2_{t-1}; omega>0, alpha,beta>=0, alpha+beta<1 (reparameterized384 as persistence s=alpha+beta in (0,1) and share w=alpha/s in (0,1))."""385 model_id = "garch11_custom"386 family = "classical"387 vol_names = ["omega", "alpha", "beta"]388389 def theta_to_vol(self, th):390 s, w = _sig(th[1]), _sig(th[2])391 return np.array([math.exp(th[0]), s * w, s * (1 - w)])392393 def vol_to_theta(self, v):394 s = v[1] + v[2]395 return np.array([math.log(v[0]), _logit(s), _logit(v[1] / s if s > 0 else 0.5)])396397 def start_vol(self, eps, s0):398 var = float(np.var(eps))399 return [np.array([var * 0.05, 0.05, 0.90]), np.array([var * 0.02, 0.10, 0.88]), np.array([var * 0.2, 0.15, 0.65])]400401 def sigma2(self, v, eps, x, s0, kappa):402 return rec_garch(eps, v[0], v[1], v[2], s0)403404 def admissible(self, v, kappa):405 return v[0] > 0 and v[1] >= 0 and v[2] >= 0 and v[1] + v[2] < 1406407 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):408 out = np.empty(h)409 out[0] = v[0] + v[1] * eps[-1] ** 2 + v[2] * s2[-1]410 for k in range(1, h):411 out[k] = v[0] + (v[1] + v[2]) * out[k - 1]412 return out413414415class IGARCH11(CustomModel):416 """Integrated GARCH (Engle & Bollerslev 1986): beta = 1 - alpha, omega > 0. Not covariance stationary; strictly417 stationary (Nelson 1990); forecasts grow linearly in h: sigma2_{t+k} = omega + sigma2_{t+k-1}."""418 model_id = "igarch11"419 family = "classical"420 vol_names = ["omega", "alpha"]421422 def theta_to_vol(self, th):423 return np.array([math.exp(th[0]), _sig(th[1])])424425 def vol_to_theta(self, v):426 return np.array([math.log(v[0]), _logit(v[1])])427428 def start_vol(self, eps, s0):429 var = float(np.var(eps))430 return [np.array([var * 0.01, 0.06]), np.array([var * 0.003, 0.03]), np.array([var * 0.03, 0.12])]431432 def sigma2(self, v, eps, x, s0, kappa):433 return rec_garch(eps, v[0], v[1], 1.0 - v[1], s0)434435 def admissible(self, v, kappa):436 return v[0] > 0 and 0 < v[1] < 1437438 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):439 out = np.empty(h)440 out[0] = v[0] + v[1] * eps[-1] ** 2 + (1 - v[1]) * s2[-1]441 for k in range(1, h):442 out[k] = v[0] + out[k - 1]443 return out444445446class GJR11(CustomModel):447 """sigma2_t = omega + alpha eps2 + gamma eps2 1{eps<0} + beta sigma2 (Glosten, Jagannathan & Runkle 1993);448 stationarity alpha + beta + gamma*kappa < 1 with kappa = E[z^2 1{z<0}] (0.5 if symmetric); gamma >= -alpha."""449 model_id = "gjr11_custom"450 family = "asymmetric"451 vol_names = ["omega", "alpha", "gamma", "beta"]452453 def theta_to_vol(self, th):454 s = _sig(th[1])455 w = np.exp(th[2:5])456 w = w / w.sum() # shares of persistence among alpha, gamma*kappa, beta457 return np.array([math.exp(th[0]), s * w[0], s * w[1] / self._kappa, s * w[2]])458459 def vol_to_theta(self, v):460 parts = np.array([max(v[1], 1e-6), max(v[2] * self._kappa, 1e-6), max(v[3], 1e-6)])461 s = parts.sum()462 return np.concatenate([[math.log(v[0]), _logit(s)], np.log(parts / s)])463464 _kappa = 0.5465466 def start_vol(self, eps, s0):467 var = float(np.var(eps))468 return [np.array([var * 0.05, 0.03, 0.08, 0.88]), np.array([var * 0.02, 0.01, 0.12, 0.90]), np.array([var * 0.1, 0.08, 0.05, 0.80])]469470 def sigma2(self, v, eps, x, s0, kappa):471 self._kappa = kappa472 return rec_gjr(eps, v[0], v[1], v[2], v[3], s0, kappa)473474 def admissible(self, v, kappa):475 return v[0] > 0 and v[1] >= 0 and v[3] >= 0 and v[1] + v[2] >= 0 and v[1] + v[3] + v[2] * kappa < 1476477 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):478 out = np.empty(h)479 e = eps[-1]480 out[0] = v[0] + v[1] * e * e + v[2] * e * e * (e < 0) + v[3] * s2[-1]481 for k in range(1, h):482 out[k] = v[0] + (v[1] + v[2] * kappa + v[3]) * out[k - 1]483 return out484485486class NAGARCH11(CustomModel):487 """Nonlinear asymmetric GARCH (Engle & Ng 1993): sigma2_t = omega + alpha (eps_{t-1} - theta sigma_{t-1})^2 + beta sigma2_{t-1}.488 Persistence alpha(1+theta^2) + beta < 1 (E z = 0, Var z = 1); forecasts sigma2_{t+k} = omega + [alpha(1+theta^2)+beta] sigma2_{t+k-1}."""489 model_id = "nagarch11"490 family = "asymmetric"491 vol_names = ["omega", "alpha", "theta", "beta"]492493 def theta_to_vol(self, th):494 s, w, theta = _sig(th[1]), _sig(th[2]), th[3]495 a_eff = s * w496 return np.array([math.exp(th[0]), a_eff / (1 + theta * theta), theta, s * (1 - w)])497498 def vol_to_theta(self, v):499 a_eff = v[1] * (1 + v[2] ** 2)500 s = a_eff + v[3]501 return np.array([math.log(v[0]), _logit(s), _logit(a_eff / s), v[2]])502503 def start_vol(self, eps, s0):504 var = float(np.var(eps))505 return [np.array([var * 0.05, 0.05, 0.5, 0.88]), np.array([var * 0.02, 0.08, 1.0, 0.85]), np.array([var * 0.1, 0.1, 0.0, 0.8])]506507 def sigma2(self, v, eps, x, s0, kappa):508 return rec_nagarch(eps, v[0], v[1], v[2], v[3], s0)509510 def admissible(self, v, kappa):511 return v[0] > 0 and v[1] >= 0 and v[3] >= 0 and v[1] * (1 + v[2] ** 2) + v[3] < 1512513 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):514 out = np.empty(h)515 d = eps[-1] - v[2] * math.sqrt(s2[-1])516 out[0] = v[0] + v[1] * d * d + v[3] * s2[-1]517 pers = v[1] * (1 + v[2] ** 2) + v[3]518 for k in range(1, h):519 out[k] = v[0] + pers * out[k - 1]520 return out521522523class QGARCH11(CustomModel):524 """Quadratic GARCH (Sentana 1995): sigma2_t = omega + alpha eps2_{t-1} + gamma eps_{t-1} + beta sigma2_{t-1}.525 Positivity: omega - gamma^2/(4 alpha) > 0; stationarity alpha + beta < 1; multi-step E[eps]=0 -> GARCH recursion."""526 model_id = "qgarch11"527 family = "asymmetric"528 vol_names = ["omega", "alpha", "gamma", "beta"]529530 def theta_to_vol(self, th):531 s, w = _sig(th[1]), _sig(th[2])532 alpha, beta = s * w, s * (1 - w)533 # gamma^2 < 4 alpha omega <=> gamma = 2 sqrt(alpha omega) * tanh(th[3])534 omega = math.exp(th[0])535 gamma = 2.0 * math.sqrt(alpha * omega) * math.tanh(th[3])536 return np.array([omega, alpha, gamma, beta])537538 def vol_to_theta(self, v):539 s = v[1] + v[3]540 g = v[2] / (2.0 * math.sqrt(max(v[1] * v[0], 1e-12)))541 return np.array([math.log(v[0]), _logit(s), _logit(v[1] / s), math.atanh(min(max(g, -0.99), 0.99))])542543 def start_vol(self, eps, s0):544 var = float(np.var(eps))545 return [np.array([var * 0.05, 0.06, -0.3 * math.sqrt(var * 0.05 * 0.06), 0.90]), np.array([var * 0.03, 0.10, 0.0, 0.87]), np.array([var * 0.1, 0.1, -0.5 * math.sqrt(var * 0.1 * 0.1), 0.8])]546547 def sigma2(self, v, eps, x, s0, kappa):548 return rec_qgarch(eps, v[0], v[1], v[2], v[3], s0)549550 def admissible(self, v, kappa):551 return v[0] > 0 and v[1] > 0 and v[3] >= 0 and v[1] + v[3] < 1 and v[2] ** 2 < 4 * v[1] * v[0]552553 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):554 out = np.empty(h)555 out[0] = max(v[0] + v[1] * eps[-1] ** 2 + v[2] * eps[-1] + v[3] * s2[-1], 1e-12)556 for k in range(1, h):557 out[k] = v[0] + (v[1] + v[3]) * out[k - 1]558 return out559560561class CGARCH11(CustomModel):562 """Component GARCH (Engle & Lee 1999): sigma2_t = q_t + alpha(eps2_{t-1} - q_{t-1}) + beta(sigma2_{t-1} - q_{t-1}),563 q_t = omega + rho q_{t-1} + phi(eps2_{t-1} - sigma2_{t-1}). Constraints 0 < alpha+beta < rho < 1, phi > 0,564 beta > phi, alpha > phi (Engle & Lee's sufficient positivity conditions); forecasts: q_{t+k} = omega + rho q_{t+k-1},565 (sigma2 - q)_{t+k} = (alpha+beta)(sigma2 - q)_{t+k-1}."""566 model_id = "cgarch11"567 family = "component"568 vol_names = ["omega", "rho", "phi", "alpha", "beta"]569570 def theta_to_vol(self, th):571 rho = _sig(th[1])572 s = rho * _sig(th[2]) # alpha + beta < rho573 w = _sig(th[3])574 alpha, beta = s * w, s * (1 - w)575 phi = min(alpha, beta) * _sig(th[4]) # phi < min(alpha, beta)576 return np.array([math.exp(th[0]), rho, phi, alpha, beta])577578 def vol_to_theta(self, v):579 s = v[3] + v[4]580 return np.array([math.log(v[0]), _logit(v[1]), _logit(s / v[1]), _logit(v[3] / s), _logit(v[2] / max(min(v[3], v[4]), 1e-9))])581582 def start_vol(self, eps, s0):583 var = float(np.var(eps))584 return [np.array([var * 0.01, 0.99, 0.02, 0.05, 0.85]), np.array([var * 0.005, 0.995, 0.03, 0.08, 0.80]), np.array([var * 0.02, 0.98, 0.01, 0.04, 0.90])]585586 def sigma2(self, v, eps, x, s0, kappa):587 s, q = rec_cgarch(eps, v[0], v[1], v[2], v[3], v[4], s0)588 self._q_last = q[-1]589 return s590591 def admissible(self, v, kappa):592 return v[0] > 0 and 0 < v[3] + v[4] < v[1] < 1 and 0 < v[2] < min(v[3], v[4])593594 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):595 s, q = rec_cgarch(eps, v[0], v[1], v[2], v[3], v[4], s2[0] if len(s2) else 1.0)596 e2 = eps[-1] ** 2597 q1 = v[0] + v[1] * q[-1] + v[2] * (e2 - s[-1])598 s1 = q1 + v[3] * (e2 - q[-1]) + v[4] * (s[-1] - q[-1])599 out = np.empty(h)600 out[0] = max(s1, 1e-12)601 qk, dk = q1, s1 - q1602 for k in range(1, h):603 qk = v[0] + v[1] * qk604 dk = (v[3] + v[4]) * dk605 out[k] = max(qk + dk, 1e-12)606 return out607608609class GARCHM11(CustomModel):610 """GARCH-in-mean (Engle, Lilien & Robins 1987): r_t = mu + delta sigma2_t + eps_t with GARCH(1,1) variance.611 Variance forecasts as GARCH(1,1); the mean forecast is mu + delta sigma2_{t+1}."""612 model_id = "garchm11"613 family = "classical"614 vol_names = ["delta", "omega", "alpha", "beta"]615616 def theta_to_vol(self, th):617 s, w = _sig(th[2]), _sig(th[3])618 return np.array([th[0], math.exp(th[1]), s * w, s * (1 - w)])619620 def vol_to_theta(self, v):621 s = v[2] + v[3]622 return np.array([v[0], math.log(v[1]), _logit(s), _logit(v[2] / s)])623624 def start_vol(self, eps, s0):625 var = float(np.var(eps))626 return [np.array([0.0, var * 0.05, 0.05, 0.90]), np.array([0.02, var * 0.02, 0.1, 0.88]), np.array([-0.02, var * 0.2, 0.15, 0.65])]627628 def _negll(self, theta, r, x, dist, s0, cache):629 kd = DIST_NPARAMS[dist]630 mu = theta[0]631 v = self.theta_to_vol(theta[1:len(theta) - kd])632 dp = _dist_theta_to_params(dist, theta[len(theta) - kd:]) if kd else np.empty(0)633 s2, eps = rec_garchm(r, mu, v[0], v[1], v[2], v[3], s0)634 if not np.all(np.isfinite(s2)) or np.any(s2 <= 0):635 return _BIG636 z = eps / np.sqrt(s2)637 val = -float(np.sum(loglik_z(dist, z, dp) - 0.5 * np.log(s2)))638 return val if np.isfinite(val) else _BIG639640 def sigma2(self, v, eps, x, s0, kappa): # eps here is r - mu; recompute with the in-mean term641 raise NotImplementedError("use filter()")642643 def filter(self, fit, r, x=None):644 mu, v, dp = self._unpack(fit)645 s2, eps = rec_garchm(np.asarray(r, float), mu, v[0], v[1], v[2], v[3], fit.aux["backcast"])646 return s2647648 def admissible(self, v, kappa):649 return v[1] > 0 and v[2] >= 0 and v[3] >= 0 and v[2] + v[3] < 1650651 def forecast(self, fit, r, x, h, seed=0, n_sim=2000):652 mu, v, dp = self._unpack(fit)653 s2, eps = rec_garchm(np.asarray(r, float), mu, v[0], v[1], v[2], v[3], fit.aux["backcast"])654 out = np.empty(h)655 out[0] = v[1] + v[2] * eps[-1] ** 2 + v[3] * s2[-1]656 for k in range(1, h):657 out[k] = v[1] + (v[2] + v[3]) * out[k - 1]658 return Forecast(self.model_id, fit.dist, out, mu + v[0] * out[0], dp, "analytic", float(s2[-1]))659660661class GAST11(CustomModel):662 """Score-driven (GAS) Student-t volatility model with inverse-information scaling (Creal, Koopman & Lucas 2013;663 equivalent to Harvey's (2013) Beta-t-GARCH up to parameterization): f_{t} = omega + alpha s_{t-1} + beta f_{t-1},664 s_t = (nu+3)/nu [ (nu+1) f_t eps2_t / ((nu-2) f_t + eps2_t) - f_t ]. For nu -> inf, s_t = eps2_t - f_t (GARCH with665 alpha_G = alpha, beta_G = beta - alpha). Stationarity: 0 < beta < 1; forecasts f_{t+k} = omega + beta f_{t+k-1}."""666 model_id = "gas_t11"667 family = "score-driven"668 vol_names = ["omega", "alpha", "beta"]669 supports = ("t",)670671 def theta_to_vol(self, th):672 return np.array([math.exp(th[0]), math.exp(th[1]), _sig(th[2])])673674 def vol_to_theta(self, v):675 return np.array([math.log(v[0]), math.log(v[1]), _logit(v[2])])676677 def start_vol(self, eps, s0):678 var = float(np.var(eps))679 return [np.array([var * 0.03, 0.05, 0.97]), np.array([var * 0.01, 0.08, 0.99]), np.array([var * 0.1, 0.03, 0.9])]680681 def _negll(self, theta, r, x, dist, s0, cache):682 mu = theta[0]683 v = self.theta_to_vol(theta[1:-1])684 dp = _dist_theta_to_params("t", theta[-1:])685 eps = r - mu686 s2 = rec_gas_t(eps, v[0], v[1], v[2], float(dp[0]), s0)687 if not np.all(np.isfinite(s2)) or np.any(s2 <= 0):688 return _BIG689 z = eps / np.sqrt(s2)690 val = -float(np.sum(loglik_z("t", z, dp) - 0.5 * np.log(s2)))691 return val if np.isfinite(val) else _BIG692693 def sigma2(self, v, eps, x, s0, kappa):694 return rec_gas_t(eps, v[0], v[1], v[2], self._nu, s0)695696 def filter(self, fit, r, x=None):697 mu, v, dp = self._unpack(fit)698 self._nu = float(dp[0])699 return rec_gas_t(np.asarray(r, float) - mu, v[0], v[1], v[2], self._nu, fit.aux["backcast"])700701 def admissible(self, v, kappa):702 return v[0] > 0 and v[1] >= 0 and 0 <= v[2] < 1703704 def forecast(self, fit, r, x, h, seed=0, n_sim=2000):705 mu, v, dp = self._unpack(fit)706 nu = float(dp[0])707 eps = np.asarray(r, float) - mu708 f = rec_gas_t(eps, v[0], v[1], v[2], nu, fit.aux["backcast"])709 e2 = eps[-1] ** 2710 sc = (nu + 3.0) / nu * ((nu + 1.0) * f[-1] * e2 / ((nu - 2.0) * f[-1] + e2) - f[-1])711 out = np.empty(h)712 out[0] = v[0] + v[1] * sc + v[2] * f[-1]713 for k in range(1, h):714 out[k] = v[0] + v[2] * out[k - 1]715 return Forecast(self.model_id, "t", out, mu, dp, "analytic", float(f[-1]))716717718class MSGARCH2(CustomModel):719 """Two-regime Markov-switching GARCH(1,1) of Haas, Mittnik & Paolella (2004): regime-specific variance recursions720 driven by the common eps_t, Markov chain with transition probabilities p11, p22, regimes ordered by unconditional721 variance (identification). Innovations normal or Student-t with common nu. Hamilton filter likelihood; multi-step722 forecasts by the exact recursion E[sigma2_{k,t+h}|F_t] (HMP 2004, eq. for the 'Haas' specification):723 E[eps2_{t+j}] = sum_k P(S_{t+j}=k|F_t) E[sigma2_{k,t+j}], E[sigma2_{k,t+j+1}] = omega_k + alpha_k E[eps2_{t+j}] + beta_k E[sigma2_{k,t+j}]."""724 model_id = "msgarch2"725 family = "regime-switching"726 vol_names = ["p11", "p22", "omega1", "alpha1", "beta1", "omega2", "alpha2", "beta2"]727 supports = ("normal", "t")728729 def theta_to_vol(self, th):730 p11, p22 = _sig(th[0]), _sig(th[1])731 s1, w1 = _sig(th[3]), _sig(th[4])732 s2, w2 = _sig(th[6]), _sig(th[7])733 return np.array([p11, p22, math.exp(th[2]), s1 * w1, s1 * (1 - w1), math.exp(th[5]), s2 * w2, s2 * (1 - w2)])734735 def vol_to_theta(self, v):736 s1, s2 = v[3] + v[4], v[6] + v[7]737 return np.array([_logit(v[0]), _logit(v[1]), math.log(v[2]), _logit(s1), _logit(v[3] / s1), math.log(v[5]), _logit(s2), _logit(v[6] / s2)])738739 def start_vol(self, eps, s0):740 var = float(np.var(eps))741 return [np.array([0.98, 0.97, var * 0.5 * 0.08, 0.03, 0.90, var * 2.0 * 0.15, 0.10, 0.80]),742 np.array([0.95, 0.95, var * 0.3 * 0.1, 0.05, 0.85, var * 3.0 * 0.1, 0.15, 0.75]),743 np.array([0.99, 0.99, var * 0.7 * 0.05, 0.02, 0.93, var * 1.5 * 0.1, 0.08, 0.85])]744745 def _ll_parts(self, theta, r, dist, s0):746 kd = DIST_NPARAMS[dist]747 mu = theta[0]748 v = self.theta_to_vol(theta[1:len(theta) - kd])749 dp = _dist_theta_to_params(dist, theta[len(theta) - kd:]) if kd else np.empty(0)750 eps = r - mu751 if dist == "t":752 from scipy.special import gammaln753 nu = float(dp[0])754 c = float(gammaln((nu + 1) / 2) - gammaln(nu / 2) - 0.5 * math.log(math.pi * (nu - 2)))755 code = 1756 else:757 nu, c, code = 1e9, 0.0, 0758 return hamilton_msgarch(eps, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], s0, code, nu, c), v, dp, eps759760 def _negll(self, theta, r, x, dist, s0, cache):761 try:762 (ll, pf, S), v, dp, eps = self._ll_parts(theta, r, dist, s0)763 except (FloatingPointError, ValueError, ZeroDivisionError):764 return _BIG765 if not np.all(np.isfinite(ll)):766 return _BIG767 # identification: regime 1 = low unconditional variance (penalize violations smoothly)768 u1 = v[2] / max(1 - v[3] - v[4], 1e-6)769 u2 = v[5] / max(1 - v[6] - v[7], 1e-6)770 pen = 0.0 if u1 <= u2 else 50.0 * (math.log(u1) - math.log(u2))771 return -float(np.sum(ll)) + pen772773 def sigma2(self, v, eps, x, s0, kappa):774 raise NotImplementedError("use filter()")775776 def filter(self, fit, r, x=None):777 kd = DIST_NPARAMS[fit.dist]778 th = np.concatenate([[fit.params[0]], self.vol_to_theta(fit.params[1:len(fit.params) - kd] if kd else fit.params[1:]),779 _dist_params_to_theta(fit.dist, fit.params[len(fit.params) - kd:]) if kd else []])780 (ll, pf, S), v, dp, eps = self._ll_parts(th, np.asarray(r, float), fit.dist, fit.aux["backcast"])781 # one-step-ahead predictive mixture variance at each t (given F_{t-1}) is pred1 S1 + (1-pred1) S2782 pred1 = np.empty(len(pf))783 pi1 = (1 - v[1]) / (2 - v[0] - v[1])784 pred1[0] = pi1785 pred1[1:] = v[0] * pf[:-1] + (1 - v[1]) * (1 - pf[:-1])786 return pred1 * S[:, 0] + (1 - pred1) * S[:, 1]787788 def admissible(self, v, kappa):789 return all(v[i] > 0 for i in (2, 5)) and 0 < v[0] < 1 and 0 < v[1] < 1790791 def forecast(self, fit, r, x, h, seed=0, n_sim=2000):792 kd = DIST_NPARAMS[fit.dist]793 th = np.concatenate([[fit.params[0]], self.vol_to_theta(fit.params[1:len(fit.params) - kd] if kd else fit.params[1:]),794 _dist_params_to_theta(fit.dist, fit.params[len(fit.params) - kd:]) if kd else []])795 (ll, pf, S), v, dp, eps = self._ll_parts(th, np.asarray(r, float), fit.dist, fit.aux["backcast"])796 p11, p22 = v[0], v[1]797 P = np.array([[p11, 1 - p11], [1 - p22, p22]]) # rows: from798 prob = np.array([pf[-1], 1 - pf[-1]]) # P(S_t | F_t)799 e2 = eps[-1] ** 2800 s_k = np.array([v[2] + v[3] * e2 + v[4] * S[-1, 0], v[5] + v[6] * e2 + v[7] * S[-1, 1]]) # sigma2_{k,t+1}801 out = np.empty(h)802 for j in range(h):803 prob = prob @ P # P(S_{t+j+1} | F_t)804 out[j] = float(prob @ s_k)805 e2f = out[j]806 s_k = np.array([v[2] + v[3] * e2f + v[4] * s_k[0], v[5] + v[6] * e2f + v[7] * s_k[1]])807 mix_sigma2 = out[0]808 fc = Forecast(self.model_id, fit.dist, out, float(fit.params[0]), dp, "analytic", float(mix_sigma2))809 fc.regime_probs = prob # type: ignore[attr-defined]810 fc.regime_var = s_k # type: ignore[attr-defined]811 return fc812813814class RealGARCH11(CustomModel):815 """Realized GARCH(1,1), log-linear specification (Hansen, Huang & Shek 2012):816 r_t = mu + sigma_t z_t; log sigma2_t = omega + beta log sigma2_{t-1} + gamma log x_{t-1};817 log x_t = xi + phi log sigma2_t + tau1 z_t + tau2 (z_t^2 - 1) + u_t, u_t ~ N(0, sigma_u^2).818 Joint likelihood ell = sum_t [ log f(z_t) - 0.5 log sigma2_t ] + [ -0.5 log(2 pi sigma_u^2) - u_t^2 / (2 sigma_u^2) ].819 x_t is the realized measure observed at the END of day t (percent^2); only x_{t-1} enters sigma2_t (no look-ahead).820 Stationarity of log sigma2: |beta + phi gamma| < 1. Multi-step forecasts by simulation (default 2000 paths)."""821 model_id = "realgarch11"822 family = "realized"823 vol_names = ["omega", "beta", "gamma", "xi", "phi", "tau1", "tau2", "log_sigma_u"]824 supports = ("normal", "t")825 needs_rm = True826 analytic_multistep = False827828 def theta_to_vol(self, th):829 return np.array([th[0], th[1], th[2], th[3], th[4], th[5], th[6], th[7]])830831 def vol_to_theta(self, v):832 return np.asarray(v, float)833834 def start_vol(self, eps, s0):835 lx = float(np.mean(np.log(np.maximum(self._x, 1e-8))))836 return [np.array([0.04 * lx, 0.55, 0.41, -0.1, 1.0, -0.07, 0.07, math.log(0.4)]),837 np.array([0.1 * lx, 0.7, 0.28, 0.0, 0.95, -0.05, 0.05, math.log(0.5)]),838 np.array([0.02 * lx, 0.4, 0.55, -0.2, 1.05, -0.1, 0.1, math.log(0.35)])]839840 def fit(self, r, x=None, dist="normal", max_restarts=3, maxiter=3000):841 self._x = np.asarray(x, float)842 return super().fit(r, x, dist, max_restarts, maxiter)843844 def sigma2(self, v, eps, x, s0, kappa):845 logx = np.log(np.maximum(np.asarray(x, float), 1e-8))846 ls, u = rec_realgarch(eps, logx, v[0], v[1], v[2], v[3], v[4], v[5], v[6], s0)847 self._u = u848 if np.any(ls > 50) or np.any(ls < -50):849 raise FloatingPointError("log variance out of range")850 return np.exp(ls)851852 def extra_loglik(self, v, eps, s2, x):853 su2 = math.exp(2 * v[7])854 return -0.5 * math.log(2 * math.pi * su2) - self._u ** 2 / (2 * su2)855856 def admissible(self, v, kappa):857 return abs(v[1] + v[4] * v[2]) < 1858859 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):860 omega, beta, gamma, xi, phi, tau1, tau2, lsu = v861 su = math.exp(lsu)862 logx_last = math.log(max(float(x[-1]), 1e-8))863 ls1 = omega + beta * math.log(s2[-1]) + gamma * logx_last864 if h == 1:865 return np.array([math.exp(ls1)])866 n = 2000867 ls = np.full(n, ls1)868 out = np.empty(h)869 out[0] = math.exp(ls1)870 for k in range(1, h):871 z = rvs(dist, dp, n, rng)872 u = rng.standard_normal(n) * su873 lx = xi + phi * ls + tau1 * z + tau2 * (z * z - 1) + u874 ls = omega + beta * ls + gamma * lx875 ls = np.clip(ls, -30, 30)876 out[k] = float(np.mean(np.exp(ls)))877 return out878879880class HEAVY11(CustomModel):881 """HEAVY model (Shephard & Sheppard 2010): h_t = omega + alpha RM_{t-1} + beta h_{t-1} (Var[r_t|F_{t-1}]) and882 mu_t = omega_R + alpha_R RM_{t-1} + beta_R mu_{t-1} (E[RM_t|F_{t-1}]). Each equation is estimated separately by883 quasi-maximum likelihood: Gaussian QL for returns, and the multiplicative-error (exponential/Gamma, i.e. QLIKE) QL884 -sum(log mu_t + RM_t/mu_t) for the realized measure. Multi-step: E[h_{t+k}] = omega + alpha E[mu_{t+k-1}] + beta E[h_{t+k-1}],885 E[mu_{t+k}] = omega_R + (alpha_R + beta_R) E[mu_{t+k-1}]. Constraints: omega, omega_R > 0, alpha, alpha_R >= 0,886 0 <= beta < 1, alpha_R + beta_R < 1."""887 model_id = "heavy11"888 family = "realized"889 vol_names = ["omega", "alpha", "beta", "omega_r", "alpha_r", "beta_r"]890 supports = ("normal", "t")891 needs_rm = True892893 def theta_to_vol(self, th):894 s = _sig(th[4]) if len(th) > 4 else 0.0895 w = _sig(th[5]) if len(th) > 5 else 0.0896 return np.array([math.exp(th[0]), math.exp(th[1]), _sig(th[2]), math.exp(th[3]), s * w, s * (1 - w)])897898 def vol_to_theta(self, v):899 s = v[4] + v[5]900 return np.array([math.log(v[0]), math.log(v[1]), _logit(v[2]), math.log(v[3]), _logit(s), _logit(v[4] / s)])901902 def start_vol(self, eps, s0):903 m = float(np.mean(self._x))904 return [np.array([0.05 * m, 0.4, 0.6, 0.05 * m, 0.4, 0.55]), np.array([0.1 * m, 0.3, 0.65, 0.1 * m, 0.3, 0.6]), np.array([0.02 * m, 0.5, 0.5, 0.02 * m, 0.5, 0.48])]905906 def fit(self, r, x=None, dist="normal", max_restarts=3, maxiter=2000):907 self._x = np.asarray(x, float)908 return super().fit(r, x, dist, max_restarts, maxiter)909910 def sigma2(self, v, eps, x, s0, kappa):911 return rec_heavy_r(np.asarray(x, float), v[0], v[1], v[2], s0)912913 def extra_loglik(self, v, eps, s2, x):914 rm = np.asarray(x, float)915 mu = rec_heavy_r(rm, v[3], v[4], v[5], float(np.mean(rm[:min(75, len(rm))])))916 return -(np.log(mu) + rm / mu)917918 def admissible(self, v, kappa):919 return v[0] > 0 and v[1] >= 0 and 0 <= v[2] < 1 and v[3] > 0 and v[4] >= 0 and v[4] + v[5] < 1920921 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):922 rm = np.asarray(x, float)923 mu = rec_heavy_r(rm, v[3], v[4], v[5], float(np.mean(rm[:min(75, len(rm))])))924 out = np.empty(h)925 out[0] = v[0] + v[1] * rm[-1] + v[2] * s2[-1]926 m = v[3] + v[4] * rm[-1] + v[5] * mu[-1] # E[RM_{t+1}]927 for k in range(1, h):928 out[k] = v[0] + v[1] * m + v[2] * out[k - 1]929 m = v[3] + (v[4] + v[5]) * m930 return out931932933class GARCHXRV11(HEAVY11):934 """GARCH-X with a lagged realized measure (Engle 2002): sigma2_t = omega + alpha eps2_{t-1} + gamma RM_{t-1} + beta sigma2_{t-1};935 the realized measure follows the HEAVY-RM equation for multi-step forecasts. Stationarity alpha + beta < 1 (given RM stationary)."""936 model_id = "garchx_rv11"937 family = "realized"938 vol_names = ["omega", "alpha", "gamma", "beta", "omega_r", "alpha_r", "beta_r"]939940 def theta_to_vol(self, th):941 s, w = _sig(th[1]), _sig(th[2])942 sr, wr = _sig(th[5]), _sig(th[6])943 return np.array([math.exp(th[0]), s * w, math.exp(th[3]), s * (1 - w), math.exp(th[4]), sr * wr, sr * (1 - wr)])944945 def vol_to_theta(self, v):946 s, sr = v[1] + v[3], v[5] + v[6]947 return np.array([math.log(v[0]), _logit(s), _logit(v[1] / s), math.log(v[2]), math.log(v[4]), _logit(sr), _logit(v[5] / sr)])948949 def start_vol(self, eps, s0):950 m = float(np.mean(self._x))951 return [np.array([0.05 * m, 0.03, 0.3, 0.6, 0.05 * m, 0.4, 0.55]), np.array([0.1 * m, 0.05, 0.2, 0.7, 0.1 * m, 0.3, 0.6]), np.array([0.02 * m, 0.1, 0.4, 0.45, 0.02 * m, 0.5, 0.48])]952953 def sigma2(self, v, eps, x, s0, kappa):954 return rec_garchx(eps, np.asarray(x, float), v[0], v[1], v[2], v[3], s0)955956 def extra_loglik(self, v, eps, s2, x):957 rm = np.asarray(x, float)958 mu = rec_heavy_r(rm, v[4], v[5], v[6], float(np.mean(rm[:min(75, len(rm))])))959 return -(np.log(mu) + rm / mu)960961 def admissible(self, v, kappa):962 return v[0] > 0 and v[1] >= 0 and v[2] >= 0 and v[3] >= 0 and v[1] + v[3] < 1 and v[4] > 0 and v[5] + v[6] < 1963964 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):965 rm = np.asarray(x, float)966 mu = rec_heavy_r(rm, v[4], v[5], v[6], float(np.mean(rm[:min(75, len(rm))])))967 out = np.empty(h)968 out[0] = v[0] + v[1] * eps[-1] ** 2 + v[2] * rm[-1] + v[3] * s2[-1]969 m = v[4] + v[5] * rm[-1] + v[6] * mu[-1]970 for k in range(1, h):971 out[k] = v[0] + (v[1] + v[3]) * out[k - 1] + v[2] * m972 m = v[4] + (v[5] + v[6]) * m973 return out974975976MODELS: dict[str, type[CustomModel]] = {m.model_id: m for m in977 [GARCH11, IGARCH11, GJR11, NAGARCH11, QGARCH11, CGARCH11, GARCHM11, GAST11, MSGARCH2,978 RealGARCH11, HEAVY11, GARCHXRV11]}979