mfe.multivariate¶

mfe.multivariate ¶

mfe.multivariate — Multivariate volatility models.

All of these are missing from the arch package (as of 2026).

Models

DCC Dynamic Conditional Correlation (Engle 2002); variants: dcc, cdcc, deco CCC Constant Conditional Correlation (Bollerslev 1990) BEKK BEKK model (Engle & Kroner 1995); variants: scalar, diagonal OGARCH Orthogonal GARCH (Alexander 2001) GOGARCH Generalized Orthogonal GARCH (van der Weide 2002); rotations: ica, moments RCC Rotated Conditional Correlation (Noureldin, Shephard & Sheppard 2014)

DCC ¶

DCC(variant: str = 'dcc')

DCC-GARCH(1,1) model (Engle 2002).

Parameters:

Name	Type	Description	Default
`variant`	`'dcc'(default) \| 'cdcc' \| 'deco'`		`'dcc'`

Source code in src/mfe/multivariate/dcc.py

def __init__(self, variant: str = "dcc") -> None:
    if variant not in ("dcc", "cdcc", "deco"):
        raise ValueError(f"variant must be 'dcc', 'cdcc', or 'deco', got '{variant}'")
    self.variant = variant

fit ¶

fit(data: FloatArray, starting_values: FloatArray | None = None) -> MultivariateVolResult

Two-step DCC estimation.

Parameters:

Name	Type	Description	Default
`data`	`(T, K) return matrix`		required
`starting_values`	`(2,) array [a, b]; if None uses [0.05, 0.90]`		`None`

Source code in src/mfe/multivariate/dcc.py

def fit(
    self,
    data: FloatArray,
    starting_values: FloatArray | None = None,
) -> MultivariateVolResult:
    """
    Two-step DCC estimation.

    Parameters
    ----------
    data : (T, K) return matrix
    starting_values : (2,) array [a, b]; if None uses [0.05, 0.90]
    """
    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape

    # Step 1: univariate GARCH
    h, z = _fit_univariate_garch(data)

    # Q_bar: sample covariance of standardized residuals
    q_bar = z.T @ z / T

    # Step 2: optimize DCC params
    x0 = np.array([0.05, 0.90]) if starting_values is None else np.asarray(starting_values)
    bounds = [(1e-6, 0.9999), (1e-6, 0.9999)]

    result = minimize(
        _dcc_log_likelihood,
        x0,
        args=(z, q_bar, h, data),
        method="L-BFGS-B",
        bounds=bounds,
        options={"maxiter": 500, "ftol": 1e-10},
    )

    if not result.success:
        warnings.warn(
            f"DCC did not converge: {result.message}",
            ConvergenceWarning,
            stacklevel=2,
        )

    a, b = float(result.x[0]), float(result.x[1])

    # Final covariances: Sigma_t = D_t R_t D_t
    Q = _dcc_recursion_numpy(z, q_bar, a, b)
    R = _q_to_correlation(Q)
    D = np.sqrt(h)  # (T, K) std dev
    sigma_t = np.empty((T, K, K), dtype=np.float64)
    for t in range(T):
        d = np.diag(D[t])
        sigma_t[t] = d @ R[t] @ d

    return MultivariateVolResult(
        params=result.x,
        log_likelihood=-result.fun,
        conditional_covariances=sigma_t,
        residuals=z,
        vcv=np.full((2, 2), np.nan),      # TODO: numeric Hessian
        vcv_robust=np.full((2, 2), np.nan),
        scores=np.zeros((T, 2)),
        converged=result.success,
        n_obs=T,
        n_params=2 + 3 * K,  # DCC + K GARCH triplets (omega, alpha, beta)
        model_name=f"DCC-GARCH(1,1) [{self.variant}]",
        diagnostics={"a": a, "b": b, "q_bar": q_bar},
    )

CCC ¶

CCC-GARCH(1,1) model (Bollerslev 1990).

No free parameters beyond the K univariate GARCH models.

fit ¶

fit(data: FloatArray) -> MultivariateVolResult

Estimate CCC-GARCH.

Parameters:

Name	Type	Description	Default
`data`	`(T, K) return matrix`		required

Source code in src/mfe/multivariate/ccc.py

def fit(self, data: FloatArray) -> MultivariateVolResult:
    """
    Estimate CCC-GARCH.

    Parameters
    ----------
    data : (T, K) return matrix
    """
    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape

    # Step 1: univariate GARCH
    h, z = _fit_univariate_garch(data)

    # Step 2: unconditional correlation of standardized residuals
    R = np.corrcoef(z.T)  # (K, K)

    # Sigma_t = D_t R D_t
    D = np.sqrt(h)  # (T, K)
    sigma_t = np.empty((T, K, K), dtype=np.float64)
    for t in range(T):
        d = np.diag(D[t])
        sigma_t[t] = d @ R @ d

    # Log-likelihood
    ll = 0.0
    for t in range(T):
        sign, logdet = np.linalg.slogdet(sigma_t[t])
        if sign <= 0:
            ll = -1e10
            break
        H_inv = np.linalg.inv(sigma_t[t])
        ll += -0.5 * (K * np.log(2 * np.pi) + logdet + float(data[t] @ H_inv @ data[t]))

    P = 3 * K  # omega, alpha, beta per asset; R is not a free param in estimation
    return MultivariateVolResult(
        params=np.array([]),  # R embedded in diagnostics
        log_likelihood=ll,
        conditional_covariances=sigma_t,
        residuals=z,
        vcv=np.zeros((0, 0)),
        vcv_robust=np.zeros((0, 0)),
        scores=np.zeros((T, 0)),
        converged=True,
        n_obs=T,
        n_params=P,
        model_name=f"CCC-GARCH(1,1), K={K}",
        diagnostics={"R": R},
    )

BEKK ¶

BEKK(variant: str | BEKKVariant = BEKKVariant.SCALAR)

BEKK-GARCH model (Engle & Kroner 1995).

Parameters:

Name	Type	Description	Default
`variant`	`'scalar' \| 'diagonal' \| 'full'`	"full" is not yet implemented.	`SCALAR`

Source code in src/mfe/multivariate/bekk.py

