GARCH and BEKK-GARCH

Definition

GARCH (Generalized Autoregressive Conditional Heteroscedasticity) models the conditional variance of a time series as a function of past squared innovations and past conditional variances. The Baba-Engle-Kraft-Kroner (BEKK) parametrization (Engle and Kroner 1995) extends GARCH to multivariate systems while guaranteeing positive definiteness of the conditional covariance matrix at every time step.

Key Ideas

Conditional heteroscedasticity — the variance of the innovation at time $t$ , conditional on the past, is itself time-varying. Unconditional variance may still be finite (covariance-stationary GARCH) or infinite (Integrated GARCH, IGARCH).
GARCH variance clusters: large shocks tend to be followed by further large shocks, small shocks by small ones. This matches the empirical behavior of financial returns and interest rate innovations.
In multivariate systems the BEKK form ensures $H_t \succ 0$ without imposing diagonal structure, allowing shock spillovers and common volatility factors across equations.
QML (Quasi Maximum Likelihood) estimation maximizes the Gaussian log-likelihood even when innovations are non-Gaussian; the estimator is consistent and asymptotically normal under mild regularity conditions.

How It Works

Univariate ARCH and GARCH

ARCH( $q$ ) — Engle (1982): Let $\varepsilon_t = \sigma_t z_t$ with $z_t \overset{iid}{\sim}(0,1)$ . The conditional variance follows

$\sigma_t^2 = \omega + \sum_{i=1}^{q} \alpha_i \varepsilon_{t-i}^2, \qquad \omega > 0,\; \alpha_i \geq 0 \tag{ARCH}$

GARCH( $p,q$ ) — Bollerslev (1986): Augments ARCH with autoregressive variance terms:

$\sigma_t^2 = \omega + \sum_{i=1}^{q} \alpha_i \varepsilon_{t-i}^2 + \sum_{j=1}^{p} \beta_j \sigma_{t-j}^2, \qquad \omega > 0,\; \alpha_i, \beta_j \geq 0 \tag{GARCH}$

Covariance stationarity. A GARCH(1,1) process has finite unconditional variance if and only if $\alpha_1 + \beta_1 < 1$ . The unconditional variance is then $\sigma^2 = \omega / (1 - \alpha_1 - \beta_1)$ .

IGARCH. When $\alpha_1 + \beta_1 = 1$ the process is integrated in variance (IGARCH): shocks to variance are permanent and the unconditional variance does not exist. IGARCH still admits well-defined multi-step conditional variance forecasts converging to $+\infty$ at rate $h$ .

Multi-step forecasting (Engle 2001). For a covariance-stationary GARCH(1,1) with $\alpha_1 + \beta_1 < 1$ , the $k$ -step-ahead conditional variance forecast has the closed-form geometric mean-reversion recursion:

$\sigma^2_{t+k|t} = V_L + (\alpha_1+\beta_1)^k\bigl(\sigma^2_t - V_L\bigr), \qquad V_L = \frac{\omega}{1-\alpha_1-\beta_1}$

where $V_L$ is the long-run (unconditional) variance. For $k=1$ this reduces to the one-step GARCH recursion; as $k\to\infty$ forecasts converge to $V_L$ at geometric rate $\alpha_1+\beta_1$ . When persistence is near unity (e.g., $\alpha_1+\beta_1 = 0.99$ ), a volatility spike takes hundreds of periods to half-revert, producing the empirically observed long memory in squared returns without requiring fractional integration.

Empirical (Engle 2001 — Dow Jones Industrial Average (DJIA) + bond portfolio, 1990–2000): $\hat\omega = 0.00054$ , $\hat\alpha_1 = 0.0461$ , $\hat\beta_1 = 0.9461$ ; $\hat\alpha_1+\hat\beta_1 = 0.9922$ ; $\hat{V}_L = 0.0054/(1-0.9922) = 0.0114$ (daily SD $\approx 1.07\%$ ). The 1% daily Value at Risk (VaR) on a $1M portfolio ranged from ≈$8,400 (low-volatility periods) to ≈$52,900 (high-volatility periods) — a six-fold range driven entirely by time-varying $\sigma_t$ .

Component GARCH (Engle-Lee 1999)

The basic GARCH(1,1) model treats the long-run variance $V_L$ as a fixed constant. Engle and Lee (1999) allow it to be time-varying by writing the GARCH(2,2) model in a component form:

$\sigma_t^2 = q_t + \alpha(\varepsilon_{t-1}^2 - q_{t-1}) + \beta(\sigma_{t-1}^2 - q_{t-1})$ $q_t = \bar\omega + \rho\, q_{t-1} + \phi(\varepsilon_{t-1}^2 - \sigma_{t-1}^2)$

where $q_t$ is the long-run (trend) variance component and $(\sigma_t^2 - q_t)$ is the short-run (transitory) component. The long-run component $q_t$ follows its own AR(1) with persistence $\rho < 1$ , while $(\alpha + \beta)$ governs the speed of mean-reversion to $q_t$ . When $\rho \approx 1$ , the long-run component moves slowly and the model approximates long memory in volatility (Andersen-Bollerslev 1997), which is routinely found in asset return variances. Gallant, Hsu, and Tauchen (1999) and Alizadeh, Brandt, and Diebold (2002) provide empirical evidence for this component structure.

Multivariate BEKK

For an $n$ -dimensional innovation vector $\varepsilon_t \in \mathbb{R}^n$ with conditional covariance $H_t = \mathrm{Var}(\varepsilon_t | \mathcal{F}_{t-1})$ , the BEKK(1,1) model (Engle and Kroner 1995) is:

$H_t = C'C + A_1' \varepsilon_{t-1} \varepsilon_{t-1}' A_1 + G_1' H_{t-1} G_1 \tag{BEKK}$

where:

$C \in \mathbb{R}^{n \times n}$ is lower triangular with positive diagonal — the "long-run" component; $C'C$ is positive definite.
$A_1, G_1 \in \mathbb{R}^{n \times n}$ — ARCH and GARCH coefficient matrices. The outer-product form $A_1'(\cdot)A_1$ guarantees $H_t \succ 0$ for all $t$ whenever $C'C \succ 0$ .

Parameter count. Full BEKK(1,1) has $n(n+1)/2 + 2n^2$ parameters. For $n=2$ : 3 + 8 = 11 parameters.

Diagonal BEKK. Restrict $A_1 = \mathrm{diag}(a_1, a_2)$ and $G_1 = \mathrm{diag}(g_1, g_2)$ . This eliminates cross-equation volatility spillovers. For $n=2$ : 3 + 4 = 7 parameters. The diagonal restriction is testable by a likelihood-ratio statistic with $\chi^2_{2n(n-1)}$ degrees of freedom.

General BEKK( $P,Q$ ):

$H_t = C'C + \sum_{i=1}^{Q} A_i' \varepsilon_{t-i} \varepsilon_{t-i}' A_i + \sum_{j=1}^{P} G_j' H_{t-j} G_j$

Covariance Stationarity Condition

The BEKK(1,1) process is covariance stationary (has finite unconditional covariance matrix) if and only if all eigenvalues of

$\mathcal{A} \equiv A_1 \otimes A_1 + G_1 \otimes G_1$

