Repository Public Research preview
1"""Post-registration extensions (labelled 'extension' in the registry; scored with the same protocol but reported separately2from the pre-registered set): smooth-transition GARCH, Beta-t-EGARCH, Realized EGARCH, GARCH-MIDAS, HYGARCH, FIAPARCH, FIEGARCH.3All share the custom MLE framework of custom.py; the fractional models share the truncated ARCH(infinity) weight machinery4(truncation 1000 lags, as in arch's FIGARCH), so that HYGARCH (amp = 1) and FIAPARCH (delta = 2, gamma = 0) coincide with5FIGARCH(1,d,1) up to the pre-sample treatment; the tests evaluate our HYGARCH likelihood at arch's FIGARCH optimum (|diff| < 56over 3,000 observations) and check the nesting inequalities. These models have no other external cross-check (exploratory)."""7from __future__ import annotations89import math1011import numpy as np12from numba import njit1314from .custom import CustomModel, GARCH11, _dist_theta_to_params, _logit, _sig, _BIG15from .distributions import DIST_NPARAMS, loglik_z, rvs1617TRUNC = 1000181920# =============================================================================== fractional weights21def frac_diff_coefs(d: float, n: int) -> np.ndarray:22 """pi_k of (1-L)^d = sum_k pi_k L^k, k = 0..n-1 (pi_0 = 1)."""23 pi = np.empty(n); pi[0] = 1.024 for k in range(1, n):25 pi[k] = pi[k - 1] * (k - 1 - d) / k26 return pi272829def figarch_weights(d: float, phi: float, beta: float, amp: float = 1.0, n: int = TRUNC) -> np.ndarray:30 """lambda_k, k = 1..n, of lambda(L) = 1 - (1 - beta L)^{-1} (1 - phi L) [1 + amp((1-L)^d - 1)] (amp = 1: FIGARCH;31 amp != 1: HYGARCH of Davidson 2004). Returns the array lambda_1..lambda_n (lambda_0 = 0 by construction)."""32 pi = frac_diff_coefs(d, n + 1)33 c = amp * pi; c[0] = 1.0 + amp * (pi[0] - 1.0) # 1 + amp((1-L)^d - 1)34 # (1 - phi L) c(L)35 a = c.copy(); a[1:] -= phi * c[:-1]36 # divide by (1 - beta L): b_k = a_k + beta b_{k-1}37 b = np.empty_like(a); b[0] = a[0]38 for k in range(1, len(a)):39 b[k] = a[k] + beta * b[k - 1]40 lam = -b[1:n + 1]41 return lam424344def fiegarch_weights(d: float, phi: float, psi: float, n: int = TRUNC) -> np.ndarray:45 """lambda_k, k = 0..n-1, of (1 - phi L)^{-1} (1 - L)^{-d} (1 + psi L)."""46 pi = frac_diff_coefs(-d, n) # (1-L)^{-d}47 a = pi.copy(); a[1:] += psi * pi[:-1] # (1 + psi L)(1-L)^{-d}48 b = np.empty_like(a); b[0] = a[0]49 for k in range(1, n):50 b[k] = a[k] + phi * b[k - 1]51 return b525354@njit(cache=True)55def rec_archinf(x, lam, const, s0):56 """s_t = const + sum_{k=1}^{K} lam_k x_{t-k}, with x_{t-k} = s0 (backcast) for t-k < 0."""57 T = x.shape[0]; K = lam.shape[0]58 out = np.empty(T)59 for t in range(T):60 acc = const61 for k in range(1, K + 1):62 j = t - k63 acc += lam[k - 1] * (x[j] if j >= 0 else s0)64 out[t] = acc65 return out666768@njit(cache=True)69def rec_fiaparch(eps, lam, const, delta, gamma, s0d):70 T = eps.shape[0]; K = lam.shape[0]71 sd = np.empty(T) # sigma^delta72 g = np.empty(T) # (|eps| - gamma eps)^delta73 for t in range(T):74 acc = const75 for k in range(1, K + 1):76 j = t - k77 acc += lam[k - 1] * (g[j] if j >= 0 else s0d)78 if acc < 1e-12:79 acc = 1e-1280 sd[t] = acc81 g[t] = (abs(eps[t]) - gamma * eps[t]) ** delta82 return sd838485@njit(cache=True)86def rec_fiegarch(eps, lam, omega, theta, gamma, ez, ls0):87 T = eps.shape[0]; K = lam.shape[0]88 ls = np.empty(T); gz = np.empty(T)89 for t in range(T):90 acc = omega91 for k in range(0, K):92 j = t - 1 - k93 acc += lam[k] * (gz[j] if j >= 0 else 0.0)94 if acc > 50.0:95 acc = 50.096 if acc < -50.0:97 acc = -50.098 ls[t] = acc99 z = eps[t] / math.exp(0.5 * acc)100 gz[t] = theta * z + gamma * (abs(z) - ez)101 return ls102103104@njit(cache=True)105def rec_stgarch(eps, omega, a1, a2, gam, beta, s0):106 T = eps.shape[0]; s = np.empty(T); s[0] = s0107 for t in range(1, T):108 e = eps[t - 1]109 a = -gam * e110 a = 50.0 if a > 50.0 else (-50.0 if a < -50.0 else a)111 F = 1.0 / (1.0 + math.exp(a))112 s[t] = omega + (a1 + a2 * F) * e * e + beta * s[t - 1]113 if s[t] < 1e-12:114 s[t] = 1e-12115 return s116117118@njit(cache=True)119def rec_betategarch(eps, omega, alpha, gamma, beta, nu, ls0):120 T = eps.shape[0]; ls = np.empty(T); ls[0] = ls0121 for t in range(1, T):122 s2 = math.exp(ls[t - 1]); z2 = eps[t - 1] * eps[t - 1] / s2123 u = (nu + 1.0) * z2 / (nu - 2.0 + z2) - 1.0124 sgn = -1.0 if eps[t - 1] > 0 else (1.0 if eps[t - 1] < 0 else 0.0)125 v = omega + beta * ls[t - 1] + alpha * u + gamma * sgn * (u + 1.0)126 ls[t] = min(max(v, -50.0), 50.0)127 return ls128129130@njit(cache=True)131def rec_realegarch(eps, logx, omega, beta, tau1, tau2, gamma, xi, phi, d1, d2, ls0):132 # log s2_t = omega + beta log s2_{t-1} + tau(z_{t-1}) + gamma u_{t-1}; log x_t = xi + phi log s2_t + delta(z_t) + u_t133 T = eps.shape[0]; ls = np.empty(T); u = np.empty(T)134 for t in range(T):135 if t == 0:136 ls[t] = ls0137 else:138 zp = eps[t - 1] / math.exp(0.5 * ls[t - 1])139 v = omega + beta * ls[t - 1] + tau1 * zp + tau2 * (zp * zp - 1.0) + gamma * u[t - 1]140 ls[t] = min(max(v, -50.0), 50.0)141 z = eps[t] / math.exp(0.5 * ls[t])142 u[t] = logx[t] - xi - phi * ls[t] - d1 * z - d2 * (z * z - 1.0)143 return ls, u144145146@njit(cache=True)147def rec_garch_midas(eps, tau, omega_g, alpha, beta, g0):148 # g_t = omega_g + alpha eps_{t-1}^2 / tau_t + beta g_{t-1}; sigma2_t = tau_t g_t (omega_g = 1 - alpha - beta);149 # tau_t is the long-run component of the month containing t, built from previous months only (EGS 2013, eq. 4)150 T = eps.shape[0]; g = np.empty(T); g[0] = g0151 for t in range(1, T):152 g[t] = omega_g + alpha * eps[t - 1] * eps[t - 1] / tau[t] + beta * g[t - 1]153 return g154155156def _simulate_paths(step_fn, h, n, rng, dist, dp):157 """Generic Monte Carlo multi-step: step_fn(state, z) -> (sigma2_next, state_next) vectorized over n paths."""158 out = np.empty(h)159 state = None160 for k in range(h):161 z = rvs(dist, dp, n, rng) if k > 0 else None162 s2, state = step_fn(state, z)163 out[k] = float(np.mean(s2))164 return out165166167# =============================================================================== models168class STGARCH11(CustomModel):169 """Logistic smooth-transition GARCH (González-Rivera 1998; Hagerud 1997):170 sigma2_t = omega + [alpha1 + alpha2 F(eps_{t-1})] eps2_{t-1} + beta sigma2_{t-1}, F(e) = 1/(1+exp(-gamma e)), gamma > 0.171 Nests GARCH(1,1) (alpha2 = 0) and approaches GJR as gamma -> infinity (alpha2 = -gamma_GJR with F -> 1{e>0}).172 Constraints: omega > 0, alpha1 >= 0, alpha1 + alpha2 >= 0, beta >= 0, alpha1 + alpha2/2 + beta < 1 (symmetric-law persistence).173 Multi-step by simulation (the transition is nonlinear in eps)."""174 model_id = "stgarch11"; family = "regime-switching"; vol_names = ["omega", "alpha1", "alpha2", "gamma", "beta"]; analytic_multistep = True175176 def theta_to_vol(self, th):177 s, w = _sig(th[1]), _sig(th[2]) # persistence s = a1 + a2/2 + beta, share w for (a1 + a2/2)178 abar = s * w; beta = s * (1 - w)179 q = math.tanh(th[3]) # a2/2 in (-abar, abar) -> a1 = abar - a2/2 >= 0 and a1 + a2 >= 0180 a2 = 2 * q * abar; a1 = abar - a2 / 2181 return np.array([math.exp(th[0]), a1, a2, math.exp(th[4]), beta])182183 def vol_to_theta(self, v):184 abar = v[1] + v[2] / 2; s = abar + v[4]185 return np.array([math.log(v[0]), _logit(s), _logit(abar / s), math.atanh(min(max(v[2] / (2 * abar), -0.99), 0.99)), math.log(v[3])])186187 def start_vol(self, eps, s0):188 var = float(np.var(eps))189 return [np.array([var * 0.05, 0.03, 0.08, 1.0, 0.88]), np.array([var * 0.03, 0.08, 0.0, 0.5, 0.88]), np.array([var * 0.1, 0.02, 0.12, 3.0, 0.8])]190191 def sigma2(self, v, eps, x, s0, kappa):192 return rec_stgarch(eps, v[0], v[1], v[2], v[3], v[4], s0)193194 def admissible(self, v, kappa):195 return v[0] > 0 and v[1] >= 0 and v[1] + v[2] >= 0 and v[4] >= 0 and v[1] + v[2] / 2 + v[4] < 1 and v[3] > 0196197 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):198 """One step exact; for k >= 2 and a symmetric innovation law, E[F(eps) eps^2 | sigma^2] = sigma^2 / 2 because199 F(e) + F(-e) = 1 (logistic), so the recursion is sigma2_{t+k} = omega + (alpha1 + alpha2/2 + beta) sigma2_{t+k-1}."""200 F = 1.0 / (1.0 + math.exp(min(max(-v[3] * eps[-1], -50.0), 50.0)))201 out = np.empty(h); out[0] = max(v[0] + (v[1] + v[2] * F) * eps[-1] ** 2 + v[4] * s2[-1], 1e-12)202 pers = v[1] + v[2] / 2 + v[4]203 for k in range(1, h):204 out[k] = v[0] + pers * out[k - 1]205 return out206207208class BetaTEGARCH11(CustomModel):209 """Beta-t-EGARCH (Harvey & Chakravarty 2008; Harvey 2013), written for the unit-variance Student-t:210 log sigma2_t = omega + beta log sigma2_{t-1} + alpha u_{t-1} + gamma sgn(-eps_{t-1})(u_{t-1} + 1),211 u_t = (nu+1) z_t^2 / (nu - 2 + z_t^2) - 1 (the score of the t density w.r.t. log sigma2, bounded in (-1, nu]).212 Stationary iff |beta| < 1; the bounded score makes the model robust to outliers (an EGARCH with Student-t innovations213 responds to |z|, which is unbounded). Multi-step by simulation."""214 model_id = "betategarch11"; family = "score-driven"; vol_names = ["omega", "alpha", "gamma", "beta"]; supports = ("t",); analytic_multistep = False215216 def theta_to_vol(self, th):217 return np.array([th[0], math.exp(th[1]), th[2], math.tanh(th[3])])218219 def vol_to_theta(self, v):220 return np.array([v[0], math.log(max(v[1], 1e-6)), v[2], math.atanh(min(max(v[3], -0.999), 0.999))])221222 def start_vol(self, eps, s0):223 lv = math.log(float(np.var(eps)))224 return [np.array([0.02 * lv, 0.05, 0.02, 0.98]), np.array([0.05 * lv, 0.1, 0.04, 0.95]), np.array([0.01 * lv, 0.03, 0.0, 0.99])]225226 def _negll(self, theta, r, x, dist, s0, cache):227 mu = theta[0]; v = self.theta_to_vol(theta[1:-1]); dp = _dist_theta_to_params("t", theta[-1:])228 eps = r - mu229 ls = rec_betategarch(eps, v[0], v[1], v[2], v[3], float(dp[0]), math.log(s0))230 s2 = np.exp(ls)231 if not np.all(np.isfinite(s2)):232 return _BIG233 val = -float(np.sum(loglik_z("t", eps / np.sqrt(s2), dp) - 0.5 * ls))234 return val if np.isfinite(val) else _BIG235236 def sigma2(self, v, eps, x, s0, kappa):237 return np.exp(rec_betategarch(eps, v[0], v[1], v[2], v[3], self._nu, math.log(s0)))238239 def filter(self, fit, r, x=None):240 mu, v, dp = self._unpack(fit); self._nu = float(dp[0])241 return np.exp(rec_betategarch(np.asarray(r, float) - mu, v[0], v[1], v[2], v[3], self._nu, math.log(fit.aux["backcast"])))242243 def admissible(self, v, kappa):244 return abs(v[3]) < 1 and v[1] >= 0245246 def forecast(self, fit, r, x, h, seed=0, n_sim=2000):247 from .base import Forecast248 mu, v, dp = self._unpack(fit); nu = float(dp[0])249 eps = np.asarray(r, float) - mu250 ls = rec_betategarch(eps, v[0], v[1], v[2], v[3], nu, math.log(fit.aux["backcast"]))251 rng = np.random.default_rng(seed); n = n_sim252253 def upd(ls_prev, e_prev):254 z2 = e_prev ** 2 / np.exp(ls_prev); u = (nu + 1) * z2 / (nu - 2 + z2) - 1255 return np.clip(v[0] + v[3] * ls_prev + v[1] * u + v[2] * np.sign(-e_prev) * (u + 1), -50, 50)256 ls1 = upd(np.full(n, ls[-1]), np.full(n, eps[-1]))257 out = np.empty(h); out[0] = float(np.mean(np.exp(ls1))); cur = ls1258 for k in range(1, h):259 e = np.exp(0.5 * cur) * rvs("t", dp, n, rng); cur = upd(cur, e); out[k] = float(np.mean(np.exp(cur)))260 return Forecast(self.model_id, "t", out, mu, dp, "simulation", float(np.exp(ls[-1])))261262263class RealEGARCH11(CustomModel):264 """Realized EGARCH(1,1) (Hansen & Huang 2016): log sigma2_t = omega + beta log sigma2_{t-1} + tau1 z_{t-1} + tau2 (z_{t-1}^2 - 1) + gamma u_{t-1};265 log x_t = xi + phi log sigma2_t + delta1 z_t + delta2 (z_t^2 - 1) + u_t, u_t ~ N(0, sigma_u^2). The measurement error u266 enters the volatility equation directly (unlike Realized GARCH, where log x_{t-1} enters). Stationary iff |beta| < 1.267 Joint likelihood; multi-step by simulation."""268 model_id = "realegarch11"; family = "realized"; needs_rm = True; analytic_multistep = False; supports = ("normal", "t")269 vol_names = ["omega", "beta", "tau1", "tau2", "gamma", "xi", "phi", "delta1", "delta2", "log_sigma_u"]270271 def theta_to_vol(self, th):272 v = np.array(th, float); v[1] = math.tanh(th[1]); return v273274 def vol_to_theta(self, v):275 th = np.array(v, float); th[1] = math.atanh(min(max(v[1], -0.999), 0.999)); return th276277 def start_vol(self, eps, s0):278 lx = float(np.mean(np.log(np.maximum(self._x, 1e-8))))279 return [np.array([0.02 * lx, 0.98, -0.03, 0.03, 0.1, 0.0, 1.0, -0.05, 0.05, math.log(0.4)]),280 np.array([0.05 * lx, 0.95, -0.05, 0.05, 0.2, 0.0, 0.9, -0.05, 0.05, math.log(0.5)]),281 np.array([0.01 * lx, 0.99, -0.02, 0.02, 0.05, 0.0, 1.0, -0.02, 0.02, math.log(0.35)])]282283 def fit(self, r, x=None, dist="normal", max_restarts=3, maxiter=3000):284 self._x = np.asarray(x, float); return super().fit(r, x, dist, max_restarts, maxiter)285286 def sigma2(self, v, eps, x, s0, kappa):287 logx = np.log(np.maximum(np.asarray(x, float), 1e-8))288 ls, u = rec_realegarch(eps, logx, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], math.log(s0))289 self._u = u290 return np.exp(ls)291292 def extra_loglik(self, v, eps, s2, x):293 su2 = math.exp(2 * v[9]); return -0.5 * math.log(2 * math.pi * su2) - self._u ** 2 / (2 * su2)294295 def admissible(self, v, kappa):296 return abs(v[1]) < 1297298 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):299 omega, beta, t1, t2, gam, xi, phi, d1, d2, lsu = v; su = math.exp(lsu); n = 2000300 logx = np.log(np.maximum(np.asarray(x, float), 1e-8))301 ls_all, u_all = rec_realegarch(eps, logx, omega, beta, t1, t2, gam, xi, phi, d1, d2, math.log(s2[0]))302 z_last = eps[-1] / math.sqrt(s2[-1]); u_last = u_all[-1]303 ls1 = omega + beta * math.log(s2[-1]) + t1 * z_last + t2 * (z_last ** 2 - 1) + gam * u_last304 out = np.empty(h); out[0] = math.exp(ls1); cur = np.full(n, ls1)305 for k in range(1, h):306 z = rvs(dist, dp, n, rng); u = rng.standard_normal(n) * su307 cur = np.clip(omega + beta * cur + t1 * z + t2 * (z * z - 1) + gam * u, -30, 30); out[k] = float(np.mean(np.exp(cur)))308 return out309310311class GARCHMIDAS(CustomModel):312 """GARCH-MIDAS (Engle, Ghysels & Sohn 2013), fixed-span version with the instrument's own monthly realized variance:313 sigma2_t = tau_t g_t; g_t = (1-alpha-beta) + alpha eps_{t-1}^2/tau_{t-1} + beta g_{t-1};314 tau_t = m + theta sum_{k=1}^{K} phi_k(w) RV^{(m)}_{t-k}, RV^{(m)} = sum of squared daily returns of calendar month,315 beta-polynomial weights phi_k(w) ~ (1 - k/K)^{w-1} (w1 = 1), K = 24 months, tau fixed within the month.316 The information set is the instrument's own daily returns (Tier A); tau positive by m > 0, theta >= 0.317 Nests GARCH(1,1) when theta = 0 (then tau = m is the unconditional variance). Analytic multi-step with tau frozen."""318 model_id = "garch_midas"; family = "component"; vol_names = ["m", "theta", "w", "alpha", "beta"]; K = 24319320 def theta_to_vol(self, th):321 s, wsh = _sig(th[3]), _sig(th[4])322 return np.array([math.exp(th[0]), math.exp(th[1]), 1.0 + math.exp(th[2]), s * wsh, s * (1 - wsh)])323324 def vol_to_theta(self, v):325 s = v[3] + v[4]326 return np.array([math.log(v[0]), math.log(max(v[1], 1e-6)), math.log(max(v[2] - 1.0, 1e-6)), _logit(s), _logit(v[3] / s)])327328 def start_vol(self, eps, s0):329 var = float(np.var(eps))330 return [np.array([var * 0.5, 0.5, 3.0, 0.06, 0.90]), np.array([var * 0.8, 0.2, 5.0, 0.04, 0.93]), np.array([var * 0.3, 0.8, 2.0, 0.08, 0.85])]331332 def _month_rv(self, r):333 """Monthly RV (sum of squared daily deviations) using the day index -> month map set in fit()."""334 rv = np.zeros(self._n_months)335 np.add.at(rv, self._month_of_day, (r - r.mean()) ** 2)336 return rv337338 def _tau(self, v, rv_m):339 K = self.K; k = np.arange(1, K + 1); w = v[2]340 phi = (1 - k / K) ** (w - 1); phi = phi / phi.sum()341 tau_m = np.full(self._n_months, np.nan)342 for mth in range(self._n_months):343 if mth >= K:344 tau_m[mth] = v[0] + v[1] * float(np.dot(phi, rv_m[mth - K:mth][::-1]))345 elif mth >= 6: # shorter history at the start of the window: use the available months (renormalized)346 kk = mth; ph = (1 - np.arange(1, kk + 1) / K) ** (w - 1); ph /= ph.sum()347 tau_m[mth] = v[0] + v[1] * float(np.dot(ph, rv_m[:mth][::-1]))348 else:349 tau_m[mth] = v[0] + v[1] * float(np.mean(rv_m[:max(mth, 1)]))350 return np.maximum(tau_m[self._month_of_day], 1e-10)351352 def fit(self, r, x=None, dist="normal", max_restarts=3, maxiter=2000, months=None):353 r = np.asarray(r, float)354 mo = np.asarray(months if months is not None else self._months_from_len(len(r)))355 _, self._month_of_day = np.unique(mo, return_inverse=True); self._n_months = int(self._month_of_day.max()) + 1356 self._rv_m = self._month_rv(r)357 f = super().fit(r, x, dist, max_restarts, maxiter)358 f.aux["months"] = [int(m) for m in mo]359 return f360361 @staticmethod362 def _months_from_len(n):363 return np.arange(n) // 21 # fallback when no calendar is supplied (21 trading days per month)364365 def sigma2(self, v, eps, x, s0, kappa):366 tau = self._tau(v, self._rv_m)367 g = rec_garch_midas(eps, tau, 1 - v[3] - v[4], v[3], v[4], 1.0)368 self._g_last, self._tau_last = g[-1], tau[-1]369 return tau * g370371 def admissible(self, v, kappa):372 return v[0] > 0 and v[1] >= 0 and v[2] >= 1 and v[3] >= 0 and v[4] >= 0 and v[3] + v[4] < 1373374 def filter(self, fit, r, x=None):375 r = np.asarray(r, float); mo = fit.aux.get("months")376 mo = np.asarray(mo) if mo is not None else self._months_from_len(len(r))377 self._month_of_day = np.unique(mo, return_inverse=True)[1]378 self._n_months = int(self._month_of_day.max()) + 1; self._rv_m = self._month_rv(r)379 mu, v, dp = self._unpack(fit)380 return self.sigma2(v, r - mu, None, fit.aux["backcast"], 0.5)381382 def forecast(self, fit, r, x, h, seed=0, n_sim=2000):383 from .base import Forecast384 s2 = self.filter(fit, r, x) # sets _rv_m, _tau_last, _g_last from the window385 mu, v, dp = self._unpack(fit)386 sv = self.step_forecast(v, np.asarray(r, float) - mu, s2, x, h, 0.5, dp, fit.dist, None)387 return Forecast(self.model_id, fit.dist, np.asarray(sv, float), mu, dp, "analytic", float(s2[-1]))388389 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):390 tau = self._tau_last; g = self._g_last391 out = np.empty(h); e2 = eps[-1] ** 2392 g1 = (1 - v[3] - v[4]) + v[3] * e2 / tau + v[4] * g; out[0] = tau * g1393 for k in range(1, h):394 g1 = (1 - v[3] - v[4]) + (v[3] + v[4]) * g1; out[k] = tau * g1395 return out396397398class HYGARCH11(CustomModel):399 """Hyperbolic GARCH (Davidson 2004): sigma2_t = omega/(1-beta) + lambda(L) eps2_t with400 lambda(L) = 1 - (1-beta L)^{-1}(1-phi L)[1 + amp((1-L)^d - 1)]; amp = 1 gives FIGARCH(1,d,1), amp < 1 a covariance-stationary401 process (amplitude parameter), amp > 1 an explosive one. ARCH(infinity) truncated at 1000 lags (as arch's FIGARCH).402 Admissible iff all truncated weights lambda_k >= 0 and omega > 0 (checked numerically). Analytic multi-step through the weights."""403 model_id = "hygarch11"; family = "fractional"; vol_names = ["omega", "d", "phi", "beta", "amp"]404405 def theta_to_vol(self, th):406 return np.array([math.exp(th[0]), _sig(th[1]), _sig(th[2]), _sig(th[3]), math.exp(th[4])])407408 def vol_to_theta(self, v):409 return np.array([math.log(v[0]), _logit(v[1]), _logit(v[2]), _logit(v[3]), math.log(max(v[4], 1e-6))])410411 def start_vol(self, eps, s0):412 var = float(np.var(eps))413 return [np.array([var * 0.1, 0.4, 0.2, 0.5, 1.0]), np.array([var * 0.05, 0.6, 0.1, 0.4, 0.9]), np.array([var * 0.2, 0.3, 0.3, 0.6, 1.1])]414415 def _weights(self, v):416 return figarch_weights(v[1], v[2], v[3], v[4], TRUNC)417418 def sigma2(self, v, eps, x, s0, kappa):419 lam = self._weights(v)420 if np.any(lam < -1e-10):421 raise FloatingPointError("negative ARCH(inf) weight")422 return rec_archinf(eps * eps, lam, v[0] / (1 - v[3]), s0)423424 def admissible(self, v, kappa):425 lam = self._weights(v)426 return v[0] > 0 and bool(np.all(lam >= -1e-10)) and v[4] < 1.0 # amp < 1: covariance stationary (Davidson 2004)427428 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):429 lam = self._weights(v); const = v[0] / (1 - v[3]); K = len(lam)430 hist = list((eps * eps)[-K:][::-1]) # most recent first431 if len(hist) < K:432 hist += [s2[0]] * (K - len(hist))433 out = np.empty(h)434 for k in range(h):435 val = const + float(np.dot(lam, np.asarray(hist[:K])))436 out[k] = max(val, 1e-12); hist.insert(0, out[k])437 return out438439440class FIAPARCH11(CustomModel):441 """FIAPARCH(1,d,1) (Tse 1998): sigma_t^delta = omega/(1-beta) + lambda(L)(|eps_t| - gamma eps_t)^delta with FIGARCH weights442 (amp = 1); delta = 2, gamma = 0 gives FIGARCH. |gamma| < 1, delta in (0.5, 4). Multi-step by simulation (power transform)."""443 model_id = "fiaparch11"; family = "fractional"; vol_names = ["omega", "d", "phi", "beta", "gamma", "delta"]; analytic_multistep = False444445 def theta_to_vol(self, th):446 return np.array([math.exp(th[0]), _sig(th[1]), _sig(th[2]), _sig(th[3]), 0.99 * math.tanh(th[4]), 0.5 + 3.5 * _sig(th[5])])447448 def vol_to_theta(self, v):449 return np.array([math.log(v[0]), _logit(v[1]), _logit(v[2]), _logit(v[3]), math.atanh(v[4] / 0.99), _logit((v[5] - 0.5) / 3.5)])450451 def start_vol(self, eps, s0):452 var = float(np.var(eps))453 return [np.array([var * 0.1, 0.4, 0.2, 0.5, 0.1, 2.0]), np.array([var ** 0.75 * 0.1, 0.5, 0.1, 0.4, 0.2, 1.5]), np.array([var * 0.2, 0.3, 0.3, 0.6, 0.0, 2.0])]454455 def sigma2(self, v, eps, x, s0, kappa):456 lam = figarch_weights(v[1], v[2], v[3], 1.0, TRUNC)457 if np.any(lam < -1e-10):458 raise FloatingPointError("negative ARCH(inf) weight")459 sd = rec_fiaparch(eps, lam, v[0] / (1 - v[3]), v[5], v[4], s0 ** (v[5] / 2))460 return sd ** (2.0 / v[5])461462 def admissible(self, v, kappa):463 return v[0] > 0 and bool(np.all(figarch_weights(v[1], v[2], v[3], 1.0, TRUNC) >= -1e-10))464465 def step_forecast(self, v, eps, s2, x, h, kappa, dp, dist, rng):466 lam = figarch_weights(v[1], v[2], v[3], 1.0, TRUNC); const = v[0] / (1 - v[3]); K = len(lam); n = 2000; delta, gam = v[5], v[4]467 g_hist = (np.abs(eps) - gam * eps) ** delta468 hist = np.tile(np.concatenate([g_hist[-K:][::-1], np.full(max(0, K - len(g_hist)), s2[0] ** (delta / 2))])[:K], (n, 1))469 out = np.empty(h)470 for k in range(h):471 sd = np.maximum(const + hist @ lam, 1e-12) # sigma^delta per path472 s2k = sd ** (2.0 / delta); out[k] = float(np.mean(s2k))473 z = rvs(dist, dp, n, rng); e = np.sqrt(s2k) * z474 hist = np.column_stack([(np.abs(e) - gam * e) ** delta, hist[:, :-1]])475 return out476477478class FIEGARCH11(CustomModel):479 """FIEGARCH(1,d,1) (Bollerslev & Mikkelsen 1996): log sigma2_t = omega + (1-phi L)^{-1}(1-L)^{-d}(1+psi L) g(z_{t-1}),480 g(z) = theta z + gamma(|z| - E|z|); weights truncated at 1000 lags; |phi| < 1, |psi| < 1, d < 0.5 for covariance stationarity of481 log sigma2 (admissible flag). d = 0 gives Nelson's EGARCH(1,1) up to the parameterization. Multi-step by simulation."""482 model_id = "fiegarch11"; family = "fractional"; vol_names = ["omega", "d", "phi", "psi", "theta", "gamma"]; analytic_multistep = False483484 def theta_to_vol(self, th):485 return np.array([th[0], 0.499 * _sig(th[1]), 0.99 * math.tanh(th[2]), 0.99 * math.tanh(th[3]), th[4], math.exp(th[5])])486487 def vol_to_theta(self, v):488 return np.array([v[0], _logit(v[1] / 0.499), math.atanh(v[2] / 0.99), math.atanh(v[3] / 0.99), v[4], math.log(max(v[5], 1e-6))])489490 def start_vol(self, eps, s0):491 lv = math.log(float(np.var(eps)))492 return [np.array([lv * 0.05, 0.3, 0.9, 0.0, -0.05, 0.15]), np.array([lv * 0.02, 0.45, 0.95, -0.2, -0.08, 0.1]), np.array([lv * 0.1, 0.2, 0.8, 0.2, 0.0, 0.2])]493494 def _ez(self, dist, dp):495 if dist == "normal":496 return math.sqrt(2 / math.pi)497 from .distributions import ppf498 u = (np.arange(20000) + 0.5) / 20000; return float(np.mean(np.abs(ppf(dist, u, dp))))499500 def _negll(self, theta, r, x, dist, s0, cache):501 kd = DIST_NPARAMS[dist]; mu = theta[0]; v = self.theta_to_vol(theta[1:len(theta) - kd]); dp = _dist_theta_to_params(dist, theta[len(theta) - kd:]) if kd else np.empty(0)502 key = tuple(np.round(dp, 6)); ez = cache.get(key)503 if ez is None:504 ez = self._ez(dist, dp); cache[key] = ez505 eps = r - mu; lam = fiegarch_weights(v[1], v[2], v[3], TRUNC)506 ls = rec_fiegarch(eps, lam, v[0], v[4], v[5], ez, math.log(s0))507 s2 = np.exp(ls)508 val = -float(np.sum(loglik_z(dist, eps / np.sqrt(s2), dp) - 0.5 * ls))509 return val if np.isfinite(val) else _BIG510511 def sigma2(self, v, eps, x, s0, kappa):512 lam = fiegarch_weights(v[1], v[2], v[3], TRUNC)513 return np.exp(rec_fiegarch(eps, lam, v[0], v[4], v[5], self._ez_cached, math.log(s0)))514515 def filter(self, fit, r, x=None):516 mu, v, dp = self._unpack(fit); self._ez_cached = self._ez(fit.dist, dp)517 return self.sigma2(v, np.asarray(r, float) - mu, None, fit.aux["backcast"], 0.5)518519 def admissible(self, v, kappa):520 return v[1] < 0.5 and abs(v[2]) < 1 and abs(v[3]) < 1521522 def forecast(self, fit, r, x, h, seed=0, n_sim=2000):523 from .base import Forecast524 mu, v, dp = self._unpack(fit); ez = self._ez(fit.dist, dp); self._ez_cached = ez525 eps = np.asarray(r, float) - mu; lam = fiegarch_weights(v[1], v[2], v[3], TRUNC); K = len(lam)526 ls = rec_fiegarch(eps, lam, v[0], v[4], v[5], ez, math.log(fit.aux["backcast"]))527 z_all = eps / np.exp(0.5 * ls); g_all = v[4] * z_all + v[5] * (np.abs(z_all) - ez)528 n = n_sim; rng = np.random.default_rng(seed)529 hist = np.tile(np.concatenate([g_all[-K:][::-1], np.zeros(max(0, K - len(g_all)))])[:K], (n, 1))530 out = np.empty(h)531 for k in range(h):532 lsk = np.clip(v[0] + hist @ lam, -50, 50); out[k] = float(np.mean(np.exp(lsk)))533 z = rvs(fit.dist, dp, n, rng); g = v[4] * z + v[5] * (np.abs(z) - ez)534 hist = np.column_stack([g, hist[:, :-1]])535 return Forecast(self.model_id, fit.dist, out, mu, dp, "simulation", float(np.exp(ls[-1])))536537538EXTENSIONS = {m.model_id: m for m in [STGARCH11, BetaTEGARCH11, RealEGARCH11, GARCHMIDAS, HYGARCH11, FIAPARCH11, FIEGARCH11]}539