def __init__(self, variant: str | BEKKVariant = BEKKVariant.SCALAR) -> None:
    self.variant = BEKKVariant(variant)
    if self.variant == BEKKVariant.FULL:
        raise NotImplementedError(
            "Full BEKK is on the Phase 2 roadmap. "
            "Use variant='scalar' or 'diagonal'."
        )

fit ¶

fit(data: FloatArray, starting_values: FloatArray | None = None, method: str = 'L-BFGS-B', options: dict | None = None) -> MultivariateVolResult

Estimate BEKK parameters via QMLE.

Parameters:

Name	Type	Description	Default
`data`	`(T, K) return matrix (demeaned)`		required

Source code in src/mfe/multivariate/bekk.py

def fit(
    self,
    data: FloatArray,
    starting_values: FloatArray | None = None,
    method: str = "L-BFGS-B",
    options: dict | None = None,
) -> MultivariateVolResult:
    """
    Estimate BEKK parameters via QMLE.

    Parameters
    ----------
    data : (T, K) return matrix (demeaned)
    """
    eps = np.asarray(data, dtype=np.float64)
    T, K = eps.shape
    H0 = _backcast(eps)

    if self.variant == BEKKVariant.SCALAR:
        return self._fit_scalar(eps, K, T, H0, starting_values, method, options)
    else:
        return self._fit_diagonal(eps, K, T, H0, starting_values, method, options)

GOGARCH ¶

GOGARCH(n_components: int | None = None, rotation: Literal['ica', 'moments'] = 'ica')

GO-GARCH (Generalized Orthogonal GARCH) — van der Weide (2002).

Extends O-GARCH by estimating an additional orthogonal rotation U from fourth-order cumulants (ICA), so that the latent factors are as close to independent as possible.

Parameters:

Name	Type	Description	Default
`n_components`	`PCA components to retain (default: all K)`		`None`
`rotation`	`Literal['ica', 'moments']`	"ica" — FastICA deflationary algorithm, kurtosis contrast "moments" — Boswijk & van der Weide (2011) cumulant minimisation	`'ica'`

Source code in src/mfe/multivariate/gogarch.py

def __init__(
    self,
    n_components: int | None = None,
    rotation: Literal["ica", "moments"] = "ica",
) -> None:
    if rotation not in ("ica", "moments"):
        raise ValueError(f"rotation must be 'ica' or 'moments', got '{rotation}'")
    self.n_components = n_components
    self.rotation = rotation

fit ¶

fit(data: FloatArray) -> GOGARCHResult

Estimate GO-GARCH.

Parameters:

Name	Type	Description	Default
`data`	`(T, K) demeaned return matrix`		required

Source code in src/mfe/multivariate/gogarch.py

def fit(self, data: FloatArray) -> GOGARCHResult:
    """
    Estimate GO-GARCH.

    Parameters
    ----------
    data : (T, K) demeaned return matrix
    """
    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape
    K_c = K if self.n_components is None else min(self.n_components, K)

    # Step 1: PCA whitening
    factors_white, W_pca, eigenvals, _ = _pca_whiten(data, n_components=K_c)

    # Step 2: estimate orthogonal rotation U from cumulants
    if self.rotation == "ica":
        U = _fastica_rotation(factors_white)
    else:
        U = _moments_rotation(factors_white)

    # Latent factors: f_t = U @ f_white_t  (rows of U = unmixing directions)
    # Convention: factors = f_white @ U.T  so that factor[t] = U @ f_white[t]
    factors = factors_white @ U.T    # (T, K_c)

    # Full mixing matrix: data_centered ≈ factors @ W'
    # W = W_pca @ U.T  (K, K_c)
    W = W_pca @ U.T

    # Step 3: fit GARCH(1,1) to each latent factor
    h_factors, garch_res = _fit_factor_garch(factors)

    # Step 4: assemble Sigma_t = W H_t W'
    sigma_t = _assemble_covariances(W, h_factors)

    ll = _gogarch_loglik(sigma_t, data)
    converged = all(r.convergence_flag == 0 for r in garch_res)

    return GOGARCHResult(
        log_likelihood=ll,
        conditional_covariances=sigma_t,
        factors=factors,
        factor_variances=h_factors,
        mixing_matrix=W,
        rotation_matrix=U,
        eigenvalues=eigenvals,
        converged=converged,
        n_obs=T,
        n_components=K_c,
        model_name=f"GO-GARCH [{self.rotation}], K={K}, K_c={K_c}",
        garch_results=garch_res,
        diagnostics={
            "eigenvalues": eigenvals,
            "W_pca": W_pca,
            "U": U,
            "rotation": self.rotation,
        },
    )

OGARCH ¶

OGARCH(n_components: int | None = None)

O-GARCH (Orthogonal GARCH) — Alexander (2001).

Factor loadings are fixed at PCA eigenvectors. Each factor follows an independent GARCH(1,1). No rotation optimisation.

Parameters:

Name	Type	Description	Default
`n_components`	`number of PCA factors to retain (default: all K)`		`None`

Source code in src/mfe/multivariate/gogarch.py

def __init__(self, n_components: int | None = None) -> None:
    self.n_components = n_components

fit ¶

fit(data: FloatArray) -> GOGARCHResult

Estimate O-GARCH.

Parameters:

Name	Type	Description	Default
`data`	`(T, K) demeaned return matrix`		required

Source code in src/mfe/multivariate/gogarch.py

def fit(self, data: FloatArray) -> GOGARCHResult:
    """
    Estimate O-GARCH.

    Parameters
    ----------
    data : (T, K) demeaned return matrix
    """
    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape
    K_c = K if self.n_components is None else min(self.n_components, K)

    # Step 1: PCA whitening
    factors, W_pca, eigenvals, _ = _pca_whiten(data, n_components=K_c)

    # Step 2: fit GARCH(1,1) to each factor
    h_factors, garch_res = _fit_factor_garch(factors)

    # Step 3: assemble Sigma_t = W_pca H_t W_pca'
    sigma_t = _assemble_covariances(W_pca, h_factors)

    ll = _gogarch_loglik(sigma_t, data)
    converged = all(r.convergence_flag == 0 for r in garch_res)

    return GOGARCHResult(
        log_likelihood=ll,
        conditional_covariances=sigma_t,
        factors=factors,
        factor_variances=h_factors,
        mixing_matrix=W_pca,
        rotation_matrix=None,
        eigenvalues=eigenvals,
        converged=converged,
        n_obs=T,
        n_components=K_c,
        model_name=f"O-GARCH, K={K}, K_c={K_c}",
        garch_results=garch_res,
        diagnostics={"eigenvalues": eigenvals},
    )