lie strictly inside the unit circle, i.e., $\rho(\mathcal{A}) < 1$ . When $\rho(\mathcal{A}) \geq 1$ the unconditional covariance matrix $\mathrm{Vech}(H) = (I - \mathcal{A})^{-1}\mathrm{Vech}(C'C)$ does not exist — the multivariate analog of IGARCH.

Empirical finding (BDV 1998). Estimating BEKK-GARCH for bivariate long/short interest rate innovations in Belgium, Germany, France, UK, and USA, Bauwens, Deprins, and Vandeuren (1998) find $\rho(\mathcal{A}) > 1$ in all five countries, implying multivariate IGARCH behavior: unconditional interest rate volatility does not exist, though conditional forecasts are well-defined for short horizons of one to three months.

IGARCH and Exponentially Weighted Moving Average

When $\alpha_1 + \beta_1 = 1$ in GARCH(1,1), shocks to variance are permanent (integrated GARCH, Engle and Bollerslev 1986). The special case with $\omega = 0$ gives the EWMA formula:

$h_t = (1-\lambda)\varepsilon_{t-1}^2 + \lambda h_{t-1}, \qquad 0 < \lambda < 1 \tag{EWMA}$

which is an exponentially weighted moving average of past squared innovations with decay factor $\lambda$ . This is the formula underlying the RiskMetrics VaR methodology. Under IGARCH/EWMA the unconditional variance does not exist, but short-horizon conditional variance forecasts converge at rate $h \cdot \alpha_1$ per step.

Strict Stationarity vs. Covariance Stationarity (Nelson 1990b)

The covariance-stationarity condition $\alpha_1+\beta_1 < 1$ ensures finite unconditional variance but is not necessary for the process to be well-defined and stationary in distribution. Nelson (1990b) established the exact condition for strict stationarity of GARCH(1,1): the process is strictly stationary (and ergodic) if and only if

$E\bigl[\ln(\alpha_1 z^2 + \beta_1)\bigr] < 0 \tag{strict stationarity}$

where $z_t \overset{iid}{\sim}(0,1)$ . Because $E[\ln(\alpha_1 z^2+\beta_1)] \leq \ln E[\alpha_1 z^2+\beta_1] = \ln(\alpha_1+\beta_1)$ by Jensen's inequality, the strict-stationarity condition is strictly weaker than $\alpha_1+\beta_1 < 1$ . In particular, IGARCH ( $\alpha_1+\beta_1 = 1$ ) can satisfy the strict-stationarity criterion.

Persistence ambiguity. IGARCH is strictly stationary yet has $E_s(\sigma^q_t) \to \infty$ as $t\to\infty$ for $q \geq 1$ . This means variance shocks never die out in expectation even though the process has a well-defined stationary distribution. Strict stationarity and moment convergence are distinct notions: IGARCH is "stationary" in one sense and "non-stationary" in another. Any empirical finding of IGARCH should therefore distinguish between these senses — a model that is strictly stationary but has diverging conditional variance forecasts is economically very different from a covariance-stationary model.

Practical implication. Near-IGARCH estimates ( $\hat\alpha_1+\hat\beta_1 \approx 0.99$ ) are common. Nelson's theorem explains why: the strict-stationarity region extends well above $\alpha_1+\beta_1 = 1$ , so such processes can be legitimate probability models even when they lack finite unconditional variance. See also the spurious-IGARCH caveat in Open Questions.

Bayesian Stationarity Analysis of GARCH (Kleibergen-Van Dijk 1993)

Kleibergen and Van Dijk (1993) sharpen the weak vs. quasi-strict distinction using Bayesian inference on US T-bill data.

Weak vs. quasi-strict. The weak condition $\alpha + \beta < 1$ is strictly more restrictive than the quasi-strict condition $E[\ln(\beta + \alpha z^2)] < 0$ : the stationary region under the weak criterion is a proper subset of the quasi-strict stationary region.

Convergence theorem. The conditional variance product $G_k = G_{k-1}(\beta + \alpha z_k^2)$ converges to 0 a.s. iff $E[\ln(\beta + \alpha z^2)] < 0$ , and diverges to $+\infty$ otherwise. The mechanism is the law of large numbers on the log scale: $k^{-1}\ln G_k \to E[\ln(\beta + \alpha z^2)]$ a.s. The weak condition incorrectly uses $E[G_k] = (\alpha+\beta)^k$ and overestimates non-stationarity.

IGARCH revisited. At $\alpha + \beta = 1$ with normal errors, $\Pr[h_{t+1} < h_t] = 0.68$ : the conditional variance decreases more often than it increases, even though $E[h_{t+k}] \to \infty$ . IGARCH satisfies the quasi-strict condition despite violating weak stationarity.

Empirical application (US 3-month T-bill, Jan 1957–Apr 1989, $T=388$ ). AR(1)-GARCH(1,1)-t model estimated by importance sampling (SISAM, Hop-Van Dijk 1992). Bayesian posterior mean: $\hat\alpha = 0.18$ , $\hat\beta = 0.72$ , $\hat\lambda = 3.45$ .

Stationarity criterion	$\Pr(\text{non-stat})$	Bayes factor $K_{01}$
Weak ( $\alpha+\beta < 1$ )	$\approx 0.99$	4.95
Quasi-strict ( $E[\ln(\beta+\alpha z^2)] < 0$ )	$\approx 0.13$	0.434

The weak criterion overstates non-stationarity by nearly 8× relative to the quasi-strict criterion.

Fat-tail / unit-root interaction. Posterior odds $K_{01}$ for a unit root in the AR mean coefficient $H_0\colon \rho=1$ :

Model	$K_{01}$	Conclusion
Constant variance, normal	0.17	Reject unit root
GARCH(1,1), normal	0.04	Strongly reject
Constant variance, t	2.30	Support unit root
GARCH(1,1), t	2.50	Support unit root

The AR and tail-thickness parameters are negatively correlated in the posterior: fat-tailed models attribute extreme observations to heavy tails rather than high AR persistence, freeing $\rho$ to concentrate near 1. The t-score function $-\partial\ln p/\partial\varepsilon^2 = (\lambda+1)/[2(\lambda+\varepsilon^2)]$ decreases for large $|\varepsilon|$ , downweighting outliers. Normal likelihoods treat the same outliers as high-volatility evidence, biasing the posterior toward stationarity in mean.

Continuous-Time Limits of ARCH Models

Different ARCH families, when sampled at increasingly fine time intervals, converge to different continuous-time diffusions. This establishes structural links between ARCH models and stochastic-volatility option-pricing models.

GARCH(1,1) limit (Nelson 1990a; Drost-Nijman 1993). As the sampling interval $h \to 0$ with GARCH parameters rescaled appropriately, the GARCH(1,1) variance process converges to the mean-reverting stochastic differential equation (SDE):

$d\sigma^2 = \theta(\omega - \sigma^2)\,dt + \sigma^2\,dW_\sigma$

where $\theta$ is the speed of mean-reversion and $dW_\sigma$ is a Brownian motion. This is a special case of the Heston (1993) stochastic-volatility model (with zero correlation between $dW_\sigma$ and the return Brownian motion). The unconditional distribution of the GARCH(1,1) process converges to an inverted-gamma mixture of normals, producing Student- $t$ tails — consistent with the fat tails documented in asset return data.

EGARCH limit. The EGARCH (Exponential GARCH) process converges to a diffusion where $\ln\sigma^2$ follows an Ornstein-Uhlenbeck process, giving a log-normal unconditional distribution for $\sigma^2$ . The resulting returns follow a normal–lognormal mixture. This provides a theoretical connection between EGARCH and the log-normal stochastic volatility model used in continuous-time finance.

ARCH-family selectivity. Different ARCH specifications (power-GARCH, TARCH, etc.) converge to different diffusions. The choice of ARCH model is therefore equivalent to choosing a specific continuous-time volatility specification, which matters for option pricing (Engle-Mustafa 1992) and for the interpretation of the unconditional return distribution.

Option pricing via diffusion limits (Bollerslev-Engle-Nelson 1994). Because the continuous-time limit of a GARCH process has a known form, the Black-Scholes PDE evaluated at the GARCH-implied diffusion provides a structural link between ARCH-based conditional variance estimates and option prices — a connection exploited in ARCH-based option pricing (Engle-Mustafa 1992; Duan 1995).

Temporal Aggregation (Drost-Nijman 1993)

A natural question is whether a GARCH model at daily frequency implies GARCH at weekly or monthly frequency. The answer depends on whether "GARCH" is interpreted in the strict or weak sense.

Weak GARCH. The weak GARCH class (Drost-Nijman 1993) defines conditional heteroscedasticity only through the second-order moment structure — it does not require the conditional distribution to be Gaussian or even i.i.d. innovations. Weak GARCH is closed under temporal aggregation: if $\{r_t\}$ is weak GARCH( $p,q$ ) at frequency $h$ , then the aggregated series $R_T = \sum_{t=(T-1)h+1}^{Th} r_t$ is weak GARCH at frequency $H = nh$ , with parameters that are explicit functions of the original parameters and $n$ .

The aggregation formula implies that the persistence $\alpha_1+\beta_1$ is preserved exactly under aggregation: if $\alpha_1+\beta_1 = 0.99$ daily, the weekly series has the same persistence. What changes is the allocation between the ARCH and GARCH components — aggregation concentrates persistence into the GARCH term $\beta$ and reduces the ARCH term $\alpha$ .

Heteroscedasticity fading. As $n \to \infty$ (very low frequency), the ARCH coefficient of the aggregated series $\to 0$ and the GARCH coefficient $\to \alpha_1+\beta_1$ of the original, so the conditional variance becomes nearly constant at the long-run value $V_L$ . If $\alpha_1+\beta_1 < 1$ , heteroscedasticity disappears at sufficiently low frequencies — consistent with the empirical observation (Akgiray 1989; BEN 1994) that daily returns show strong GARCH while monthly returns are approximately i.i.d. normal.

Strict GARCH is not closed. If the daily process is strict GARCH (conditional distribution Gaussian), the weekly aggregate is generally not strict GARCH. Aggregation introduces leptokurtosis at the weekly level that cannot be fully captured by a GARCH model with Gaussian innovations. This is why Gaussian GARCH tends to underestimate kurtosis at intermediate frequencies.

Optimal ARCH Filter (Nelson-Foster 1994)

Given an underlying continuous-time diffusion $d\sigma^2_t = \mu_t\,dt + v_t\,dW_t$ , what discrete-time ARCH model minimises the asymptotic mean-squared error of the resulting volatility estimate?

Optimality criterion. Nelson and Foster (1994) frame this as a filter-design problem: among all ARCH processes of the form $h_t = f(\varepsilon_{t-1}, h_{t-1}, \theta)$ , find the one for which the asymptotic estimation error variance $E[(h_t - \sigma^2_t)^2]$ is minimised. The answer depends on the assumed shape of the conditional distribution.

Result: absolute-value GARCH. When the true conditional distribution of $\varepsilon_t/\sigma_t$ has heavy tails (e.g., Student- $t$ or any distribution with excess kurtosis), the optimal filter uses absolute-value residuals $|\varepsilon_{t-1}|$ rather than squared residuals $\varepsilon^2_{t-1}$ . This provides a formal theoretical justification for TARCH (Zakoïan 1991) and NARCH/Power-ARCH (Higgins-Bera 1992) over standard GARCH when innovations are fat-tailed.

Intuition. Large squared residuals from heavy-tailed distributions carry excessive noise relative to their information content. Absolute-value residuals weight large observations less aggressively, reducing the estimator's variance. In the extreme (Cauchy innovations), squared residuals are meaningless, and only absolute values carry useful information about scale.

Empirical relevance. BEN (1994) and Hentschel (1995) find estimated power parameters $\hat\nu \approx 1.5$ on US equity data — intermediate between absolute-value ( $\nu=1$ ) and squared ( $\nu=2$ ) — consistent with Nelson-Foster optimality for moderately fat-tailed distributions.

Co-Persistence in Variance (Bollerslev and Engle 1993)

Just as I(1) series can be co-integrated (a linear combination is stationary in the mean), IGARCH series can be co-persistent in variance: a portfolio $\gamma'y_t$ has no variance persistence even though each individual series does. The co-persistent vector $\gamma$ plays the exact role of the co-integrating vector.

Formal definition. Let $H_t^*(s) = E_t(\text{vech}(H_t)) - E_s(\text{vech}(H_t))$ measure the influence of time- $s$ shocks on future covariance forecasts. The process is persistent in variance if $\limsup_{t\to\infty} H_t^*(s) \neq 0$ a.s. It is co-persistent if a nonzero $\gamma$ exists such that $\limsup_{t\to\infty} |E_t(\gamma'H_t\gamma) - E_s(\gamma'H_t\gamma)| = 0$ a.s.

Theorem 2. Vector GARCH( $p,q$ ) is co-persistent iff $\text{vec}_2(\gamma)'v_i = 0$ for all right eigenvectors $v_i$ corresponding to $|\lambda_i| \geq 1$ . The co-persistent portfolio eliminates the explosive subspace.

Factor GARCH connection (Lemma 2). In a $K$ -factor GARCH model with $r < K$ integrated factors and $K-r$ stationary factors, any portfolio orthogonal to the integrated factors is co-persistent. Co-persistence exists whenever the number of persistent variance factors is less than $N$ .

Empirical finding. Daily DM and BP vs. USD (1980–1985, 1,245 obs): both near-IGARCH ( $\hat\alpha+\hat\beta \approx 0.97$ ); estimated co-persistent vector $\hat\gamma \approx (1,-1)$ — the DM/BP bilateral rate. Univariate GARCH for DM/BP: t-stat for $\hat\alpha+\hat\beta=1$ equals 6.033 — strongly rejects IGARCH. Interpretation: persistence in USD exchange rate volatility is dollar-specific news; the bilateral European rate shows no variance persistence.

Empirical Evidence — Stock Returns (Akgiray 1989)

One of the first rigorous GARCH applications. Data: CRSP value-weighted daily index returns, Jan 1963 – Dec 1986 (T = 6,030), split into four 6-year sub-periods.

Key diagnostic finding. After AR(1) pre-filtering ( $R_t = \phi_0 + \phi_1 R_{t-1} + e_t$ , $\phi_1 = 0.18$ – $0.31$ ), standard independence tests on $\{e_t\}$ fail to reject white noise. But $\{|e_t|\}$ and $\{e_t^2\}$ remain significantly autocorrelated at all lags up to 60 — conclusive evidence of nonlinear (non-strict-white-noise) dependence. This is the diagnostic argument for GARCH over linear models.

Estimation. GARCH(1,1) strictly dominates ARCH( $p^*$ ) ( $p^* = 2$ – $5$ ) by log-likelihood in every period. Representative GARCH(1,1) estimates: $\alpha_1 \approx 0.05$ – $0.09$ , $\beta \approx 0.76$ – $0.92$ , $\alpha_1+\beta \approx 0.96$ – $0.99$ . Dickey-Fuller rejects $\alpha_1+\beta=1$ in 3 of 4 periods (covariance-stationary, but near-IGARCH). $\beta \gg \alpha_1$ in all periods — past conditional variance dominates new shocks, producing persistent volatility. Standardized residuals pass normality tests; the GARCH model fully absorbs the excess kurtosis of raw returns.

The autocorrelation function (ACF) of $\{e_t^2\}$ closely follows the GARCH(1,1) recursion $c_n = (\alpha_1+\beta) \cdot c_{n-1}$ — 48 of 60 autocorrelations lie within $\pm 10\%$ of the implied values.

Volatility forecasting comparison. 24 rolling out-of-sample monthly variance forecasts per period, compared across four methods:

Method	Description
Historical mean	Unconditional sample variance over estimation window
EWMA	Exponentially weighted moving average, smoothing constant w estimated from data (w ≈ 0.11–0.24)
ARCH	h-step-ahead forecast from fitted ARCH(p) via iterated expectations
GARCH	Same formula, $A = \alpha_1+\beta$ , from GARCH(1,1)

GARCH wins on mean error (ME), root-mean-square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE) in all four periods. Gains are largest in high-volatility periods (1969-74, 1975-80). Even so, minimum MAPE > 30% across all methods — variance forecasting is inherently difficult.

