Repository Public Research preview
1"""Models estimated with the `arch` package (Sheppard, version-locked 8.0.0): constant variance, ARCH(q), GARCH(p,q),2GJR-GARCH, TARCH/ZARCH (Zakoïan threshold on the conditional standard deviation), AVGARCH (Taylor/Schwert),3APARCH (Ding, Granger & Engle), EGARCH (Nelson), FIGARCH (Baillie, Bollerslev & Mikkelsen; arch's truncation 1000 lags),4EWMA (RiskMetrics 1994, lambda fixed 0.94), RiskMetrics2006 (Zumbach). Distributions: normal, t, ged, skewt.56Multi-step variance forecasts: arch's analytic recursion where it exists (ARCH, GARCH/GJR/TARCH/AVGARCH/APARCH with7power 2 only, FIGARCH, EWMA, RM2006); otherwise simulation with 2,000 antithetic paths and a fixed seed8(EGARCH, APARCH, TARCH/AVGARCH for h > 1). The method actually used is recorded per forecast.9"""10from __future__ import annotations1112import time13import warnings1415import numpy as np16from arch.univariate import (APARCH, ARCH, EGARCH, EWMAVariance, FIGARCH, GARCH, ConstantMean, ConstantVariance,17 GeneralizedError, Normal, RiskMetrics2006, SkewStudent, StudentsT)1819from .base import FitResult, Forecast2021DISTS = {"normal": Normal, "t": StudentsT, "ged": GeneralizedError, "skewt": SkewStudent}2223SPECS: dict[str, dict] = {24 "constvar": dict(vol=lambda: ConstantVariance(), analytic=True),25 "arch5": dict(vol=lambda: ARCH(p=5), analytic=True),26 "garch11": dict(vol=lambda: GARCH(p=1, o=0, q=1), analytic=True),27 "garch21": dict(vol=lambda: GARCH(p=2, o=0, q=1), analytic=True),28 "garch12": dict(vol=lambda: GARCH(p=1, o=0, q=2), analytic=True),29 "gjr11": dict(vol=lambda: GARCH(p=1, o=1, q=1), analytic=True),30 "tarch11": dict(vol=lambda: GARCH(p=1, o=1, q=1, power=1.0), analytic=False),31 "avgarch11": dict(vol=lambda: GARCH(p=1, o=0, q=1, power=1.0), analytic=False),32 "aparch11": dict(vol=lambda: APARCH(p=1, o=1, q=1), analytic=False),33 "egarch11": dict(vol=lambda: EGARCH(p=1, o=1, q=1), analytic=False),34 "egarch11c": dict(vol=lambda: EGARCHConstrained(p=1, o=1, q=1), analytic=False),35 "figarch11": dict(vol=lambda: FIGARCH(p=1, q=1), analytic=True),36 "ewma094": dict(vol=lambda: EWMAVariance(lam=0.94), analytic=True, fixed=True),37 "rm2006": dict(vol=lambda: RiskMetrics2006(), analytic=True, fixed=True),38}39N_SIM = 2000404142class EGARCHConstrained(EGARCH):43 """EGARCH(1,1) with the sign/stationarity constraints of the textbook parameterization (A-06, extension 'egarch11c'):44 alpha >= 0 (volatility rises with |z|), 0 <= beta <= 0.995. Unconstrained arch fits with alpha < 0 and beta -> 1 are45 locally well-fitting but dynamically explosive out of sample (large sigma -> tiny |z| -> larger sigma)."""4647 def bounds(self, resids):48 b = super().bounds(resids)49 names = self.parameter_names()50 out = []51 for n, (lo, hi) in zip(names, b):52 if n.startswith("alpha"):53 out.append((0.0, hi))54 elif n.startswith("beta"):55 out.append((0.0, 0.995))56 else:57 out.append((lo, hi))58 return out596061def _scalar_backcast(am, res):62 try:63 b = am.volatility.backcast(res.resid)64 return float(b) if np.ndim(b) == 0 else None65 except Exception: # noqa: BLE00166 return None676869def _nit(res) -> int:70 o = getattr(res, "optimization_result", None)71 if o is None:72 return 073 for k in ("nit", "nfev"):74 try:75 return int(o[k]) if isinstance(o, dict) else int(getattr(o, k))76 except (KeyError, AttributeError, TypeError):77 continue78 return 0798081class ArchModel:82 def __init__(self, model_id: str):83 self.model_id = model_id84 self.spec = SPECS[model_id]85 self.family = {"constvar": "baseline", "ewma094": "baseline", "rm2006": "baseline", "arch5": "classical",86 "garch11": "classical", "garch21": "classical", "garch12": "classical", "gjr11": "asymmetric",87 "tarch11": "asymmetric", "avgarch11": "power", "aparch11": "power", "egarch11": "asymmetric", "egarch11c": "asymmetric",88 "figarch11": "fractional"}[model_id]8990 def _build(self, r: np.ndarray, dist: str, seed: int | None = None):91 am = ConstantMean(np.asarray(r, float))92 am.volatility = self.spec["vol"]()93 am.distribution = DISTS[dist](seed=seed) # seeded generator -> reproducible simulation forecasts (A-01)94 return am9596 def _alt_starts(self, am, r: np.ndarray, dist: str) -> list[np.ndarray]:97 """Alternative starting vectors (A-05): arch's own volatility start with two persistence regimes, unconditional98 variance matched to the sample variance, distribution parameters at their defaults."""99 resid = np.asarray(r, float) - float(np.mean(r))100 var = float(np.var(resid))101 try:102 sv = np.asarray(am.volatility.starting_values(resid), float)103 except Exception: # noqa: BLE001104 return []105 names = list(am.volatility.parameter_names())106 dstart = {"normal": [], "t": [8.0], "ged": [1.5], "skewt": [8.0, -0.1]}[dist]107 outs = []108 for alpha, beta in ((0.10, 0.85), (0.03, 0.96), (0.10, 0.98)):109 v = sv.copy()110 na = max(1, sum(1 for m in names if m.startswith("alpha"))); nb = max(1, sum(1 for m in names if m.startswith("beta")))111 for i, n in enumerate(names):112 if n.startswith("alpha"):113 v[i] = alpha / na114 elif n.startswith("beta"):115 v[i] = beta / nb116 elif n.startswith("gamma"):117 v[i] = -0.05 if self.model_id == "egarch11" else 0.05118 if "omega" in names:119 i = names.index("omega")120 if self.model_id in ("egarch11", "egarch11c"):121 v[i] = (1 - beta) * np.log(max(var, 1e-8))122 elif self.model_id in ("tarch11", "avgarch11"):123 v[i] = np.sqrt(var) * max(1 - alpha - beta, 0.02)124 elif self.model_id == "figarch11":125 v[i] = var * 0.05126 else:127 v[i] = var * max(1 - alpha - beta, 0.02)128 outs.append(np.concatenate([[float(np.mean(r))], v, dstart]))129 return outs130131 def _reference_loglik(self, r: np.ndarray, dist: str) -> float | None:132 """GARCH(1,1) with the same innovation law: a likelihood floor for models at least as flexible (A-05)."""133 if self.model_id in ("garch11", "constvar", "arch5", "ewma094", "rm2006"):134 return None135 try:136 am = ConstantMean(np.asarray(r, float)); am.volatility = GARCH(1, 0, 1); am.distribution = DISTS[dist]()137 return float(am.fit(disp="off").loglikelihood)138 except Exception: # noqa: BLE001139 return None140141 DEGENERATE_TOL = 50.0142143 def fit(self, r: np.ndarray, x=None, dist: str = "normal", max_restarts: int = 2, maxiter: int = 1000) -> FitResult:144 t0 = time.perf_counter()145 am = self._build(r, dist)146 with warnings.catch_warnings(record=True) as w:147 warnings.simplefilter("always")148 # A-05: multi-start — arch default plus two alternative starting vectors; keep the best converged optimum149 cands = []150 n_iter, n_rest = 0, 0151 for k, sv in enumerate([None] + (self._alt_starts(am, r, dist) if not self.spec.get("fixed") else [])):152 try:153 res_k = am.fit(disp="off", starting_values=sv, options={"maxiter": maxiter})154 except Exception: # noqa: BLE001155 continue156 n_iter += _nit(res_k); n_rest += (k > 0)157 cands.append(res_k)158 if not cands:159 raise RuntimeError("arch fit failed for every starting value")160 good = [c for c in cands if c.convergence_flag == 0] or cands161 res = max(good, key=lambda c: c.loglikelihood)162 conv = res.convergence_flag == 0163 # A-05: likelihood floor — a flexible model far below GARCH(1,1) with the same law is a degenerate optimum164 ref = self._reference_loglik(r, dist)165 degenerate = ref is not None and np.isfinite(ref) and res.loglikelihood < ref - self.DEGENERATE_TOL166 if degenerate:167 for i in range(max_restarts):168 sv = res.params.values * (0.8 + 0.4 * np.random.default_rng(i).uniform(size=len(res.params)))169 try:170 res2 = am.fit(disp="off", starting_values=sv, options={"maxiter": maxiter * 2})171 except Exception: # noqa: BLE001172 continue173 n_rest += 1; n_iter += _nit(res2)174 if res2.loglikelihood > res.loglikelihood:175 res = res2; conv = res.convergence_flag == 0176 if res.loglikelihood >= ref - self.DEGENERATE_TOL:177 degenerate = False178 break179 if degenerate:180 conv = False # reported as a failure (status 'failed'); the driver keeps the last good fit181 notes = ";".join(sorted({str(x.message)[:80] for x in w if "ConvergenceWarning" in str(x.category) or "DataScaleWarning" in str(x.category)}))182 if degenerate:183 notes = (notes + ";" if notes else "") + f"degenerate: loglik {res.loglikelihood:.1f} < GARCH(1,1) ref {ref:.1f} - {self.DEGENERATE_TOL:.0f}"184 params = res.params.values.astype(float)185 names = list(res.params.index)186 adm = True187 vp = self._vol_params(names, params)188 if self.model_id in ("garch11", "garch21", "garch12", "arch5"):189 adm = sum(v for k, v in vp.items() if k.startswith("alpha") or k.startswith("beta")) < 1.0190 elif self.model_id == "gjr11":191 adm = vp.get("alpha[1]", 0) + vp.get("beta[1]", 0) + 0.5 * vp.get("gamma[1]", 0) < 1.0192 elif self.model_id == "figarch11":193 adm = 0.0 <= vp.get("d", 0) <= 1.0194 elif self.model_id in ("egarch11", "egarch11c"):195 adm = abs(vp.get("beta[1]", 0)) < 1.0 and vp.get("alpha[1]", 0) >= 0.0 # A-06: alpha < 0 flagged (explosive feedback)196 return FitResult(self.model_id, dist, names, params, float(res.loglikelihood), bool(conv), int(res.nobs), n_iter,197 n_rest, time.perf_counter() - t0, bool(adm), notes,198 {"backcast": _scalar_backcast(am, res),199 "arch_version": "8.0.0"})200201 @staticmethod202 def _vol_params(names, params):203 return {n: float(p) for n, p in zip(names, params) if n not in ("mu", "nu", "eta", "lambda")}204205 def filter(self, fit: FitResult, r: np.ndarray, x=None) -> np.ndarray:206 am = self._build(r, fit.dist)207 res = am.fix(fit.params)208 return np.asarray(res.conditional_volatility, float) ** 2209210 def forecast(self, fit: FitResult, r: np.ndarray, x, h: int, seed: int = 0, n_sim: int = N_SIM) -> Forecast:211 am = self._build(r, fit.dist, seed=seed)212 res = am.fix(fit.params)213 s2 = np.asarray(res.conditional_volatility, float) ** 2214 if self.spec["analytic"] or h == 1:215 method = "analytic"216 with warnings.catch_warnings():217 warnings.simplefilter("ignore")218 fc = res.forecast(horizon=h, reindex=False, method="analytic")219 else:220 method = "simulation" # draws come from the seeded distribution generator (all innovation laws)221 fc = res.forecast(horizon=h, reindex=False, method="simulation", simulations=n_sim)222 sv = np.asarray(fc.variance.values[0], float)223 mu = float(fc.mean.values[0][0])224 kd = {"normal": 0, "t": 1, "ged": 1, "skewt": 2}[fit.dist]225 dp = fit.params[len(fit.params) - kd:] if kd else np.empty(0)226 return Forecast(self.model_id, fit.dist, sv, mu, dp, method, float(s2[-1]))227