GOGARCHResult `dataclass` ¶

GOGARCHResult(log_likelihood: float, conditional_covariances: FloatArray, factors: FloatArray, factor_variances: FloatArray, mixing_matrix: FloatArray, rotation_matrix: FloatArray | None, eigenvalues: FloatArray, converged: bool, n_obs: int, n_components: int, model_name: str, garch_results: list = list(), diagnostics: dict = dict())

Extended result for GO-GARCH / O-GARCH models.

factor_correlations ¶

factor_correlations() -> FloatArray

(K, K) unconditional correlation of the K latent factors with the original returns (mixing matrix scaled to unit-variance factors).

Source code in src/mfe/multivariate/gogarch.py

def factor_correlations(self) -> FloatArray:
    """
    (K, K) unconditional correlation of the K latent factors with the
    original returns (mixing matrix scaled to unit-variance factors).
    """
    return self.mixing_matrix / np.sqrt(self.eigenvalues)[None, :]

RCC ¶

RCC(rotation: str = 'symmetric')

Rotated Conditional Correlation (RCC) model.

Noureldin, Shephard & Sheppard (2014). Scalar parameterisation only (full RARCH — fully parametric A, B matrices — is left as a future extension).

Parameters:

Name	Type	Description	Default
`rotation`	`'symmetric'(default) \| 'cholesky'`	How to compute P^{1/2}: "symmetric" — symmetric (spectral) square root. PSD-preserving, recommended. "cholesky" — lower-triangular Cholesky. Faster but ordering-dependent.	`'symmetric'`

Source code in src/mfe/multivariate/rcc.py

def __init__(self, rotation: str = "symmetric") -> None:
    if rotation not in ("symmetric", "cholesky"):
        raise ValueError(f"rotation must be 'symmetric' or 'cholesky', got '{rotation}'")
    self.rotation = rotation

fit ¶

fit(data: FloatArray, starting_values: FloatArray | None = None, options: dict | None = None) -> RCCResult

Estimate RCC by two-step QML.

Parameters:

Name	Type	Default
`data`	`FloatArray`	required
`starting_values`	`FloatArray \| None`	`None`
`options`	`dict \| None`	`None`

Returns:

Type	Description
`RCCResult`

Source code in src/mfe/multivariate/rcc.py

def fit(
    self,
    data: FloatArray,
    starting_values: FloatArray | None = None,
    options: dict | None = None,
) -> RCCResult:
    """
    Estimate RCC by two-step QML.

    Parameters
    ----------
    data             : (T, K) return matrix (demeaned)
    starting_values  : (2,) [a0, b0]; if None uses [0.05, 0.90]
    options          : passed to scipy.optimize.minimize

    Returns
    -------
    RCCResult
    """
    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape

    # Step 1: Rotation
    P = data.T @ data / T  # unconditional covariance

    if self.rotation == "symmetric":
        P_half, P_inv_half = _rotation_matrix(P)
    else:
        try:
            L = np.linalg.cholesky(P)
            P_half = L
            P_inv_half = np.linalg.inv(L)
        except np.linalg.LinAlgError:
            P_half, P_inv_half = _rotation_matrix(P)

    u = data @ P_inv_half.T   # (T, K) rotated residuals; u_t = P^{-T/2} r_t
    # Verify: u.T @ u / T ≈ I_K
    # (P_inv_half.T @ P @ P_inv_half = I if symmetric, approx otherwise)

    # Step 2: Optimize (a, b) in rotated space
    x0 = np.array([0.05, 0.90]) if starting_values is None else np.asarray(starting_values)
    bounds = [(1e-6, 0.9999), (1e-6, 0.9999)]

    result = minimize(
        _rcc_loglik,
        x0,
        args=(u,),
        method="L-BFGS-B",
        bounds=bounds,
        options=options or {"maxiter": 500, "ftol": 1e-10},
    )

    if not result.success:
        warnings.warn(
            f"RCC did not converge: {result.message}",
            ConvergenceWarning,
            stacklevel=2,
        )

    a, b = float(result.x[0]), float(result.x[1])

    # Final G_t in rotated space
    if _HAS_CYTHON:
        G_t = np.asarray(_dcc_q_recursion(
            np.ascontiguousarray(u, dtype=np.float64),
            np.ascontiguousarray(np.eye(K), dtype=np.float64),
            a, b,
        ))
    else:
        G_t = _rcc_recursion_numpy(u, a, b)

    # Reconstruct Sigma_t = P^{1/2} G_t P^{1/2}
    sigma_t = np.empty((T, K, K), dtype=np.float64)
    for t in range(T):
        sigma_t[t] = P_half @ G_t[t] @ P_half.T

    # Log-likelihood (step 1 + step 2)
    # Step 1: Gaussian log-lik with const cov P
    sign, logdet_P = np.linalg.slogdet(P)
    ll_step1 = -0.5 * T * (K * np.log(2 * np.pi) + logdet_P) - 0.5 * float(np.sum(u ** 2))
    ll_step2 = -result.fun  # stored as negative by _rcc_loglik
    ll_total = ll_step1 + ll_step2

    return RCCResult(
        params=result.x,
        log_likelihood=ll_total,
        conditional_covariances=sigma_t,
        G_t=G_t,
        u_t=u,
        P=P,
        P_half=P_half,
        converged=result.success,
        n_obs=T,
        n_vars=K,
        diagnostics={"a": a, "b": b, "rotation": self.rotation},
    )

RCCResult `dataclass` ¶

RCCResult(params: FloatArray, log_likelihood: float, conditional_covariances: FloatArray, G_t: FloatArray, u_t: FloatArray, P: FloatArray, P_half: FloatArray, converged: bool, n_obs: int, n_vars: int, diagnostics: dict = dict())