Frequency dependence. GARCH effects are a daily phenomenon: monthly returns are approximately i.i.d. normal (no significant autocorrelation in {e²} at monthly frequency), likely reflecting a central limit theorem (CLT) for sums of dependent daily returns.

ARCH-in-Mean

The ARCH-M model (Engle, Lilien, and Robins 1987) augments the conditional mean with a function of conditional variance:

$r_t = \mu + \delta\, g(h_t) + \varepsilon_t, \qquad h_t = \omega + \alpha_1\varepsilon_{t-1}^2 + \beta_1 h_{t-1} \tag{ARCH-M}$

Common choices: $g(h_t) = h_t^{1/2}$ (standard deviation), $g(h_t) = h_t$ (variance), or $g(h_t) = \ln h_t$ . The parameter $\delta > 0$ captures the risk premium: investors demand higher expected returns when conditional risk is higher. Engle et al. (1987) found a significant $\delta$ for excess returns of 6-month US Treasury bills. This model is the first formal link between ARCH volatility and the expected return level.

Conditional Variance ARMA Structure (Fiorentini-Sentana 1998)

Fiorentini and Sentana (1998) show that when $[I-A(L)]\mathbf{x}_t = [I-B(L)]\boldsymbol{\epsilon}_t$ , the conditional mean $\boldsymbol{\mu}_t = E_{t-1}(\mathbf{x}_t)$ itself follows a vector linear process with the same AR polynomial:

