Repository Public Research preview
1"""Tier C — covariance forecasts for a fixed basket of N daily return series (percent), h = 1.23Models (model_id): rolling_cov252 (sample covariance of the last 252 days), ewma_cov094 (RiskMetrics, lambda 0.94),4ccc_garch11 (Bollerslev 1990: univariate GARCH(1,1)-normal per asset + constant correlation of standardized residuals),5dcc11 (Engle 2002, two-step QMLE: GARCH(1,1) per asset, then DCC(1,1) correlation targeting with the sample correlation6of the standardized residuals; Q_t = (1-a-b) S + a z_{t-1} z_{t-1}' + b Q_{t-1}, R_t = diag(Q_t)^{-1/2} Q_t diag(Q_t)^{-1/2}),7adcc11 (Cappiello, Engle & Sheppard 2006: + g n_{t-1} n_{t-1}' with n = z 1{z<0}, targeting (1-a-b) S - g N with8N = E[n n']), scalar_bekk11 (Engle & Kroner 1995, scalar: H_t = C C' + a eps_{t-1} eps_{t-1}' + b H_{t-1} with covariance9targeting C C' = (1-a-b) Sigma_bar, Gaussian QMLE). All one-step forecasts are symmetric positive definite by construction10(checked: minimum eigenvalue > 0 is recorded). Parameters are estimated by Gaussian QMLE with scipy; standardized11residuals of the first step come from wp46.models.arch_models GARCH(1,1).12"""13from __future__ import annotations1415import math16import time1718import numpy as np19from scipy import optimize2021from .arch_models import ArchModel222324def _corr_from_cov(H):25 d = 1.0 / np.sqrt(np.diag(H))26 return H * np.outer(d, d)272829def _univariate_garch(R: np.ndarray):30 """Fit GARCH(1,1)-normal to each column; return (sigma2 matrix T x N, one-step variance forecasts, fits)."""31 T, N = R.shape32 S2 = np.empty((T, N)); f1 = np.empty(N); fits = []33 m = ArchModel("garch11")34 for j in range(N):35 f = m.fit(R[:, j], dist="normal")36 S2[:, j] = m.filter(f, R[:, j])37 f1[j] = m.forecast(f, R[:, j], None, 1).step_var[0]38 fits.append(f)39 return S2, f1, fits404142def _dcc_negll(theta, Z, S, asym, Nbar):43 a = 1 / (1 + math.exp(-theta[0])); b = (1 - a) / (1 + math.exp(-theta[1])); g = 0.044 if asym:45 g = (1 - a - b) / (1 + math.exp(-theta[2])) * 0.99946 T, N = Z.shape47 Q = S.copy(); ll = 0.048 Qbar = (1 - a - b) * S - (g * Nbar if asym else 0.0)49 try:50 for t in range(T):51 if t > 0:52 z = Z[t - 1]53 Q = Qbar + a * np.outer(z, z) + b * Q54 if asym:55 n = z * (z < 0)56 Q = Q + g * np.outer(n, n)57 Rt = _corr_from_cov(Q)58 L = np.linalg.cholesky(Rt)59 w = np.linalg.solve(L, Z[t])60 ll += -0.5 * (2 * np.sum(np.log(np.diag(L))) + w @ w - Z[t] @ Z[t])61 except np.linalg.LinAlgError:62 return 1e1063 return -ll646566class CovModel:67 def __init__(self, model_id: str):68 self.model_id = model_id6970 def fit_forecast(self, R: np.ndarray) -> dict:71 """R: T x N returns (window). Returns dict(H1 = one-step covariance forecast, params, converged, time_s, min_eig)."""72 t0 = time.perf_counter()73 T, N = R.shape74 Rc = R - R.mean(axis=0)75 mid = self.model_id76 params = {}; conv = True77 if mid == "rolling_cov252":78 H1 = np.cov(Rc[-252:].T, bias=True)79 elif mid == "ewma_cov094":80 lam = 0.94; H = np.cov(Rc[:50].T, bias=True)81 for t in range(T):82 H = lam * H + (1 - lam) * np.outer(Rc[t], Rc[t])83 H1 = H84 elif mid in ("ccc_garch11", "dcc11", "adcc11"):85 S2, f1, fits = _univariate_garch(R)86 Z = Rc / np.sqrt(S2)87 S = np.corrcoef(Z.T)88 D1 = np.sqrt(f1)89 if mid == "ccc_garch11":90 R1 = S91 else:92 asym = mid == "adcc11"93 Nn = Z * (Z < 0); Nbar = (Nn.T @ Nn) / T94 th0 = np.array([math.log(0.02 / 0.98), math.log(0.95 / 0.05)] + ([math.log(0.02 / 0.98)] if asym else []))95 res = optimize.minimize(_dcc_negll, th0, args=(Z, S, asym, Nbar), method="L-BFGS-B", options={"maxiter": 500})96 conv = bool(res.success) and res.fun < 1e997 a = 1 / (1 + math.exp(-res.x[0])); b = (1 - a) / (1 + math.exp(-res.x[1])); g = (1 - a - b) / (1 + math.exp(-res.x[2])) * 0.999 if asym else 0.098 params = {"a": a, "b": b, "g": g, "loglik": -float(res.fun)}99 Q = S.copy(); Qbar = (1 - a - b) * S - g * Nbar100 for t in range(1, T + 1):101 z = Z[t - 1]102 Q = Qbar + a * np.outer(z, z) + b * Q103 if asym:104 n = z * (z < 0); Q = Q + g * np.outer(n, n)105 R1 = _corr_from_cov(Q)106 H1 = R1 * np.outer(D1, D1)107 params.update({f"garch_{j}": list(map(float, f.params)) for j, f in enumerate(fits)})108 elif mid == "scalar_bekk11":109 Sig = np.cov(Rc.T, bias=True)110111 def negll(th):112 a = 1 / (1 + math.exp(-th[0])); b = (1 - a) / (1 + math.exp(-th[1]))113 C = (1 - a - b) * Sig; H = Sig.copy(); ll = 0.0114 try:115 for t in range(T):116 if t > 0:117 H = C + a * np.outer(Rc[t - 1], Rc[t - 1]) + b * H118 L = np.linalg.cholesky(H); w = np.linalg.solve(L, Rc[t])119 ll += -0.5 * (2 * np.sum(np.log(np.diag(L))) + w @ w)120 except np.linalg.LinAlgError:121 return 1e10122 return -ll123 res = optimize.minimize(negll, np.array([math.log(0.05 / 0.95), math.log(0.9 / 0.1)]), method="L-BFGS-B", options={"maxiter": 500})124 conv = bool(res.success) and res.fun < 1e9125 a = 1 / (1 + math.exp(-res.x[0])); b = (1 - a) / (1 + math.exp(-res.x[1]))126 C = (1 - a - b) * Sig; H = Sig.copy()127 for t in range(1, T + 1):128 H = C + a * np.outer(Rc[t - 1], Rc[t - 1]) + b * H129 H1 = H; params = {"a": a, "b": b, "loglik": -float(res.fun)}130 else:131 raise KeyError(mid)132 H1 = 0.5 * (H1 + H1.T)133 return {"H1": H1, "params": params, "converged": conv, "time_s": time.perf_counter() - t0,134 "min_eig": float(np.linalg.eigvalsh(H1).min())}135136137def frobenius_loss(Y: np.ndarray, H: np.ndarray) -> float:138 return float(np.sum((Y - H) ** 2))139140141def mv_qlike(Y: np.ndarray, H: np.ndarray) -> float:142 """Multivariate QLIKE (Laurent, Rombouts & Violante 2013): log|H| + tr(H^{-1} Y) - log|Y| - N (>= 0 iff Y, H SPD)."""143 N = Y.shape[0]144 Hi = np.linalg.inv(H)145 sy, ldy = np.linalg.slogdet(Y); sh, ldh = np.linalg.slogdet(H)146 if sy <= 0 or sh <= 0:147 return float("nan")148 return float(ldh + np.trace(Hi @ Y) - ldy - N)149150151def portfolio_variance_loss(Y: np.ndarray, H: np.ndarray, w: np.ndarray) -> float:152 """QLIKE of the equally weighted (or given w) portfolio variance."""153 y = float(w @ Y @ w); h = float(w @ H @ w)154 return y / h - math.log(y / h) - 1.0155