RCC model estimation result.

conditional_correlations ¶

conditional_correlations() -> FloatArray

Extract (T, K, K) conditional correlation matrices from Sigma_t.

Source code in src/mfe/multivariate/rcc.py

def conditional_correlations(self) -> FloatArray:
    """
    Extract (T, K, K) conditional correlation matrices from Sigma_t.
    """
    T, K, _ = self.conditional_covariances.shape
    R = np.empty_like(self.conditional_covariances)
    for t in range(T):
        S = self.conditional_covariances[t]
        d = np.sqrt(np.diag(S))
        d_inv = np.where(d > 0, 1.0 / d, 0.0)
        R[t] = d_inv[:, None] * S * d_inv[None, :]
    return R

base ¶

Abstract base class for multivariate volatility models.

Design principles: - Results are dataclasses, not mutable objects, to avoid state mutation bugs. - Convergence failures surface as ConvergenceWarning, not silent bad params. - Robust VCV (sandwich) is computed by default alongside the Hessian-only VCV.

ConvergenceWarning ¶

Bases: UserWarning

Raised when the optimizer did not fully converge.

MultivariateVolResult `dataclass` ¶

MultivariateVolResult(params: FloatArray, log_likelihood: float, conditional_covariances: FloatArray, residuals: FloatArray, vcv: FloatArray, vcv_robust: FloatArray, scores: FloatArray, converged: bool = True, n_obs: int = 0, n_params: int = 0, model_name: str = '', diagnostics: dict = dict())

Container for multivariate volatility estimation results.

MultivariateVolatilityProcess ¶

Bases: ABC

Base class for all multivariate volatility models.

Subclasses must implement: _log_likelihood(params, data) -> float _compute_covariances(params, data) -> (T, K, K) _starting_values(data) -> FloatArray _parameter_bounds(data, K) -> list[tuple[float, float]] _parameter_names(K) -> list[str]

fit ¶

fit(data: FloatArray, starting_values: FloatArray | None = None, method: str = 'L-BFGS-B', options: dict | None = None) -> MultivariateVolResult

Fit the model via quasi-maximum likelihood.

Parameters:

Name	Type	Default
`data`	`FloatArray`	required
`starting_values`	`(P,) starting parameter vector; if None, uses heuristic`	`None`
`method`	`str`	`'L-BFGS-B'`
`options`	`dict \| None`	`None`

Returns:

Type	Description
`MultivariateVolResult`

Source code in src/mfe/multivariate/base.py

def fit(
    self,
    data: FloatArray,
    starting_values: FloatArray | None = None,
    method: str = "L-BFGS-B",
    options: dict | None = None,
) -> MultivariateVolResult:
    """
    Fit the model via quasi-maximum likelihood.

    Parameters
    ----------
    data            : (T, K) return matrix
    starting_values : (P,) starting parameter vector; if None, uses heuristic
    method          : scipy.optimize.minimize method
    options         : passed to minimize

    Returns
    -------
    MultivariateVolResult
    """
    from scipy.optimize import minimize

    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape

    if starting_values is None:
        x0 = self._starting_values(data)
    else:
        x0 = np.asarray(starting_values, dtype=np.float64)

    bounds = self._parameter_bounds(data, K)

    result = minimize(
        self._log_likelihood,
        x0,
        args=(data,),
        method=method,
        bounds=bounds,
        options=options or {"maxiter": 1000, "ftol": 1e-9},
    )

    if not result.success:
        warnings.warn(
            f"{self.__class__.__name__} did not converge: {result.message}. "
            "Use results with caution.",
            ConvergenceWarning,
            stacklevel=2,
        )

    params = result.x
    sigma_t = self._compute_covariances(params, data)

    # Score matrix via finite differences (TODO: analytic scores per model)
    eps = 1e-6
    P = len(params)
    scores = np.zeros((T, P), dtype=np.float64)
    # ... (placeholder: per-obs score via finite differences)

    # Hessian from scipy (numerical)
    try:
        from scipy.optimize import approx_fprime
        hessian = np.zeros((P, P), dtype=np.float64)
        for i in range(P):
            ei = np.zeros(P)
            ei[i] = eps
            g_plus = approx_fprime(params + ei, self._log_likelihood, eps, data)
            g_minus = approx_fprime(params - ei, self._log_likelihood, eps, data)
            hessian[i] = (g_plus - g_minus) / (2 * eps)
        hessian = 0.5 * (hessian + hessian.T)
        vcv = np.linalg.inv(hessian) if np.linalg.det(hessian) != 0 else np.full((P, P), np.nan)
    except Exception:
        vcv = np.full((P, P), np.nan)

    vcv_robust = sandwich(scores, hessian) if not np.any(np.isnan(scores)) else vcv

    return MultivariateVolResult(
        params=params,
        log_likelihood=-result.fun,
        conditional_covariances=sigma_t,
        residuals=data,  # caller should standardize if needed
        vcv=vcv,
        vcv_robust=vcv_robust,
        scores=scores,
        converged=result.success,
        n_obs=T,
        n_params=P,
        model_name=self.__class__.__name__,
        diagnostics={"optimizer_result": result},
    )

bekk ¶

BEKK-GARCH models (Engle & Kroner 1995).

Three variants: Scalar BEKK: H_t = C'C + a^2 * eps_{t-1} eps_{t-1}' + b^2 * H_{t-1} Diagonal BEKK: H_t = C'C + A' * eps eps' * A + B' * H_{t-1} * B (A,B diagonal) Full BEKK: H_t = C'C + A' * eps eps' * A + B' * H_{t-1} * B (A,B full K×K)

Estimation via QMLE (two-step or direct). We implement the scalar and diagonal variants first; full BEKK is numerically expensive (O(K^4) per recursion step) and relegated to Phase 2.

Key performance note: - Scalar BEKK inner loop: O(K^2 * T) — fast enough in numpy for K <= 20 - Diagonal BEKK: O(K^2 * T) — same - Full BEKK: O(K^4 * T) — needs Cython for K > 5

References

Engle, R.F. & Kroner, K.F. (1995): "Multivariate Simultaneous Generalized ARCH", Econometric Theory.