$[I-A(L)]\boldsymbol{\mu}_t = [A(L)-B(L)]\boldsymbol{\epsilon}_t$

Applied to multivariate GARCH $(p,q)$ , the conditional variance $\text{vech}(\Sigma_t)$ is a stationary ARMA process for past squared innovations and past variances; the Fiorentini-Sentana result pins down the ACF structure. For GARCH $(p,q)$ the conditional variance $\text{vech}(\Sigma_t)$ displays the ACF of a vector autoregressive moving-average (VARMA $(m, q-1)$ ) process with the same AR polynomial and $m = \max(p,q)$ . The GARCH(1,1) special case is particularly clean:

$\text{vech}(\Sigma_{t+1}) = \alpha_0 + (\alpha_1 + \beta_1)\,\text{vech}(\Sigma_t) + \alpha_1\,\mathbf{v}_t$

where $\mathbf{v}_t = \text{vech}(\epsilon_t\epsilon_t' - \Sigma_t)$ is the innovation to squared returns. This is a vector autoregression (VAR(1)) for the conditional covariance matrix — the one-period persistence is $\alpha_1 + \beta_1$ exactly, consistent with the GARCH recursion.

GARCH-M Nonneg-ACF Constraint (Fiorentini-Sentana 1998)

The GARCH-M model augments the conditional mean with conditional variance: $x_t = \delta\sigma_t^2 + \epsilon_t$ . Because $\mu_t = \delta\sigma_t^2$ is proportional to conditional variance, the autocorrelation structure of $x_t$ is entirely determined by the autocorrelation structure of $\sigma_t^2$ , which is always nonneg:

$\text{Cov}(x_t, x_{t-k}) = \delta^2\,\text{Cov}(\sigma_t^2, \sigma_{t-k}^2) \geq 0, \qquad \forall\, k \geq 1$

This is a testable structural restriction: GARCH-M processes can only generate nonneg autocorrelations in the observable series $x_t$ . Financial returns routinely exhibit near-zero or weakly negative autocorrelations, which explains why GARCH-M models tend to fit financial time series poorly despite their theoretical appeal as risk-return models. The constraint is exact (not approximate) and applies to all parametrisations of the variance function $g(h_t)$ as long as $x_t = \delta g(h_t) + \epsilon_t$ with $g > 0$ .

Multivariate vech-GARCH

Bollerslev, Engle, and Wooldridge (1988) proposed the vech model for an $n$ -vector of returns $r_t$ :

$\mathrm{vech}(H_t) = c + \sum_{j=1}^{q} A_j\, \mathrm{vech}(\varepsilon_{t-j}\varepsilon_{t-j}') + \sum_{j=1}^{p} B_j\, \mathrm{vech}(H_{t-j}) \tag{vech}$

where $\mathrm{vech}(\cdot)$ stacks the lower triangle of a symmetric matrix ( $n(n+1)/2$ elements). Diagonal restrictions on $A_1$ and $B_1$ reduce the parameter count but eliminate cross-variance dynamics. With $n=3$ (bills, bonds, stocks), BEW show that Capital Asset Pricing Model (CAPM) $\beta$ -coefficients are time-varying — the conditional covariance with the market portfolio fluctuates with $H_t$ .

A practical limitation: the vech model does not automatically guarantee $H_t \succ 0$ ; extra parameter constraints are required. This motivated the BEKK form (see above).

Factor-ARCH

For large $n$ (many assets), the BEKK and vech parameter counts become prohibitive. The factor-ARCH model (Engle 1987; Engle, Ng, and Rothschild 1990) imposes a factor structure:

$r_t = B\xi_t + \varepsilon_t, \qquad \mathrm{Var}(\varepsilon_t|\mathcal{F}_{t-1}) = \Psi, \qquad \mathrm{Var}(\xi_t|\mathcal{F}_{t-1}) = \Lambda_t$

where $B$ ( $n\times k$ ) is a loadings matrix, $k \ll n$ is the number of common volatility factors, $\Psi$ is constant ( $n\times n$ ), and $\Lambda_t = \mathrm{diag}(\lambda_{1t},\ldots,\lambda_{kt})$ with each $\lambda_{jt}$ following an ARCH process. The conditional covariance matrix is:

$H_t = \Psi + B\Lambda_t B' = \Psi + \sum_{j=1}^{k} \lambda_{jt}\,\beta_j\beta_j'$

Common volatility factors drive all cross-asset covariances simultaneously; only $k$ ARCH equations (plus $\Psi$ and $B$ ) need to be estimated. Applied to Diebold-Nerlove (1989) daily exchange rates and to Treasury bills of different maturities (Engle et al. 1990).

Constant Conditional Correlation (CCC)

Bollerslev (1990) proposed a computationally tractable multivariate GARCH by restricting conditional correlations to be constant through time:

$H_t = D_t \Gamma D_t$

where $D_t = \mathrm{diag}(\sigma_{1t},\ldots,\sigma_{nt})$ contains individual GARCH standard deviations and $\Gamma$ is a time-invariant positive definite correlation matrix. The conditional covariance between series $i$ and $j$ is $h_{ijt} = \rho_{ij}\sigma_{it}\sigma_{jt}$ — time-varying level driven by each series' own GARCH volatility, but constant relative coherence.

Estimation. The ML estimate of $\Gamma$ is the sample correlation of standardized residuals: $\hat\Gamma = T^{-1}\sum_t \tilde\varepsilon_t\tilde\varepsilon_t'$ where $\tilde\varepsilon_t = D_t^{-1}e_t$ , positive definite by construction. This concentrates $\Gamma$ out of the likelihood, reducing $T$ matrix inversions to 1. The concentrated log-likelihood reduces to:

$\ell_c = \mathrm{const} - \frac{1}{2}\sum_{t=1}^T \Bigl[2\log|D_t| + \log|\hat\Gamma| + \tilde\varepsilon_t'\hat\Gamma^{-1}\tilde\varepsilon_t\Bigr]$

SUR connection. With constant correlations the model is an extension of Seemingly Unrelated Regressions (SUR) allowing GARCH heteroskedasticity — a framing that clarifies why $\hat\Gamma$ has the SUR sample analogue form.

Empirical finding (Bollerslev 1990). Applied to weekly DM, FF, IL, SF, BP returns vs. USD: all conditional correlations significantly higher in the European Monetary System (EMS) period (1979–1985) than in the pre-EMS float (1973–1979). E.g. DM–FF: 0.607 → 0.932. LR test for parameter constancy: 388 ~ $\chi^2(30)$ — highly significant. Even non-EMS currencies (BP, SF) showed higher correlations, suggesting a common USD factor beyond EMS policy coordination.

Limitation. Constant correlations is a strong assumption relaxed by DCC (Engle 2002) below.

Dynamic Conditional Correlation (DCC)

Engle (2002a) introduced the Dynamic Conditional Correlation (DCC)-GARCH model as a computationally tractable multivariate GARCH with time-varying correlations. Write $r_t = H_t^{1/2} z_t$ and decompose $H_t = D_t R_t D_t$ where $D_t = \mathrm{diag}(h_{11t}^{1/2},\ldots,h_{nnt}^{1/2})$ contains individual GARCH standard deviations and $R_t$ is the conditional correlation matrix. $R_t$ evolves as:

$Q_t = \bar Q(1-a-b) + a\, u_{t-1}u_{t-1}' + b\, Q_{t-1}, \qquad R_t = \mathrm{diag}(Q_t)^{-1/2} Q_t\, \mathrm{diag}(Q_t)^{-1/2}$

where $u_t = D_t^{-1} r_t$ are standardized residuals and $\bar Q = E[u_t u_t']$ is the unconditional correlation. DCC extends the CCC model (Bollerslev 1990) where $R_t = \Gamma$ is fixed. DCC is estimated in two stages: (1) fit $n$ univariate GARCH models for the variances; (2) estimate $(a, b, \bar Q)$ from the standardized residuals.

Multivariate GARCH: Taxonomy and Properties (Bauwens-Laurent-Rombouts 2006)

Bauwens, Laurent, and Rombouts (2006) survey the MGARCH literature by organizing all models into three families and documenting their structural properties. See Bauwens-Laurent-Rombouts (2006).

Family 1 — Direct Generalizations

VEC(1,1) (Bollerslev-Engle-Wooldridge 1988). The full vector model is $h_t = c + A\eta_{t-1} + Gh_{t-1}$ where $h_t = \mathrm{vech}(H_t)$ , $\eta_t = \mathrm{vech}(\varepsilon_t\varepsilon_t')$ , and $A$ , $G$ are $[N(N+1)/2] \times [N(N+1)/2]$ matrices. Parameter count for the full model: $[N(N+1)(N(N+1)+1)]/2$ — for $N=3$ this is 78, already unidentifiable in practice. Diagonal VEC restricts $A$ and $G$ to diagonal, giving $[N(N+5)]/2$ parameters.

RiskMetrics. Scalar diagonal VEC with pre-set $\lambda = 0.94$ (daily) or $\lambda = 0.97$ (monthly): $H_t = (1-\lambda)\varepsilon_{t-1}\varepsilon_{t-1}' + \lambda H_{t-1}$ . Not estimated — calibrated by J.P. Morgan.

BEKK(1,1,K) (Engle-Kroner 1995). $H_t = C^{*\prime}C^* + \sum_{k=1}^K A_k^{*\prime}\varepsilon_{t-1}\varepsilon_{t-1}'A_k^* + \sum_{k=1}^K G_k^{*\prime}H_{t-1}G_k^*$ . PD by construction. For scalar $K=1$ : $[N(5N+1)]/2$ parameters. Covariance-stationary iff all eigenvalues of $A\otimes A + G\otimes G$ lie strictly inside the unit disk (BLR notation).

Factor GARCH. BEKK with rank-1 matrices $A_k = a_k b_k'$ , $G_k = g_k b_k'$ ; each factor $h_{kt} = \omega_k + \alpha_k (b_k'\varepsilon_{t-1})^2 + \beta_k h_{k,t-1}$ is a univariate GARCH(1,1). Common factors drive all cross-asset covariances. For $K=1$ factor: $[N(N+5)]/2$ parameters (same as diagonal VEC). All factors share persistence implied by the respective $\alpha_k + \beta_k$ .

FF-GARCH (Vrontos-Dellaportas-Politis 2003). $H_t = W\Sigma_t W'$ where $W$ is lower-triangular with ones on the diagonal and $\Sigma_t = \mathrm{diag}(\sigma_{1t}^2,\ldots,\sigma_{Nt}^2)$ contains individual GARCH variances. Always PD. Generalizes the Cholesky-based Pourahmadi (1999) reparametrization.

Family 2 — Linear Combinations of Univariate GARCH

O-GARCH (principal component GARCH). Extract $m \leq N$ eigenvectors from the sample correlation matrix of returns; fit univariate GARCH to each principal component. The conditional covariance is $H_t = PF_tP'$ where $P$ is the $N\times m$ eigenvector matrix and $F_t$ is the $m\times m$ diagonal matrix of factor variances. Simple and scalable, but orthogonality of components is a restrictive assumption.

GO-GARCH (van der Weide 2002). Generalizes O-GARCH by allowing the $N\times N$ mixing matrix to be $\Lambda = PL^{1/2}U$ where $P$ is the $N\times N$ eigenvector matrix, $L$ is the diagonal eigenvalue matrix, and $U$ is an orthogonal matrix parametrized by $N(N-1)/2$ Givens rotation angles. Two-step estimation: first estimate $PL^{1/2}$ by PCA, then estimate $U$ by ML. Nests O-GARCH as $U = I$ .

Family 3 — Nonlinear Combinations

CCC (Bollerslev 1990) — recap. $H_t = D_t R D_t$ with constant $R$ . The constant-correlation assumption is testable via the LM test of Tse (2000): $\mathcal{LM} \sim \chi^2(N(N-1)/2)$ under $H_0$ of constant $R$ . CCC is not invariant to linear transformations of $\varepsilon_t$ — multiplying by a non-diagonal matrix changes the model structure.

DCC-T (Tse-Tsui 2002). $R_t = (1-\theta_1-\theta_2)R + \theta_1\Psi_{t-1} + \theta_2 R_{t-1}$ where $\Psi_{t-1}$ is the $M$ -period rolling correlation matrix of $\{u_{t-j}\}_{j=1}^M$ and $u_t = D_t^{-1}\varepsilon_t$ . Requires $M \geq N$ for $\Psi_{t-1}$ to be invertible. The baseline $R$ and two scalars $\theta_1, \theta_2$ are free.

DCC-E (Engle 2002) — extended version. Scalar DCC uses two parameters $(\alpha, \beta)$ ; extended DCC uses $N\times N$ matrices $A$ , $B$ so $Q_t = \bar{Q}\odot(11' - A - B) + A\odot u_{t-1}u_{t-1}' + B\odot Q_{t-1}$ , preserving PD via Hadamard products. Scalar DCC has $[(N+1)(N+4)]/2$ parameters including the $N(N+1)/2$ entries of $\bar{Q}$ .

GDC (Kroner-Ng 1998). $H_t = D_t R_t D_t + \Phi\odot\Theta_t$ where $D_t$ contains GARCH standard deviations, $R_t$ is a dynamic correlation matrix, and $\Phi\odot\Theta_t$ adds element-wise cross-product terms. The ADC (Asymmetric Dynamic Conditional Covariance) extends GDC with leverage: the $ij$ -element is augmented by $b_i'v_{t-1}v_{t-1}'b_j$ where $v_t = \max(0,-\varepsilon_t)$ (component-wise). GDC nests DCC, CCC, DVEC, BEKK, and Factor GARCH; $[N(7N-1)+4]/2$ free parameters.

Copula-GARCH (Patton 2000, 2002; Jondeau-Rockinger 2001). Sklar's (1959) theorem: any multivariate distribution can be written as $F(\varepsilon_1,\ldots,\varepsilon_N) = C(F_1(\varepsilon_1),\ldots,F_N(\varepsilon_N))$ where $C:[0,1]^N\to[0,1]$ is the copula and $F_i$ are marginals. Copula-GARCH fits $N$ univariate GARCH models for the margins and a (possibly time-varying) copula for the dependence structure. Time-varying copula parameters evolve as ARMA-type recursions in $u_t = F_t^{-1}(\varepsilon_t)$ . Two-step MLE naturally decomposes marginal and dependence likelihoods.

Structural Properties

Invariance. A model is invariant if multiplying $\varepsilon_t$ by a non-singular matrix $F$ yields a model of the same type for $F\varepsilon_t$ . General VEC and BEKK are invariant; diagonal VEC and BEKK are NOT (cross-covariance dynamics change). CCC is NOT invariant because $FD_t$ is generally not diagonal.

Marginalization. A bivariate VEC(1,1) model implies each univariate margin is at most a weak GARCH(3,3) — not necessarily a strong GARCH. Diagonal VEC(1,1) implies each margin is a strong GARCH(1,1). CCC, DCC, and GDC marginalization properties are unknown.

Temporal aggregation. Weak multivariate GARCH is closed under temporal aggregation (Hafner 2003): if the daily process is weak MGARCH, the weekly aggregate is also weak MGARCH with parameters that are explicit functions of the daily parameters. This generalizes the Drost-Nijman (1993) univariate result.

QML Theory and Two-Step Estimation

QML consistency (Jeantheau 1998). The Gaussian quasi-log-likelihood $L_T(\theta) = -\frac{1}{2}\sum_t[\log|H_t(\theta)| + \varepsilon_t'H_t(\theta)^{-1}\varepsilon_t]$ produces a consistent estimator whenever the conditional mean and covariance are correctly specified, regardless of the true innovation distribution. The analytical score (eq. 56 of BLR) avoids numerical differentiation.

QML asymptotic normality. Proved only for BEKK by Comte and Lieberman (2003). The asymptotic distribution theory for CCC, DCC, and GDC remains an open problem.

Variance targeting (Engle-Mezrich 1996). Replace the long-run intercept of $H_t$ (e.g., $C^{*\prime}C^*$ in BEKK or $\bar{Q}$ in DCC) by the sample unconditional covariance $\hat{H} = T^{-1}\sum_t\varepsilon_t\varepsilon_t'$ . Reduces the free-parameter count substantially without affecting consistency; makes the remaining parameters identifiable from the correlation dynamics alone.

DCC two-step. The full log-likelihood decomposes into $L_T = QL_1 + QL_2$ where: $QL_1 = \sum_{i=1}^N \ell_i(\theta_i) = -\frac{1}{2}\sum_t\sum_i\bigl[\log h_{ii,t} + \varepsilon_{it}^2/h_{ii,t}\bigr]$ $QL_2 = -\frac{1}{2}\sum_t\bigl[\log|R_t(\theta_c)| + u_t'R_t(\theta_c)^{-1}u_t - u_t'u_t\bigr]$ Step 1 maximizes $QL_1$ for each margin separately; step 2 maximizes $QL_2$ for correlation parameters $\theta_c$ given the step-1 estimates. The two-step estimator is consistent but loses efficiency relative to joint maximization; one Newton-Raphson step at the two-step estimates achieves asymptotic efficiency.

CCC-MGARCH with Structural Break: Korean Currency Crisis (Kim-Tsurumi 2000)

Kim and Tsurumi (2000) apply Bollerslev's (1990) CCC-MGARCH to four daily Korean financial returns — KOSPI spot (KS), KOSPI futures (KF), Won/Dollar spot (WS), and 3-month NDF (NF) — over 10 October 1996 – 1 April 1998. Individual GARCH(1,1) estimates yield $\hat\alpha_i + \hat\beta_i \approx 1$ (IGARCH) for all four series. Full-sample $\hat\rho_{KF,KS} = 0.8597$ ; pooled log-likelihood = −2100.77.

Structural break date $t^*$ is treated as an unknown parameter. The Bayesian log-posterior is approximated via Laplace: $\ln g(t^* \mid \hat\theta, y) \propto \sum_{j=1}^{2} \Bigl(L(\hat\theta_j) - \tfrac{1}{2}\ln\bigl|H(t^*_j)\bigr|\Bigr)$ where $L(\hat\theta_j)$ and $H(t^*_j)$ are split-sample log-likelihoods and Hessians at the respective MLE. Grid search over 5 July – 10 December 1997 identifies the MGARCH posterior break mode at 20 October 1997 — three weeks before the 8 November official devaluation (pre-break LL = −889.14; post-break LL = −1144.72). Individual-series posteriors find earlier breaks for derivative markets: KF = 11 August 1997; NF = 14 August 1997 — consistent with informed trading anticipating the crisis. The IGARCH finding is consistent with the Hillebrand (2004) / Lamoureux-Lastrapes (1990) result that structural breaks in GARCH parameters can induce spurious near-unit-root variance persistence. See Kim-Tsurumi (2000) and Structural Break Testing.

Diagnostics for Multivariate GARCH

Given standardized residuals $\hat{u}_t = \hat{H}_t^{-1/2}\hat{\varepsilon}_t$ and $\hat{z}_t = \mathrm{vech}(\hat{u}_t\hat{u}_t' - I_N)$ :

Moment conditions A/B/C: Test $E[\hat{u}_t]=0$ , $E[\hat{u}_t\hat{u}_t'-I]=0$ , $E[\hat{z}_t\otimes\hat{z}_t-I_{N^*(N^*+1)/2}]=0$ where $N^* = N(N+1)/2$ .
$HM(M)$ (Hosking 1980): Multivariate Ljung-Box for $\hat{u}_t$ : $HM(M) \sim \chi^2(K^2 M)$ under $H_0$ of no serial correlation.
$LL(M)$ (Ling-Li 1997): Distribution-free portmanteau for $\hat{u}_t^2$ ; $LL(M) \sim \chi^2(M)$ . Does not assume a specific innovation distribution.
$RB(M)_1$ and $RB(M)_2$ (Tse 2002): Residual-based tests on the off-diagonal and diagonal elements of $\hat{u}_t\hat{u}_t'$ ; each $\sim \chi^2(M)$ .
LM test for constant correlation (Tse 2000): $\mathcal{LM} \sim \chi^2(N(N-1)/2)$ ; tests $H_0:$ CCC against DCC-type alternatives.
DCC test (Engle-Sheppard 2001): Auxiliary regression of cross-products of standardized residuals on lagged products; tests whether $R_t$ varies over time.

Dimensionality Reduction: Flex-GARCH and Base-Asset Factor (Andersen et al. 2004)

For large N, the DCC model is the most flexible tractable option (N+2 parameters) but still requires N univariate GARCH estimations. An intermediate approach is Flex-GARCH (Ledoit, Santa-Clara, and Wolf 2003): estimate N(N+1)/2 bivariate GARCH(1,1) models with parameter constraints, and "paste" them together into the coefficient matrices A, B, C of a scalar GARCH. Specific transformations of bivariate parameters ensure psd covariance forecasts. Flex-GARCH becomes cumbersome above N ≈ 30; beyond that, the base-asset factor approach is preferred.

The base-asset mapping sets: $r_t = \Lambda f_t + \varepsilon_t$ where $f_t \in \mathbb{R}^{N_F}$ ( $N_F \ll N$ ) are returns on $N_F$ liquid base assets (e.g., equity indices, FX rates, benchmark rates), $\Lambda$ is an $N \times N_F$ loading matrix (from regression, pricing model, or expert matching), and $\varepsilon_t$ captures idiosyncratic risk. A realized $N_F \times N_F$ covariance matrix for the base assets is feasible because $N_F < 1/D$ (where $D^{-1}$ is the intraday sampling frequency, e.g., 48 for 30-minute data). Filtered Historical Simulation (see Value at Risk) then resamples factor and idiosyncratic innovations separately.

The R² Puzzle and Volatility Forecast Evaluation (Andersen-Bollerslev 1998)

Standard forecast evaluation regresses ex-post squared daily returns $r_t^2$ on the GARCH one-day-ahead conditional variance forecast $\hat\sigma_t^2$ . The resulting R² values of 0.02–0.10 led many researchers to conclude GARCH forecasts are useless. Andersen and Bollerslev (1998) showed this conclusion is wrong.

Analytical population R². Under correctly specified GARCH(1,1) with excess kurtosis $\kappa$ , the population R² of the Mincer-Zarnowitz regression ( $r_t^2 = a + b\hat\sigma_t^2 + u_t$ ) is:

$R^2(m=1) = \frac{\mathrm{Var}(\hat\sigma_t^2)}{\mathrm{Var}(r_t^2)} = \frac{\alpha^2/(1-\beta)^2}{\alpha^2/(1-\beta)^2 + \kappa\sigma^4/T + \cdots}$

For DM- $($ \alpha=0.068 $,$ \beta=0.898 $,$ \kappa=2.8 $) this gives$ R^2=0.064$, matching the empirical 0.047. The "disappointing" R² is exactly what a correctly specified model predicts — not evidence of misspecification.

Measurement error decomposition. The actual forecast MSE using $r_t^2$ decomposes as:

$\mathrm{MSE}(r_t^2) = \mathrm{MSE}^*(\sigma_t^2) + \underbrace{\mathrm{Var}(r_t^2 - \sigma_t^2)}_{\text{measurement noise}}$

For DM- $: actual MSE = 1.221 vs. ideal MSE = 0.084 — measurement noise accounts for 93% of the residual variance. Evaluating GARCH against$ r_t^2$ is almost entirely measuring noise, not forecast quality.

Realized volatility as the solution. By quadratic variation theory (Karatzas-Shreve 1988), $RV_t^{(m)} = \sum_{j=1}^m r_{t,j/m}^2 \to \int_0^1 \sigma_{t+s}^2\,ds$ as $m\to\infty$ . The noise terms shrink by a factor $m$ . At 5-minute frequency (m=288), R² rises from ≈0.05 to ≈0.48 (theory 0.483) for DM- $— a 10× improvement. The theoretical R² under continuous sampling converges to ≈0.495, determined solely by the persistent component$ (\alpha+\beta)$.

Continuous-time GARCH limit (Nelson 1990). GARCH(1,1) converges in distribution to the SDE $d\sigma^2 = \theta(\omega-\sigma^2)\,dt + (2\lambda\theta)^{1/2}\sigma^2\,dW_\sigma$ , establishing the structural link between discrete GARCH parameters and continuous-time integrated variance. Temporal aggregation uses the weak GARCH framework of Drost-Nijman (1993) and Drost-Werker (1996). See Andersen-Bollerslev (1998).

Realized Variance and Covariance (Andersen-Bollerslev-Diebold-Labys 2001–2003)

For assets with liquid intraday data, the realized variance on day $t$ is:

$RV_t = \sum_{j=1}^{D} r^2_{t,j/D}$

where $r_{t,j/D}$ are returns sampled at the $j$ -th interval on day $t$ (e.g., 30-minute intervals, $D = 48$ ). As $D \to \infty$ (continuous sampling), $RV_t$ converges to the integrated variance $\int_0^1 \sigma^2_{t+u}\,du$ , the true daily variance of the underlying continuous-time process. At 15–30 min intervals, measurement error from market microstructure noise is small enough to treat $RV_t$ as approximately observed.

Key empirical regularities:

$\log(RV_t)$ is approximately lognormally distributed, implying that modeling log-realized variance is more natural than modeling realized variance directly
Sample autocorrelations of $RV_t$ are significantly positive out to 100+ lags — long memory — motivating Autoregressive Fractionally Integrated Moving Average (ARFIMA) rather than short-memory autoregressive moving average (ARMA) models for realized variance
Returns standardized by $\sqrt{RV_t}$ are approximately standard normal (Zhou 1996; ABDL 2001), in sharp contrast to GARCH-standardized returns which still show fat left tails

Realized covariance generalizes to the multivariate case: $RC_t = \sum_{j=1}^{D} r_{t,j/D}\, r_{t,j/D}'$

$RC_t$ is positive definite as long as $N < D$ (portfolio dimension smaller than the number of intraday observations per day). For large $N$ this fails; the base-asset factor approach (above) resolves this. Non-synchronous trading can bias cross-asset covariances; correction methods (analogous to Scholes-Williams 1977 for CAPM betas) add lead and lag cross-products.

The log-normal/normal mixture model connects realized and GARCH-based approaches: if $\log(\sigma_t) | \mathcal{F}_{t-1} \sim \mathcal{N}(\mu,\tau^2)$ and $r_t|\sigma_t \sim \mathcal{N}(0,\sigma_t^2)$ , then daily returns follow a scale mixture of normals. ABDL show this fits FX returns well in terms of VaR coverage rates. See Value at Risk for the conditional distribution modeling context.

News Impact Curve

The news impact curve (Engle and Ng 1993) quantifies how the conditional variance at time $t$ responds to the single most recent innovation $\varepsilon_{t-1}$ , holding all prior information fixed at its unconditional level:

$\mathrm{NIC}(\varepsilon_{t-1}) = h_t\bigl|_{\text{lagged } h = \bar\sigma^2}$

For standard GARCH(1,1): $\mathrm{NIC}(\varepsilon) = A + \alpha_1 \varepsilon^2$ — symmetric in $\varepsilon$ , opening upward as a parabola. For EGARCH (Nelson 1991): $\ln h_t = \omega + \beta \ln h_{t-1} + \gamma z_{t-1} + \alpha(|z_{t-1}| - E|z|)$ — the NIC is asymmetric: negative shocks ( $z_{t-1}<0$ ) increase variance more than positive shocks of the same magnitude (leverage effect). For power-GARCH (Ding, Granger, and Engle 1993): $h_t^\delta = \omega + \alpha|\varepsilon_{t-1}|^{2\delta} + \beta h_{t-1}^\delta$ — estimated $\hat\delta = 1.43$ on daily S&P 500 (1928–1991), rejecting both $\delta=1$ and $\delta=2$ .

Engle and Ng (1993) also provide misspecification tests for GARCH models based on whether the standardized squared residuals correlate with sign and magnitude of past innovations.

Hentschel (1995) Family of GARCH Models

Hentschel (1995) nests eight widely-used GARCH specifications under a single variance equation governed by four parameters: a Box-Cox power $\lambda$ applied to the conditional standard deviation, a shock exponent $\nu$ , a shift $b$ , and a rotation $c$ :

$\frac{\sigma_t^\lambda - 1}{\lambda} = \omega + \alpha \sigma_{t-1}^{\lambda-\nu} f^\nu(\varepsilon_t) + \beta\frac{\sigma_{t-1}^\lambda - 1}{\lambda}, \qquad f(\varepsilon_t) = |\varepsilon_t - b| - c(\varepsilon_t - b) \tag{Hentschel}$

The function $f$ is a shifted and rotated absolute value. When $\lambda \to 0$ the left side reduces to $\ln \sigma_t$ (the EGARCH parameterization).

Nesting table — standard models as parameter restrictions on $(\lambda, \nu, b, c)$ :

Model	$\lambda$	$\nu$	$b$	$c$
EGARCH	0	1	0	free, $\|c\|<1$
TGARCH	1	1	0	free
AGARCH	1	1	free	0
GARCH	2	2	0	0
NA-GARCH	2	2	free	0
GJR-GARCH	2	2	0	free
NARCH	$\nu$	$\nu$	0	0
A-PARCH	$\nu$	$\nu$	0	free

Two types of asymmetry. The shift $b$ moves the news impact curve minimum rightward, dominating the asymmetric response for small shocks. The rotation $c$ changes the slopes on either side of the minimum, dominating for large shocks. Prior GARCH models allow only one type at a time, conflating two empirically distinct dimensions of the news impact curve.

Empirical result (CRSP 1926–1990, $T = 17{,}486$ daily excess returns). All symmetric and standard asymmetric models are rejected by LR statistics exceeding 200. The unconstrained optimum is $\hat\lambda \approx 1$ , $\hat\nu \approx 1.5$ — intermediate between the Taylor-Schwert absolute-value GARCH and Bollerslev squared-residual GARCH, closer to the AGARCH family. The shift $b$ is significant in every maintained specification; the rotation $c$ is not. Asymmetry in U.S. equity volatility is primarily a small-shock phenomenon: the leverage effect conventionally attributed to large negative returns is driven instead by the shift of the news impact curve minimum away from zero.

Mixture-of-Normals GARCH (MN-GARCH)

The standard GARCH(1,1) imposes a single conditional distribution on returns. The MN-GARCH model (Haas, Mittnik, and Paolella 2004a; Bauwens and Rombouts 2005) instead writes the conditional density as a $K$ -component normal mixture where each component has its own GARCH(1,1) variance process:

$f(y_t|\mathcal{F}_{t-1}) = \sum_{k=1}^K \pi_k\, \frac{\phi\!\bigl((y_t - \mu_k)/\sqrt{h_{k,t}}\bigr)}{\sqrt{h_{k,t}}}, \qquad h_{k,t} = \omega_k + \alpha_k y_{t-1}^2 + \beta_k h_{k,t-1}$

Zero-mean constraint: For the unconditional mean to equal zero, $\mu_K = -\sum_{k=1}^{K-1}(\pi_k \mu_k)/\pi_K$ .

Covariance stationarity: The mixture is weakly stationary even when individual components have $\alpha_k + \beta_k > 1$ (explosive), because the stationarity condition involves a weighted combination across all components. This permits the model to capture a heavy-tailed, high-variance regime (explosive component with small weight) alongside a persistent but stationary normal regime (dominant component).

Bayesian estimation (Bauwens-Rombouts 2005): A 4-block Gibbs sampler uses data augmentation via latent state indicators $S_t \in \{1,\ldots,K\}$ :

States $S_T$ : multinomial draw $P(S_t = k) \propto \pi_k \phi(y_t; \mu_k, h_{k,t})$
Weights $\pi$ : Dirichlet conjugate: $\text{Di}(a_{k0} + \#\{t: S_t=k\})$
Means $\mu$ : $(K-1)$ free means drawn jointly from a multivariate Normal enforcing the zero-mean constraint
GARCH parameters $(\omega_k, \alpha_k, \beta_k)$ : griddy-Gibbs over a 3-D grid (non-conjugate)

Marginal likelihood for $K$ selection: A Laplace approximation around the posterior mode $\hat\Psi$ gives $\log m(y) \approx \log L(\hat\Psi|y) + \log \phi(\hat\Psi) + \frac{\dim(\Psi)}{2}\log(2\pi) + \frac{1}{2}\log|\hat\Sigma|$ .

Near-IGARCH as misspecification artifact: Applied to S&P 500 returns (Jan 1994–Jun 2005, $T=3047$ ), $K=2$ is strongly preferred (marginal log-likelihood $-4091$ vs. $-4139$ for $K=1$ ). The two-component solution has: Component 1 weight $\approx 0.79$ , $\hat\alpha_1 + \hat\beta_1 \approx 0.98$ ; Component 2 weight $\approx 0.21$ , $\hat\alpha_2 + \hat\beta_2 \approx 1.14$ (explosive). The standard GARCH(1,1) estimate $\hat\alpha + \hat\beta \approx 0.996$ is a weighted compromise between the two components — the near-IGARCH finding is therefore a misspecification artifact from imposing a single distribution, not evidence of a unit root in variance. See Bauwens-Rombouts (2005).

Markov-Switching GARCH (Das-Yoo 2004)

Standard GARCH assumes a single volatility regime. The Markov Switching GARCH (MS-GARCH) model allows the level of the conditional variance to shift across a discrete set of regimes governed by a hidden Markov chain. Das and Yoo (2004) develop the first workable Bayesian Markov Chain Monte Carlo (MCMC) estimator for this model.

The path-dependence problem. In the model $\sigma_t^2 = \mu_0 + \mu_1 S_t + \sum_j \alpha_j \varepsilon_{t-j}^2 + \sum_j \beta_j \sigma_{t-j}^2$ with $S_t \in \{0,1\}$ Markov, the variance recursion implies $\sigma_t^2$ depends on the entire past history of states $S_1,\ldots,S_{t-1}$ . The exact likelihood therefore requires summing over $2^T$ state paths — infeasible for any realistic sample size. This is the fundamental obstacle to MLE for MS-GARCH; it also prevents the multi-move forward–backward state-sampling algorithm (used in simpler Markov switching models) from functioning correctly under GARCH dynamics.

Bayesian solution. Treating the full state vector $S = (S_1,\ldots,S_T)$ as a latent parameter and conditioning on it at every MCMC step eliminates path-dependence: given $S$ , the variance recursion is deterministic and the likelihood factorizes across $t$ as a product of Gaussian densities. This is the core idea that makes the Bayesian approach computationally tractable where MLE is not.

MCMC algorithm (four blocks):

State vector $S_t$ : Single-move sampling (one $S_t$ at a time). For each $t$ , the conditional probability $\Pr(S_t = 0 \mid Y, S_{\neq t}, \Theta)$ is computed by running the variance recursion forward from $t$ to $T$ , including all future contributions $\prod_{s=t}^T (\sigma_s^2)^{-1/2} \exp(-\varepsilon_s^2/2\sigma_s^2)$ ; a Bernoulli draw completes the update.
Transition probabilities $(p_{00}, p_{11})$ : Conjugate Beta posteriors using the state transition count statistics extracted from the current draw of $S$ .
Regression coefficients $\gamma$ : Metropolis–Hastings with a tailored Normal proposal from generalized least squares (GLS).
GARCH parameters $(\mu_0, \mu_1, \alpha, \beta)$ : Metropolis–Hastings using Nakatsuma's (2000) ARMA linearization $\varepsilon_t^2 \approx \mu_0 + \mu_1 S_t + \sum(\alpha_j+\beta_j)\varepsilon_{t-j}^2 + w_t - \sum \beta_j w_{t-j}$ where $w_t = \varepsilon_t^2 - \sigma_t^2$ ; $\beta$ is additionally linearized via a first-order Taylor expansion around a nonlinear least squares (NLS) pilot estimate.

Relation to switching ARCH (Kaufmann–Frühwirth-Schnatter 2002). Earlier Bayesian work on regime-switching volatility was limited to switching ARCH (no lagged variance term $\beta$ ) precisely because multi-move state sampling fails under GARCH path-dependence. Das and Yoo demonstrate that single-move sampling sidesteps this limitation at the cost of slower Markov chain mixing for the state sequence.

Numerical results. MS-GARCH(1,1), $T=500$ : posterior means track true values well; MH acceptance rates 13–19%; the GARCH persistence parameter $\beta$ shows chain autocorrelation $\hat\rho \approx 0.98$ ; the regime indicator $\mu_1$ is imprecisely estimated (posterior mean 0.150 vs. true 0.050). See Das-Yoo (2004).

QML Estimation

Let $\varepsilon_t = y_t - \mu_t(\theta)$ where $\mu_t(\theta)$ is the conditional mean (e.g., the Vector Error Correction Model (VECM) component). The quasi log-likelihood is

$\ell(\theta) = -\frac{Tn}{2}\log(2\pi) - \frac{1}{2}\sum_{t=1}^T \bigl[\log|H_t(\theta)| + \varepsilon_t'H_t(\theta)^{-1}\varepsilon_t\bigr]$

maximized over $\theta$ by nonlinear numerical optimization. QML is consistent and asymptotically normal even when $z_t$ is non-Gaussian (Weiss 1986 univariate; Bollerslev and Wooldridge 1992; Tuncer 1996 and Bauwens-Vandeuren 1996 for BEKK), provided the conditional mean and variance are correctly specified.

Standard errors. Under non-normality, the usual inverse Hessian SE is invalid. Robust (sandwich) SE use $\hat\Sigma = \mathcal{H}^{-1}\mathcal{G}\mathcal{H}^{-1}$ where $\mathcal{H}$ is the Hessian and $\mathcal{G}$ is the outer-product-of-gradients matrix.

Value at Risk Application

GARCH-based VaR (see Value at Risk) uses the estimated conditional variance $h_t$ to compute a time-varying quantile. Under conditional normality, the 99% one-day VaR for a portfolio of value $W$ is $W \times z_{0.01} \sqrt{h_t}$ where $z_{0.01} \approx 2.326$ . Nobel Committee (2003) S&P 500 example: the 99% VaR on a $1M portfolio ranged from $12,400 (September 1995, minimum volatility) to $61,500 (July 2002, maximum volatility) — illustrating the practical importance of time-varying volatility modeling.

Why It Matters

Volatility modeling: GARCH captures the well-documented volatility clustering in asset returns and interest rates, improving risk measurement (VaR, Conditional Value at Risk (CVaR)) and option pricing.
Efficiency gains: Ignoring GARCH heteroscedasticity in the residuals of a VECM amounts to using inefficient GLS weights; incorporating GARCH into the likelihood directly improves coefficient efficiency.
Test power: Cointegration tests (Johansen LR) under homoscedastic assumptions can fail to detect long-run relationships when GARCH innovations amplify short-run noise. BDV (1998) demonstrate this concretely for UK and German interest rates.
Common volatility factors: Off-diagonal BEKK elements capture cross-market volatility spillovers (e.g., between short and long interest rates), which are economically meaningful for yield curve dynamics and monetary policy transmission.

Open Questions

Distinguishing true IGARCH from near-unit-root variance in finite samples is difficult; formal tests for an exact unit root in variance (analogous to Dickey-Fuller for means) have low power.
Spurious IGARCH from structural breaks. Hillebrand (2004) shows that unaccounted-for parametric regime changes in GARCH models cause estimated $\alpha+\beta$ to converge to one even when the true process is not IGARCH — the direct variance analogue of Perron's (1989) result that ignoring structural breaks in the mean generates spurious unit roots. Any finding of IGARCH should therefore be preceded by a structural break test on the GARCH parameters; see Structural Break Testing.
The positive-definiteness guarantee of BEKK comes at the cost of a non-convex parameter space, making numerical optimization sensitive to starting values.
Bayesian inference for BEKK (e.g., via Gibbs sampling) requires specialized MCMC; see Bauwens and Lubrano (1994b) for an approach on the univariate GARCH model.