Noureldin, D., Shephard, N. & Sheppard, K. (2012): "Multivariate High-Frequency-Based Volatility (HEAVY) Models", JoE.

BEKK ¶

BEKK(variant: str | BEKKVariant = BEKKVariant.SCALAR)

BEKK-GARCH model (Engle & Kroner 1995).

Parameters:

Name	Type	Description	Default
`variant`	`'scalar' \| 'diagonal' \| 'full'`	"full" is not yet implemented.	`SCALAR`

Source code in src/mfe/multivariate/bekk.py

def __init__(self, variant: str | BEKKVariant = BEKKVariant.SCALAR) -> None:
    self.variant = BEKKVariant(variant)
    if self.variant == BEKKVariant.FULL:
        raise NotImplementedError(
            "Full BEKK is on the Phase 2 roadmap. "
            "Use variant='scalar' or 'diagonal'."
        )

fit ¶

fit(data: FloatArray, starting_values: FloatArray | None = None, method: str = 'L-BFGS-B', options: dict | None = None) -> MultivariateVolResult

Estimate BEKK parameters via QMLE.

Parameters:

Name	Type	Description	Default
`data`	`(T, K) return matrix (demeaned)`		required

Source code in src/mfe/multivariate/bekk.py

def fit(
    self,
    data: FloatArray,
    starting_values: FloatArray | None = None,
    method: str = "L-BFGS-B",
    options: dict | None = None,
) -> MultivariateVolResult:
    """
    Estimate BEKK parameters via QMLE.

    Parameters
    ----------
    data : (T, K) return matrix (demeaned)
    """
    eps = np.asarray(data, dtype=np.float64)
    T, K = eps.shape
    H0 = _backcast(eps)

    if self.variant == BEKKVariant.SCALAR:
        return self._fit_scalar(eps, K, T, H0, starting_values, method, options)
    else:
        return self._fit_diagonal(eps, K, T, H0, starting_values, method, options)

ccc ¶

CCC-GARCH (Constant Conditional Correlation).

Bollerslev, T. (1990): "Modelling the Coherence in Short-Run Nominal Exchange Rates: A Multivariate Generalized ARCH Model", Review of Economics and Statistics.

H_t = D_t R D_t

Where D_t = diag(sigma_{1,t}, ..., sigma_{K,t}) from K independent GARCH(1,1) and R = unconditional correlation matrix (constant).

Two-step estimation: Step 1: Fit GARCH(1,1) to each series independently. Step 2: Compute R = sample correlation of standardized residuals.

This is just DCC with a=b=0 (zero dynamics in the correlation), but it's worth having as an explicit model for testing (CCC vs DCC likelihood ratio test).

CCC ¶

CCC-GARCH(1,1) model (Bollerslev 1990).

No free parameters beyond the K univariate GARCH models.

fit ¶

fit(data: FloatArray) -> MultivariateVolResult

Estimate CCC-GARCH.

Parameters:

Name	Type	Description	Default
`data`	`(T, K) return matrix`		required

Source code in src/mfe/multivariate/ccc.py

def fit(self, data: FloatArray) -> MultivariateVolResult:
    """
    Estimate CCC-GARCH.

    Parameters
    ----------
    data : (T, K) return matrix
    """
    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape

    # Step 1: univariate GARCH
    h, z = _fit_univariate_garch(data)

    # Step 2: unconditional correlation of standardized residuals
    R = np.corrcoef(z.T)  # (K, K)

    # Sigma_t = D_t R D_t
    D = np.sqrt(h)  # (T, K)
    sigma_t = np.empty((T, K, K), dtype=np.float64)
    for t in range(T):
        d = np.diag(D[t])
        sigma_t[t] = d @ R @ d

    # Log-likelihood
    ll = 0.0
    for t in range(T):
        sign, logdet = np.linalg.slogdet(sigma_t[t])
        if sign <= 0:
            ll = -1e10
            break
        H_inv = np.linalg.inv(sigma_t[t])
        ll += -0.5 * (K * np.log(2 * np.pi) + logdet + float(data[t] @ H_inv @ data[t]))

    P = 3 * K  # omega, alpha, beta per asset; R is not a free param in estimation
    return MultivariateVolResult(
        params=np.array([]),  # R embedded in diagnostics
        log_likelihood=ll,
        conditional_covariances=sigma_t,
        residuals=z,
        vcv=np.zeros((0, 0)),
        vcv_robust=np.zeros((0, 0)),
        scores=np.zeros((T, 0)),
        converged=True,
        n_obs=T,
        n_params=P,
        model_name=f"CCC-GARCH(1,1), K={K}",
        diagnostics={"R": R},
    )

dcc ¶

DCC-GARCH (Dynamic Conditional Correlation) model.

Engle, R. (2002): "Dynamic Conditional Correlations — A Simple Class of Multivariate Generalized Autoregressive Conditional Heteroskedasticity Models", JBES.

Two-step estimation: Step 1: Fit univariate GARCH(1,1) to each asset. Extract standardized residuals. Step 2: Estimate DCC parameters (a, b) by maximizing the correlation likelihood.

Also implements: - cDCC (Aielli 2013): consistent DCC (avoids bias in Q_bar estimation) - DECO (Engle & Kelly 2012): equicorrelation restricted DCC

Key fix vs. the MATLAB mfe-toolbox: The MATLAB dcc.m computes Q_bar = mean(z_t z_t') once at the start and inside the log-likelihood. For long panels this is fine. We pre-compute it and pass it as a constant to the inner loop to avoid the recomputation on every likelihood call.

DCC ¶

DCC(variant: str = 'dcc')

DCC-GARCH(1,1) model (Engle 2002).

Parameters:

Name	Type	Description	Default
`variant`	`'dcc'(default) \| 'cdcc' \| 'deco'`		`'dcc'`

Source code in src/mfe/multivariate/dcc.py

def __init__(self, variant: str = "dcc") -> None:
    if variant not in ("dcc", "cdcc", "deco"):
        raise ValueError(f"variant must be 'dcc', 'cdcc', or 'deco', got '{variant}'")
    self.variant = variant

fit ¶

fit(data: FloatArray, starting_values: FloatArray | None = None) -> MultivariateVolResult

Two-step DCC estimation.

Parameters:

Name	Type	Description	Default
`data`	`(T, K) return matrix`		required
`starting_values`	`(2,) array [a, b]; if None uses [0.05, 0.90]`		`None`

Source code in src/mfe/multivariate/dcc.py

def fit(
    self,
    data: FloatArray,
    starting_values: FloatArray | None = None,
) -> MultivariateVolResult:
    """
    Two-step DCC estimation.

    Parameters
    ----------
    data : (T, K) return matrix
    starting_values : (2,) array [a, b]; if None uses [0.05, 0.90]
    """
    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape

    # Step 1: univariate GARCH
    h, z = _fit_univariate_garch(data)

    # Q_bar: sample covariance of standardized residuals
    q_bar = z.T @ z / T

    # Step 2: optimize DCC params
    x0 = np.array([0.05, 0.90]) if starting_values is None else np.asarray(starting_values)
    bounds = [(1e-6, 0.9999), (1e-6, 0.9999)]

    result = minimize(
        _dcc_log_likelihood,
        x0,
        args=(z, q_bar, h, data),
        method="L-BFGS-B",
        bounds=bounds,
        options={"maxiter": 500, "ftol": 1e-10},
    )

    if not result.success:
        warnings.warn(
            f"DCC did not converge: {result.message}",
            ConvergenceWarning,
            stacklevel=2,
        )

    a, b = float(result.x[0]), float(result.x[1])

    # Final covariances: Sigma_t = D_t R_t D_t
    Q = _dcc_recursion_numpy(z, q_bar, a, b)
    R = _q_to_correlation(Q)
    D = np.sqrt(h)  # (T, K) std dev
    sigma_t = np.empty((T, K, K), dtype=np.float64)
    for t in range(T):
        d = np.diag(D[t])
        sigma_t[t] = d @ R[t] @ d

    return MultivariateVolResult(
        params=result.x,
        log_likelihood=-result.fun,
        conditional_covariances=sigma_t,
        residuals=z,
        vcv=np.full((2, 2), np.nan),      # TODO: numeric Hessian
        vcv_robust=np.full((2, 2), np.nan),
        scores=np.zeros((T, 2)),
        converged=result.success,
        n_obs=T,
        n_params=2 + 3 * K,  # DCC + K GARCH triplets (omega, alpha, beta)
        model_name=f"DCC-GARCH(1,1) [{self.variant}]",
        diagnostics={"a": a, "b": b, "q_bar": q_bar},
    )

gogarch ¶

GO-GARCH: Generalized Orthogonal GARCH.

Two closely related models:

O-GARCH (Alexander 2001): Factors = PCA rotation of returns. Factor loadings fixed at eigenvectors. Each factor follows an independent GARCH(1,1). Sigma_t = W diag(h_{1,t}, ..., h_{K,t}) W' where W = eigenvectors of the unconditional covariance, scaled to unit-variance factors.

GO-GARCH (van der Weide 2002): Extends O-GARCH. The mixing matrix is W = P * U where: P = PCA whitening matrix (from unconditional covariance) U = K×K orthogonal matrix estimated from higher-order cumulants (ICA-style; we use the FastICA / JADE approach) Sigma_t = (PU) diag(h_{1,t}, ..., h_{K,t}) (PU)'

References

Alexander, C. (2001): "Orthogonal GARCH", in Mastering Risk, vol. 2, FT Prentice Hall.

van der Weide, R. (2002): "GO-GARCH: A Multivariate Generalized Orthogonal GARCH Model", Journal of Applied Econometrics, 17(5), 549-564.

Boswijk, H.P. & van der Weide, R. (2011): "Method of Moments Estimation of GO-GARCH Models", Journal of Econometrics, 163(1), 118-126.

Implementation notes

O-GARCH is exact PCA-GARCH: W is fixed from eigendecomposition, no additional optimisation.
GO-GARCH adds an orthogonal rotation U estimated via one of: "moments" — Boswijk & van der Weide (2011) GMM on fourth-order cumulants "ica" — FastICA (deflationary, kurtosis contrast) The MATLAB mfe-toolbox gogarch.m uses a direct numerical optimisation over U. That approach has a memory-leak issue (anonymous function closes over volData in a loop). We avoid this entirely by using the closed-form cumulant matching.
Factor GARCH: each factor uses GARCH(1,1) from the arch package.
The full conditional covariance is assembled from factor variances in O(K^2 * T).

MATLAB bugs fixed

Memory leak: the MATLAB version closes over volData in a nested fmincon call inside a loop. We avoid closures entirely — all state is explicit.
No convergence warning: MATLAB silently used non-converged parameters. We raise ConvergenceWarning on any non-converged factor GARCH.
Inconsistent factor ordering: MATLAB returns factors in eigenvalue order (descending). We follow the same convention but document it explicitly.

GOGARCHResult `dataclass` ¶

GOGARCHResult(log_likelihood: float, conditional_covariances: FloatArray, factors: FloatArray, factor_variances: FloatArray, mixing_matrix: FloatArray, rotation_matrix: FloatArray | None, eigenvalues: FloatArray, converged: bool, n_obs: int, n_components: int, model_name: str, garch_results: list = list(), diagnostics: dict = dict())

Extended result for GO-GARCH / O-GARCH models.

factor_correlations ¶

factor_correlations() -> FloatArray

(K, K) unconditional correlation of the K latent factors with the original returns (mixing matrix scaled to unit-variance factors).

Source code in src/mfe/multivariate/gogarch.py

def factor_correlations(self) -> FloatArray:
    """
    (K, K) unconditional correlation of the K latent factors with the
    original returns (mixing matrix scaled to unit-variance factors).
    """
    return self.mixing_matrix / np.sqrt(self.eigenvalues)[None, :]

OGARCH ¶

OGARCH(n_components: int | None = None)

O-GARCH (Orthogonal GARCH) — Alexander (2001).

Factor loadings are fixed at PCA eigenvectors. Each factor follows an independent GARCH(1,1). No rotation optimisation.

Parameters:

Name	Type	Description	Default
`n_components`	`number of PCA factors to retain (default: all K)`		`None`

Source code in src/mfe/multivariate/gogarch.py

def __init__(self, n_components: int | None = None) -> None:
    self.n_components = n_components

fit ¶

fit(data: FloatArray) -> GOGARCHResult

Estimate O-GARCH.

Parameters:

Name	Type	Description	Default
`data`	`(T, K) demeaned return matrix`		required

Source code in src/mfe/multivariate/gogarch.py

def fit(self, data: FloatArray) -> GOGARCHResult:
    """
    Estimate O-GARCH.

    Parameters
    ----------
    data : (T, K) demeaned return matrix
    """
    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape
    K_c = K if self.n_components is None else min(self.n_components, K)

    # Step 1: PCA whitening
    factors, W_pca, eigenvals, _ = _pca_whiten(data, n_components=K_c)

    # Step 2: fit GARCH(1,1) to each factor
    h_factors, garch_res = _fit_factor_garch(factors)

    # Step 3: assemble Sigma_t = W_pca H_t W_pca'
    sigma_t = _assemble_covariances(W_pca, h_factors)

    ll = _gogarch_loglik(sigma_t, data)
    converged = all(r.convergence_flag == 0 for r in garch_res)

    return GOGARCHResult(
        log_likelihood=ll,
        conditional_covariances=sigma_t,
        factors=factors,
        factor_variances=h_factors,
        mixing_matrix=W_pca,
        rotation_matrix=None,
        eigenvalues=eigenvals,
        converged=converged,
        n_obs=T,
        n_components=K_c,
        model_name=f"O-GARCH, K={K}, K_c={K_c}",
        garch_results=garch_res,
        diagnostics={"eigenvalues": eigenvals},
    )

GOGARCH ¶

GOGARCH(n_components: int | None = None, rotation: Literal['ica', 'moments'] = 'ica')

GO-GARCH (Generalized Orthogonal GARCH) — van der Weide (2002).

Extends O-GARCH by estimating an additional orthogonal rotation U from fourth-order cumulants (ICA), so that the latent factors are as close to independent as possible.

Parameters:

Name	Type	Description	Default
`n_components`	`PCA components to retain (default: all K)`		`None`
`rotation`	`Literal['ica', 'moments']`	"ica" — FastICA deflationary algorithm, kurtosis contrast "moments" — Boswijk & van der Weide (2011) cumulant minimisation	`'ica'`

Source code in src/mfe/multivariate/gogarch.py

def __init__(
    self,
    n_components: int | None = None,
    rotation: Literal["ica", "moments"] = "ica",
) -> None:
    if rotation not in ("ica", "moments"):
        raise ValueError(f"rotation must be 'ica' or 'moments', got '{rotation}'")
    self.n_components = n_components
    self.rotation = rotation

fit ¶

fit(data: FloatArray) -> GOGARCHResult

Estimate GO-GARCH.

Parameters:

Name	Type	Description	Default
`data`	`(T, K) demeaned return matrix`		required

Source code in src/mfe/multivariate/gogarch.py

def fit(self, data: FloatArray) -> GOGARCHResult:
    """
    Estimate GO-GARCH.

    Parameters
    ----------
    data : (T, K) demeaned return matrix
    """
    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape
    K_c = K if self.n_components is None else min(self.n_components, K)

    # Step 1: PCA whitening
    factors_white, W_pca, eigenvals, _ = _pca_whiten(data, n_components=K_c)

    # Step 2: estimate orthogonal rotation U from cumulants
    if self.rotation == "ica":
        U = _fastica_rotation(factors_white)
    else:
        U = _moments_rotation(factors_white)

    # Latent factors: f_t = U @ f_white_t  (rows of U = unmixing directions)
    # Convention: factors = f_white @ U.T  so that factor[t] = U @ f_white[t]
    factors = factors_white @ U.T    # (T, K_c)

    # Full mixing matrix: data_centered ≈ factors @ W'
    # W = W_pca @ U.T  (K, K_c)
    W = W_pca @ U.T

    # Step 3: fit GARCH(1,1) to each latent factor
    h_factors, garch_res = _fit_factor_garch(factors)

    # Step 4: assemble Sigma_t = W H_t W'
    sigma_t = _assemble_covariances(W, h_factors)

    ll = _gogarch_loglik(sigma_t, data)
    converged = all(r.convergence_flag == 0 for r in garch_res)

    return GOGARCHResult(
        log_likelihood=ll,
        conditional_covariances=sigma_t,
        factors=factors,
        factor_variances=h_factors,
        mixing_matrix=W,
        rotation_matrix=U,
        eigenvalues=eigenvals,
        converged=converged,
        n_obs=T,
        n_components=K_c,
        model_name=f"GO-GARCH [{self.rotation}], K={K}, K_c={K_c}",
        garch_results=garch_res,
        diagnostics={
            "eigenvalues": eigenvals,
            "W_pca": W_pca,
            "U": U,
            "rotation": self.rotation,
        },
    )

rcc ¶

Rotated Conditional Correlation (RCC) model.

Noureldin, D., Shephard, N. & Sheppard, K. (2014): "Multivariate Rotated ARCH Models", Journal of Econometrics, 179(1), 16-30.

Core idea

Standard DCC operates on the correlation structure of standardised returns. RCC instead:

Rotates raw returns by the inverse square-root of their unconditional covariance: u_t = P^{-1/2} r_t where P = E[r_t r_t'] Under the true DGP, u_t has unconditional covariance I_K.
Fits a BEKK-type process to u_t u_t' in the rotated space: G_t = (I_K - A - B) + A * u_{t-1} u_{t-1}' * A + B * G_{t-1} * B With covariance targeting, the unconditional mean of G_t is I_K by construction, which removes the free intercept.
The DCC-type RCC sets A = aI_K, B = bI_K (scalar): G_t = (1 - a - b) I_K + a u_{t-1}u_{t-1}' + b G_{t-1}
Reconstructs the conditional covariance of original returns: Sigma_t = P^{1/2} G_t P^{1/2}

Advantages over DCC

Covariance targeting is exact by construction (no approximate initialisation).
Estimation is more stable for large K: scalar RCC has only 2 free parameters (a, b) regardless of dimension, compared to DCC which also has 2 but its inner Q-bar computation can be numerically unstable.
The rotated space has a cleaner likelihood because the rotated residuals u_t have unconditional identity covariance, so the step-2 likelihood simplifies.

Two-step estimation

Step 1: Estimate P = sample covariance of returns. Compute P^{1/2} (Cholesky or spectral decomp). Rotate: u_t = P^{-1/2} r_t. Step 2: Estimate scalar (a, b) by maximising the correlation log-likelihood using the G_t recursion.

This matches the MFE MATLAB rcc.m implementation.

RCCResult `dataclass` ¶

RCCResult(params: FloatArray, log_likelihood: float, conditional_covariances: FloatArray, G_t: FloatArray, u_t: FloatArray, P: FloatArray, P_half: FloatArray, converged: bool, n_obs: int, n_vars: int, diagnostics: dict = dict())

RCC model estimation result.

conditional_correlations ¶

conditional_correlations() -> FloatArray

Extract (T, K, K) conditional correlation matrices from Sigma_t.

Source code in src/mfe/multivariate/rcc.py

def conditional_correlations(self) -> FloatArray:
    """
    Extract (T, K, K) conditional correlation matrices from Sigma_t.
    """
    T, K, _ = self.conditional_covariances.shape
    R = np.empty_like(self.conditional_covariances)
    for t in range(T):
        S = self.conditional_covariances[t]
        d = np.sqrt(np.diag(S))
        d_inv = np.where(d > 0, 1.0 / d, 0.0)
        R[t] = d_inv[:, None] * S * d_inv[None, :]
    return R

RCC ¶

RCC(rotation: str = 'symmetric')

Rotated Conditional Correlation (RCC) model.

Noureldin, Shephard & Sheppard (2014). Scalar parameterisation only (full RARCH — fully parametric A, B matrices — is left as a future extension).

Parameters:

Name	Type	Description	Default
`rotation`	`'symmetric'(default) \| 'cholesky'`	How to compute P^{1/2}: "symmetric" — symmetric (spectral) square root. PSD-preserving, recommended. "cholesky" — lower-triangular Cholesky. Faster but ordering-dependent.	`'symmetric'`

Source code in src/mfe/multivariate/rcc.py

def __init__(self, rotation: str = "symmetric") -> None:
    if rotation not in ("symmetric", "cholesky"):
        raise ValueError(f"rotation must be 'symmetric' or 'cholesky', got '{rotation}'")
    self.rotation = rotation

fit ¶

fit(data: FloatArray, starting_values: FloatArray | None = None, options: dict | None = None) -> RCCResult

Estimate RCC by two-step QML.

Parameters:

Name	Type	Default
`data`	`FloatArray`	required
`starting_values`	`FloatArray \| None`	`None`
`options`	`dict \| None`	`None`

Returns:

Type	Description
`RCCResult`

Source code in src/mfe/multivariate/rcc.py

def fit(
    self,
    data: FloatArray,
    starting_values: FloatArray | None = None,
    options: dict | None = None,
) -> RCCResult:
    """
    Estimate RCC by two-step QML.

    Parameters
    ----------
    data             : (T, K) return matrix (demeaned)
    starting_values  : (2,) [a0, b0]; if None uses [0.05, 0.90]
    options          : passed to scipy.optimize.minimize

    Returns
    -------
    RCCResult
    """
    data = np.asarray(data, dtype=np.float64)
    T, K = data.shape

    # Step 1: Rotation
    P = data.T @ data / T  # unconditional covariance

    if self.rotation == "symmetric":
        P_half, P_inv_half = _rotation_matrix(P)
    else:
        try:
            L = np.linalg.cholesky(P)
            P_half = L
            P_inv_half = np.linalg.inv(L)
        except np.linalg.LinAlgError:
            P_half, P_inv_half = _rotation_matrix(P)

    u = data @ P_inv_half.T   # (T, K) rotated residuals; u_t = P^{-T/2} r_t
    # Verify: u.T @ u / T ≈ I_K
    # (P_inv_half.T @ P @ P_inv_half = I if symmetric, approx otherwise)

    # Step 2: Optimize (a, b) in rotated space
    x0 = np.array([0.05, 0.90]) if starting_values is None else np.asarray(starting_values)
    bounds = [(1e-6, 0.9999), (1e-6, 0.9999)]

    result = minimize(
        _rcc_loglik,
        x0,
        args=(u,),
        method="L-BFGS-B",
        bounds=bounds,
        options=options or {"maxiter": 500, "ftol": 1e-10},
    )

    if not result.success:
        warnings.warn(
            f"RCC did not converge: {result.message}",
            ConvergenceWarning,
            stacklevel=2,
        )

    a, b = float(result.x[0]), float(result.x[1])

    # Final G_t in rotated space
    if _HAS_CYTHON:
        G_t = np.asarray(_dcc_q_recursion(
            np.ascontiguousarray(u, dtype=np.float64),
            np.ascontiguousarray(np.eye(K), dtype=np.float64),
            a, b,
        ))
    else:
        G_t = _rcc_recursion_numpy(u, a, b)

    # Reconstruct Sigma_t = P^{1/2} G_t P^{1/2}
    sigma_t = np.empty((T, K, K), dtype=np.float64)
    for t in range(T):
        sigma_t[t] = P_half @ G_t[t] @ P_half.T

    # Log-likelihood (step 1 + step 2)
    # Step 1: Gaussian log-lik with const cov P
    sign, logdet_P = np.linalg.slogdet(P)
    ll_step1 = -0.5 * T * (K * np.log(2 * np.pi) + logdet_P) - 0.5 * float(np.sum(u ** 2))
    ll_step2 = -result.fun  # stored as negative by _rcc_loglik
    ll_total = ll_step1 + ll_step2

    return RCCResult(
        params=result.x,
        log_likelihood=ll_total,
        conditional_covariances=sigma_t,
        G_t=G_t,
        u_t=u,
        P=P,
        P_half=P_half,
        converged=result.success,
        n_obs=T,
        n_vars=K,
        diagnostics={"a": a, "b": b, "rotation": self.rotation},
    )

mfe.multivariate¶