Parsimonious and Flexible Parameterizations for Realized Covariance Matrix Models

Introduction

The availability of detailed financial time-series data, including intraday observations recorded at millisecond frequencies, has fostered the development of high-quality measures of realized variances and covariances [1]. These measures are essential for assessing financial market risk and analyzing correlations among assets. Recent approaches have focused on modeling $(n \times n)$ realized covariance matrices, including matrices of relatively high dimension.

Two main technical challenges arise when modeling covariance matrices (see, e.g., Bauwens, Laurent, and Rombouts, 2006): (1) the curse of dimensionality, since the number of unknown parameters increases rapidly with $n$, given that a covariance matrix has $n(n + 1)/2$ distinct elements; and (2) the need to ensure that the model produces positive definite (PD) matrices.

The models proposed in the literature generally extend univariate conditional-variance models, such as GARCH models [10], using parameterizations designed to address these two challenges. Examples include scalar and diagonal specifications: in scalar models, all elements of the covariance matrix share the same set of coefficients, whereas in diagonal models, each element follows its own dynamics independently of the others, allowing equation-by-equation estimation. These approaches are commonly adopted in BEKK models [20]. However, they offer limited flexibility and do not account for potential interactions among the elements of the covariance matrix. Another widely used approach is the dynamic conditional correlation (DCC) model of Engle [18], which involves two estimation stages. In the first stage, conditional variances are estimated using $n$ univariate models; in the second, a scalar GARCH-type model is fitted to the correlation matrix.

More recent approaches have been developed specifically for realized measures. One example is the vech-HAR model of Chiriac and Voev [14], which extends the HAR model for realized volatility proposed by Corsi [17], and has a parametric structure similar to that of the diagonal BEKK model. Another general framework for realized covariance matrices is the conditional autoregressive Wishart (CAW) model proposed by Golosnoy, Gribisch, and Liesenfeld [23]. This framework uses the Wishart distribution as the conditional distribution of realized covariance matrices and exploits its properties to obtain consistent estimates. More recent contributions by Bauwens and Otranto [8, 9] extend the approach of Bauwens and Otranto [7] by introducing the Hadamard exponential parameterization. This parameterization provides greater flexibility by allowing the coefficient associated with each element of the covariance matrix to depend on a small number of unknown parameters, thereby substantially reducing the number of parameters to be estimated. An alternative approach, closer to the conventional normal-distribution framework, uses a matrix-variate log-normal distribution (see, e.g., Gupta and Nagar [24]). However, the resulting covariance matrix estimator is biased and requires a correction, based either on the ratio of the medians of the actual and estimated standard deviations [4] or on a second-order Taylor expansion [35].

In this paper, we develop a general family of models within the CAW framework, based on several alternative parameterizations. The proposed parameterizations ensure the positive definiteness of the estimated covariance matrices. Moreover, the coefficient matrices are time-varying, with distinct coefficients for each element of the covariance matrix, while depending on only a small number of unknown parameters. This parsimonious structure mitigates the curse of dimensionality and makes the models suitable for applications involving relatively high-dimensional matrices. The proposed parameterizations are based on multivariate logistic functions, multivariate linear functions, Hadamard exponential functions, and combinations of these approaches. We impose parameter constraints to prevent explosive dynamics and ensure the positive definiteness of the resulting conditional covariance matrices. We illustrate the approach within the BEKK framework, although it can readily be extended to alternative models by appropriately parameterizing the coefficient matrix associated with the lagged realized covariance matrix. We conduct Monte Carlo experiments to investigate the finite-sample properties of the estimators and their behavior as the time-series length increases. Finally, an empirical application shows that combining the proposed parameterizations can improve model performance. The remainder of the paper is organized as follows. § 2 introduces the CAW framework and the proposed family of models. § 3 presents the feasible parameterizations. § 4 investigates the properties of the estimators through Monte Carlo simulations. § 5 provides an empirical application using the 30 assets included in the Dow Jones Industrial Average (DJIA). § 6 concludes.

The Model Framework

Consider a multivariate time series of $n \times n$ realized covariance matrices, denoted by $\Gamma_t$, where $t$ is the time index. Let $\mathcal{F}_t = (\Gamma_t, \Gamma_{t-1}, \ldots)$ denote the information set at time $t$. Conditional on $\mathcal{F}_{t-1}$, $\Gamma_t$ follows a Wishart distribution with $\nu> n$ degrees of freedom and scale matrix $S_t/\nu$, where $S_t$ is symmetric and positive definite. Under this parameterization, $E(\Gamma_t \mid\mathcal{F}_{t-1}) = S_t$. In what follows, we model the conditional expectation $S_t$ using a BEKK representation [20] with time-varying coefficients, although the proposed approach can readily be extended to other model specifications. The general model is given by:

where $\odot$ is the symbol for the element-by-element (Hadamard) product, $K$, $A_t = \{a_{ij,t}\}$, $B_t = \{b_{ij,t}\}$ are $n \times n$ symmetric matrices of coefficients. The curse of dimensionality is apparent from the fact that these matrices involve $3n(n + 1)/2$ different parameters, of which $n(n + 1)$ can vary at any time $t$.

The number of coefficients can be reduced by replacing the matrix $B_t$ with a scalar $0 \le b < 1$; several works show empirical evidence of this result (see, for example, Bauwens and Otranto, 2025).¹

A further reduction in the number of unknown coefficients can be achieved through the covariance targeting specification; this is a simple extension of Engle and Mezrich [21]’s variance targeting idea to the multivariate case [see also 33]. It consists in reparameterizing the matrix $K$ so that the unconditional expectation of $S_t$ is equal to the unconditional expectation of $\Gamma_t$, which we will call $\bar{\Gamma}$. For model (1), this could be achieved by replacing $K$ with

where $J_n$ is an $n \times n$ matrix of ones. This substitution has two drawbacks: first, this parameterization does not provide a constant matrix $K$; second, and most importantly, (2) is not guaranteed to be PD. One possibility is to adopt the Hafner and Franses [25] approximation, replacing $A_t$ in (2) with a scalar obtained from the average of all its elements at each instant. We propose an alternative parameterization by introducing a new unknown parameter $\kappa$ and setting $K = \kappa\bar{\Gamma}$, with $\kappa> 0$. As in Hafner and Franses [25], it does not provide exact covariance targeting because it implies that $E(S_t) \ne\bar{\Gamma}$, but adds flexibility to the constant part of the model. The unconditional mean $\bar{\Gamma}$ can be consistently estimated in a preliminary step as the sample mean of the realized covariance matrices.

Considering the new parameterization of model (2.1):

the following proposition holds:

Proof. Obviously, from (3) $S_t$ is PSD and, by definition, all the $\Gamma_t$ covariance matrices are PSD. By the Schur product theorem, the Hadamard product of a PD matrix and a PSD matrix results in a PSD matrix, so $A_t \odot\Gamma_{t-1}$ is PSD. If at least one covariance matrix $\Gamma_t$ is PD, then the sample mean $\hat{\bar{\Gamma}}$ is PD, and also $\kappa\hat{\bar{\Gamma}}$, being $\kappa> 0$. It follows that the right-hand side of (3) is a sum of a PD matrix and two PSD matrices, thus yielding a PD matrix $S_t$, with different coefficients for all the variances and covariances, but with a small set of unknown coefficients.

The presence of time-varying coefficients implies that (probably) the first and second moments of the process (3) are not constant; however, in the framework of time-varying coefficients, non-stationarity can be modeled (see Lütkepohl, [31]). What we need is to derive some (uniform contraction) condition (Kawano and Hosoe, [30]) to rule out explosive dynamics arising from the linear recursion.

Proof. Applying the vech operator to (3), which stacks the lower triangular elements of a symmetric matrix into a vector of dimension $n^{*} = n(n + 1)/2$, we obtain

where $s_t$, $\bar{\gamma}$, and $\gamma_t$ denote the vectors obtained by applying the vech operator to the corresponding matrices denoted with capital letters, and $D_{A_t} = \operatorname{diag}(\operatorname{vech}(A_t))$, where $\operatorname{diag}(\bullet)$ is the operator applied to a vector to create a diagonal matrix with diagonal entries equal to the elements of the same vector.

Adding and subtracting $\gamma_t$ and $b\gamma_{t-1}$ on the right-hand side of (4), and using $\eta_t = \gamma_t - s_t$, yields the following VARMA-type representation with time-varying autoregressive coefficients:

where $I_{n^{*}}$ denotes the identity matrix of dimension $n(n + 1)/2$.

Let $D_t = D_{A_t} + bI_{n_*}$. Since $D_t$ is diagonal, its diagonal entries are $a_{ij,t} + b$. Hence, the matrix norm induced by the vector infinity norm satisfies

Suppose that $\sup_t \lVert D_t\rVert_{\infty} \leq\rho< 1$. By repeated substitution, the coefficient matrix multiplying $\gamma_{t-h-1}$ in the associated homogeneous recursion is $\prod_{k=0}^{h} D_{t-k}$. By submultiplicativity of the induced matrix norm,

Therefore, the homogeneous component of the associated time-varying linear recursion decays geometrically, establishing uniform contraction. By (2.6), this is verified when $\sup_t \max_{i,j} \left\lvert a_{ij,t} + b\right\rvert< 1$.

Finally, since the moving-average polynomial in (2.5) is $(I_{n_*} - bI_{n_*}L)$, where $L$ denotes the lag operator, the innovation part is invertible when $\lvert b\rvert< 1$.

In the next section, we will provide some different feasible parameterizations for the matrix $A_t$ that guarantee their PSDness and parsimony. Furthermore, to ensure the positivity of the variances, we impose $0 \leq b < 1$.

Feasible Parameterizations

We propose four parameterizations of the coefficient matrix $A_t$ that ensure its positive semidefiniteness while allowing all its elements to vary over time through a small number of parameters. The proposed specifications are based on logistic, linear, and Hadamard exponential (HE) functions. positive semidefiniteness of $A_t$ implies that its diagonal elements are nonnegative, the resulting conditional covariance matrix $S_t$ has nonnegative diagonal elements.

Logistic Parameterization

In (2.3), $A_t$ is reparameterized as follows:

where $a_0 > 0$, $a_1 \geq0$, and $(a_0 + a_1 + b) < 1$. $G_t$ is a multivariate logistic function, defined as:

with $\phi> 0$ and $c$ a scalar coefficient, $I_n$ denotes the $n$-dimensional identity matrix. Logistic functions have been used to obtain time-varying coefficients by Teräsvirta [36], introducing smooth-transition autoregressive models in a univariate context. The coefficient $c$ acts as a threshold value, while $\phi$ influences the shape of the logistic function, which can vary from a constant straight line to an S-shape to a step function. Potential identification problems arise because, as $\phi$ increases, and after a certain point, the shape of the function does not change (logistic function saturation). We discuss this issue in greater detail in § 4.

The parameterization (3.1) ensures the positive definiteness of $S_t$ in (2.3) and, under the constraints on the coefficients, the stability and invertibility of the process.

Proof. Let $0_n$ denote the $n \times n$ matrix of zeros. Considering the spectral decomposition of $\Gamma_t$:

with $\Lambda$ a diagonal matrix containing the $n$ eigenvalues $\lambda_i \ge0$ and $V_t$ the corresponding matrix of eigenvectors, for the spectral definition of matrix functions (see, for example, [28]) we obtain the spectral decomposition of $G_t$:

with $\Lambda^{*}$ a diagonal matrix with elements $0 < (1 + \exp(-\phi(\lambda_i - c)))^{-1} < 1$, so that $0_n \prec G_{X,t} \prec I_n$, where the symbol $\prec$ denotes the Loewner order, so that $X \prec Y$ if and only if $Y-X$ is positive definite.² This implies the positive definiteness of $G_t$ and, as a consequence, $0_n \prec A_t \preceq a_0 + a_1 I_n$.

Proof. Let $a_0 + g_{ij,t}$ be the $(i,j)$ entry of $A_t$ in (3.1). From Proposition 2, the eigenvalues of the matrix $D_t$ are $(a_0 + g_{ij,t} + b)$.

Proposition 3 establishes that $0_n \prec A_t \preceq a_0 + a_1 I_n$; this implies that $|g_{ij,t}| \le a_1$ for each $(i,j)$. Thus, uniform contraction is verified if every element on the diagonal of $D_t$ is less than, in magnitude, $(a_0 + a_1 + b) < 1$ for every $t$. Due to the non-negativity constraints on the coefficients $a_0$, $a_1$, $b$, the previous result also implies $0 \le b < 1$, which is the invertibility condition derived in Proposition 2.

We call model (2.3) with parameterization (3.1) the Log-BEKK model.

Linear Parameterization

An alternative (simpler) specification can be obtained by parameterizing matrix $A_t$ in (2.3) as:

where $C$ is a PD matrix with all diagonal elements (variance coefficients) equal to $\alpha_v$ and all off-diagonal elements (covariance coefficients) equal to $\alpha_c$ ($\alpha_v > \alpha_c \ge0$; $\alpha_a \ge0$).

This parameterization introduces a dependence of $S_t$ on the elementwise product of the lagged realized covariance matrix and a coefficient matrix that varies with $\Gamma_{t-1}$.

The coefficient $\alpha_a$ captures a negative effect on $S_t$ and is consistent with the idea of an attenuation effect proposed by Bollerslev, Patton, and Quaedvlieg [11] in the univariate case, and subsequently extended to the multivariate case by Bollerslev, Patton, and Quaedvlieg [12].³ In these works, the authors parameterize the lagged realized variance coefficients as a linear function of the lagged realized quarticity. Cipollini, Gallo, and Otranto [16] demonstrates that realized variance is a valid substitute for quarticity, with the advantage of being more readily available. Therefore, the parameterization (9) can be considered an extension to the multivariate case of the univariate parameterization proposed by Cipollini et al. [16], reducing the effect of possible measurement errors and providing a time-varying matrix of coefficients $A_t$, all different from each other.

The following propositions establish conditions under which $A_t$ is positive semidefinite.

Proof. From (3.4) we immediately note that the largest eigenvalue of $C$ is equal to $\alpha_v + (n - 1)\alpha_c$ and the remaining $n - 1$ eigenvalues are all equal to $\alpha_v - \alpha_c$. $A_t$ is PSD if $C \succeq\alpha_a\Gamma_{t-1}$, that is, if the eigenvalues of this difference are non negative for every $t$. A sufficient condition is that $\alpha_a \max(\xi) \leq(\alpha_v - \alpha_c)$, yielding the constraint for $\alpha_a$ expressed in Proposition 6.

Proof. Denoting $\gamma_{rs,t}$ as the generic element of the covariance matrix $\Gamma_t$, the parameterization (3.3) implies that every diagonal element of $D_t$ in Proposition 2 is given by $(\alpha_v - \alpha_a\gamma_{rr,t} + b)$ (when referring to the diagonal elements of $A_t$) and $(\alpha_c - \alpha_a\gamma_{rs,t} + b)$ (when referring to the off-diagonal elements of $A_t$); $r, s = 1, \ldots, n$; $r \ne s$. Considering that $\alpha_a\gamma_{rr,t} \geq0$

For the off-diagonal elements $|\gamma_{rs,t}| < \max(\xi)$, and

The condition derived in Proposition 2 is therefore satisfied by $(\alpha_v + b) < 1$.

The parameter constraints for the linear specification are computed using the maximum eigenvalue observed in the estimation sample. They therefore ensure positive semidefiniteness of the coefficient matrices over that sample, but do not by themselves guarantee the same property for future observations.

We call model (2.3) with parameterization (3.3) the Lin-BEKK model.

Hadamard Exponential Parameterization

Bauwens and Otranto [8] introduce the following parameterization of $A_t$ for modeling realized covariance matrices based on the HE operator:

where $\exp^{\odot}$ represents the HE operator (element-by-element exponential), $P_t$ is a correlation matrix. The coefficient $a$ is strictly positive, and $\varphi$ is nonnegative ($\varphi= 0$ provides a constant coefficient $a$ for all variances and covariances).

The parameterization (3.6) implies that the diagonal elements of $A_t$ are equal, while the off-diagonal elements are different. In practice, the variance coefficient is constant and equal for all, while the covariance coefficients are specific and time-varying, increasing when there is a greater correlation, at the previous time, between the corresponding assets. Bauwens and Otranto [7] prove that $W_t$ is PD. Furthermore, considering that all the elements of $W_t$ are nonnegative and that the highest value is $a$, a sufficient condition for uniform contraction of the homogeneous component of the associated time-varying linear recursion is $(a + b) < 1$.

We call model (2.3) with parameterization (3.6) the HE-BEKK model.

To allow for variance coefficients that are not constant and equal across assets, we can enrich the parametrization (3.6), for the variance part, with the logistic and linear parameterizations previously described.

By combining the Log-BEKK for the variance coefficients and the HE-BEKK for the covariance coefficients, we obtain the Log-HE-BEKK model:

where $V_{t-1}$ is a diagonal matrix containing the realized variances at time $t-1$. Both $W_t$ and $G_t^*$ are PD, $a_1$ is nonnegative, which implies positive definiteness of $A_t$ in (13). Based on the previous results on the stability of the Log-HE model (Proposition 4), it is trivial to show that a sufficient condition for uniform contraction of the homogeneous component of the associated time-varying linear recursion is $(a+a_1+b)<1$.

By combining the Lin-BEKK for the variance coefficients and the HE-BEKK for the covariance coefficients, we obtain the Lin-HE-BEKK model: (1):

with $\alpha_a \ge0$. Again, $\alpha_a$ represents the attenuation effect described in subsection 3.2 (in this case referring only to the variances and not to the covariances). The parameterization (14) can be considered a special case belonging to the family of models developed in Bauwens and Otranto [9].

Proof. Both $W_t$ and $V_{t-1}$ are PD matrices. We rewrite $A_t$ in (3.6) as:

$A_t$ is PSD if and only if

Consequently, a sufficient and necessary condition for PSDness of $A_t$ is that $\alpha_a$ must be less than or equal to the minimum eigenvalue of $V_{t-1}^{-1/2}W_tV_{t-1}^{-1/2}$ for every $t$, i.e. $\alpha_a \le\min(\xi^*)$. Considering that $W_t$ depends on the two unknown coefficients $a$ and $\varphi$, this condition is difficult to enforce in the estimation routine. We can rewrite:

Considering that:

where $\odot h$ represents the $h$-th Hadamard (element-wise) power. By Schur’s theorem, every element on the right of (16) is PSD and

Consequently:

From (15), we can derive as a sufficient condition for the PSDness of $A_t$:

so that the constraint guaranteeing the condition is:

Proof. The diagonal elements of $A_t$ are $a-\alpha_a\gamma_{rr}<a$, while the off-diagonal elements are less than $a$. Consequently every element of the diagonal matrix $D_t$ in Proposition 2 is not greater than $(a+b)$, so uniform contraction is obtained if $(a+b)<1$.

Estimator Properties

The estimation of model (2.3) carried out using a two-step M-estimation procedure. In the first step, the matrix $\bar{\Gamma}$ is estimated; in the second step, all remaining parameters are estimated conditional on the first-step estimate. Under standard regularity conditions, replacing nuisance quantities with consistent first-step estimators preserves the consistency of the second-step estimator [32, 38]. The matrix $\bar{\Gamma}$ can be consistently estimated by the sample mean of the realized covariance matrices. In the second step, $\kappa$, $b$, and the reduced set of parameters defining $A_t$ are estimated by quasi maximum likelihood (QML).

Asai and So [2] developed QML estimation for CAW-BEKK models and established the consistency and asymptotic normality of the resulting estimators under suitable regularity conditions. See also Asai and So [3] for the covariance-targeting case.

Let us decompose $\Gamma_t$ as follows:

where $S_t$ is parameterized as in (2.3). Under the following assumptions:

1. $S_t$ is $\mathcal{F}_t$-measurable;

2. $E(\mathcal{F}_t \mid\mathcal{F}_{t-1}) = S_t$

3. $U_t$ is independent on $\mathcal{F}_{t-1}$ and $U_t \sim W(\nu^{-1}I_n,\nu)$

the second-step estimation can be performed by maximizing the log-likelihood:

with

The log-likelihood (18) is obtained without specifying the conditional distribution of the covariance matrices and can be interpreted as QML estimator. Furthermore, the estimates of the model parameters do not depend on $\nu$, so we can exclude $\nu$ from the estimation procedure. This will not affect the estimates of the other parameters, which are consistent.⁴ This result is also supported by the simulation results in Asai and So [2], who found that, setting different values for $\nu$, the estimates of the other parameters do not change.

The properties established for the CAW-BEKK model cannot be readily extended to the proposed models. The two-step procedure may entail some loss of efficiency, and the process need not be stationary. Moreover, the asymptotic properties of the estimators require further investigation because the model combines approximate covariance targeting, based on a first-step estimate, with nonlinear specifications of $A_t$. A formal analytical study of these properties is beyond the scope of this paper. Instead, we conduct Monte Carlo experiments to provide empirical evidence on the estimators’ behavior in finite samples and as the sample size increases.

Monte Carlo Evidence

We conducted a Monte Carlo experiment with a data generating process represented by (3) with $n = 30$, using the four alternative parameterizations for $A_t$. Additionally, we generated sets of covariance matrices constructed, as a benchmark, from a simple scalar BEKK model (sc-BEKK):

and from the HE-BEKK (12) of Bauwens and Otranto [8]. The data generating parameters are the estimated values, reported in Table 2, for the application described in § 5. For each specification, we generated 500 replications for each sample size, $T = 100, 500, 1000, 2000$. The covariance matrices $\Gamma_t$ were generated from the Wishart distribution with parameters $S_t$ and $\nu= 90$;⁵ The initial value of $S_0$ was extracted from the realized covariance series of the application in § 5.

Parameters	sc	Log	Lin	HE	Log-HE	Lin-HE
$\kappa$	0.089	0.144	0.111	0.109	0.131	0.109
	(0.018)	(0.028)	(0.021)	(0.019)	(0.026)	(0.020)
$a$	0.381	0.312		0.492	0.419	0.492
	(0.051)	(0.024)		(0.068)	(0.040)	(0.069)
$b$	0.554	0.355	0.381	0.362	0.350	0.362
	(0.055)	(0.098)	(0.094)	(0.086)	(0.085)	(0.087)
$a_1$		0.163			0.069
		(0.049)			(0.026)
$\phi$		10.318			15.062
		(3.662)			(5.217)
$c$		0.780			1.994
		(0.078)			(0.150)
$\varphi$				0.763	0.491	0.764
				(0.117)	(0.044)	(0.116)
$\alpha_v$			0.477
			(0.074)
$\alpha_c$			0.338
			(0.031)
$\alpha_a$			1.439	2.880		0.128
			(0.494)	(0.871)		(0.083)
Log-Likelihood	-89064.6	-88490.9	-88578.9	-88527.1	-88480.4	-88526.8
AIC	64.45	64.04	64.10	64.06	64.03	64.06
BIC	64.46	64.05	64.11	64.07	64.04	64.07

Table 1 reports the bias, root mean squared error (RMSE), and the smallest sample size for which normality is not rejected (based on a Jarque-Bera test at 1% significance). The results are consistent with the consistency of the estimators across all models, with a marked reduction in bias as $T$ increases; in particular, the bias of the $\kappa$ estimator converges to zero very quickly for all models. For $T = 2000$ the maximum bias is 0.008 ($\alpha_a$ in the Lin-BEKK model). For small samples ($T = 100$, 500),⁶ the bias is, in absolute terms, less than 0.06, excluding the coefficients multiplied by the covariance matrices constructed to give the time variability of $A_t$: $\phi$ in Log-BEKK and Log-HE-BEKK, and $\alpha_a$ in Lin-BEKK.

Parameter	True	$T = 100$ Bias	$T = 100$ RMSE	$T = 500$ Bias	$T = 500$ RMSE	$T = 1000$ Bias	$T = 1000$ RMSE	$T = 2000$ Bias	$T = 2000$ RMSE	Normality from $T=$
sc-BEKK
$\kappa$	0.089	-0.002	0.003	-0.001	0.001	0.000	0.001	0.000	0.000	500
$a$	0.381	0.027	0.029	0.005	0.006	0.003	0.003	0.001	0.002	500
$b$	0.554	-0.024	0.026	-0.004	0.005	-0.002	0.003	-0.001	0.002	100
Log-BEKK
$\kappa$	0.144	0.005	0.007	0.001	0.002	0.000	0.002	0.000	0.001	500
$a_0$	0.312	0.019	0.021	0.004	0.005	0.002	0.003	0.001	0.002	500
$a_1$	0.163	0.023	0.028	0.005	0.006	0.002	0.003	0.001	0.001	500
$\phi$	10.318	0.663	10.928	-0.056	1.414	-0.002	0.797	0.001	0.007
$c$	0.780	0.011	0.048	0.001	0.016	0.000	0.011	0.000	0.006	500
$b$	0.355	-0.050	0.058	-0.010	0.012	-0.005	0.006	-0.002	0.002	500
Lin-BEKK
$\kappa$	0.111	0.002	0.419	0.000	0.168	0.000	0.113	0.000	0.047	100
$\alpha_v$	0.477	0.039	4.465	0.008	0.943	0.004	0.498	0.002	0.262	500
$\alpha_c$	0.338	0.022	2.366	0.004	0.516	0.002	0.274	0.001	0.169	500
$\alpha_a$	1.439	0.174	22.91	0.036	4.921	0.019	2.701	0.008	1.317	500
$b$	0.381	-0.044	5.131	-0.009	1.064	-0.004	0.545	-0.002	0.253	500
HE-BEKK
$\kappa$	0.109	0.006	0.008	0.001	0.002	0.001	0.002	0.000	0.001	100
$a$	0.492	0.048	0.053	0.010	0.012	0.005	0.006	0.002	0.003	500
$\varphi$	0.763	0.024	0.041	0.005	0.011	0.003	0.007	0.001	0.004	500
$b$	0.365	-0.059	0.066	-0.013	0.015	-0.006	0.008	-0.003	0.004	500
Log-HE-BEKK
$\kappa$	0.131	0.006	0.009	0.001	0.002	0.001	0.002	0.000	0.001	100
$a$	0.419	0.028	0.036	0.007	0.009	0.004	0.002	0.002	0.003	500
$\varphi$	0.491	-0.002	0.080	0.003	0.018	0.002	0.010	0.001	0.006
$a_1$	0.069	0.016	0.040	0.002	0.007	0.001	0.003	0.000	0.002	500
$\phi$	15.062	240.8	1271.1	0.570	9.472	-0.006	0.007	0.008	0.065
$c$	1.994	-0.023	0.213	0.003	0.034	0.003	0.021	0.001	0.015	500
$b$	0.350	-0.054	0.062	-0.011	0.014	-0.006	0.007	-0.003	0.004	500
Lin-HE-BEKK
$\kappa$	0.109	0.006	0.008	0.001	0.002	0.001	0.001	0.000	0.001	100
$a$	0.492	0.048	0.054	0.010	0.012	0.005	0.006	0.002	0.003	1000
$\varphi$	0.764	0.025	0.040	0.005	0.011	0.002	0.007	0.001	0.004	500
$\alpha_a$	0.128	0.013	0.015	0.003	0.003	0.001	0.002	0.001	0.001	2000
$b$	0.362	0.013	0.015	-0.013	0.015	-0.006	0.008	-0.003	0.004	1000

The results are also encouraging with respect to efficiency; despite the two-stage estimation procedure, the parameter estimators of the new class of models show a similar reduction in RMSE compared to the benchmark scalar BEKK. The main exceptions are again $\phi$ in Log-BEKK and $\alpha_a$ in Lin-BEKK; the former maintains a high RMSE up to $T = 1000$, then decreases dramatically for $T = 2000$, while the latter shows the highest RMSE (equal to 1.317) for $T = 2000$. The coefficient $\phi$ in Log-HE-BEKK starts from a very high RMSE (240.8) for $T = 100$, but decreases very quickly, reaching the value of 0.065 for $T = 2000$. The Log-HE-BEKK model is the most flexible specification, featuring distinct nonlinear parameterizations for the variance and covariance parameters; this greater flexibility may come at the cost of increased estimation variability in small samples.

In general, normality is verified even for small samples; most estimators achieve normality starting from $T = 500$; this result is obtained at $T = 100$ for the estimator of $\kappa$ in Lin-BEKK, Log-HE-BEKK and Lin-HE-BEKK, as for the HE-BEKK, while, in Lin-HE-BEKK, for $a$ and $b$ starting from $T = 1000$ and for $\alpha_a$ from $T = 2000$. The estimator of $\phi$ in Log-BEKK and Log-HE-BEKK does not achieve normality, nor does $\varphi$ in Log-HE-BEKK. The p-value of the Jarque-Bera statistic of $\varphi$ increases with $T$ and is unlikely to reject normality for values of $T$ greater than 2000, while the p-value of $\phi$ is always zero for both Log-BEKK and Log-HE-BEKK. A plausible explanation is the magnitude of the estimated parameter. As stated in subsection 3.1, $\phi$ determines the shape of the logistic function, which does not change for large values of this parameter (saturation of the logistic matrix function). Figure 1 shows the likelihood profile with respect to $\phi$, computed on the real 30-variable dataset that will be analyzed in § 5 (DJIA dataset). It is clear that, for values of $\phi$ approximately greater than 7, the likelihood profile is flat, resulting in weak local identification.

Empirical Application

We consider daily time series of annualized covariance matrices constructed from the 30 stocks included in the Dow Jones Industrial Average (DJIA), covering the period from January 2, 2015, to May 29, 2026. We use observations through the end of 2025 ($T = 2765$) for in-sample estimation and model comparison, reserving the 102 observations from 2026 for out-of-sample evaluation.⁷ The latter is based on one-step-ahead predictions, with parameter estimates held fixed at their in-sample values. The data are freely available from the VOLARE platform (volare.unime.it; Cipollini, Cruciani, Gallo, Insana, Otranto, and Spagnolo, 2026).

We estimate six BEKK models: sc-BEKK, Log-BEKK, Lin-BEKK, HE-BEKK, Log-HE-BEKK, and Lin-HE-BEKK. The estimates are reported in Table 2. We use QML estimation and calculate standard errors using the robust sandwich covariance estimator proposed by White [37].

In terms of log-likelihood, the proposed models outperform sc-BEKK. Log-HE-BEKK achieves the highest log-likelihood, which includes more parameters. However, the AIC and BIC results likewise indicate that Log-HE-BEKK performs well relative to the alternative specifications.

All parameter estimates are statistically significant at conventional levels, except for $\alpha_a$ in Lin-HE-BEKK, which is significant at the 10% level. The estimated impact of new information—represented by the scalar $a$ in sc-BEKK and by the matrix $A_t$ in the other models—differs substantially between the scalar and time-varying specifications. The persistence coefficient $b$ is notably larger in sc-BEKK, indicating a greater contribution from the lagged conditional covariance matrix $S_{t-1}$ relative to the lagged realized covariance matrix $\Gamma_{t-1}$, It is useful to compare the coefficient $a$ in Log-BEKK and Log-

HE-BEKK. In the former, it represents the minimum value of the coefficient of lagged realized variances and covariances, augmented by the values of the multiple logistic function; in the latter, it represents the maximum value (in absolute terms) of the correlation coefficients, corrected (downward) by the HE function of lagged correlations, and the minimum value of the coefficient of realized variances, corrected (upward) by the values of lagged realized variance via the logistic function.

The dynamics of the time-varying coefficients (for the lagged realized variance of the Apple Inc.-AAPL stock and the lagged covariance between AAPL and Amgen Inc.-AMGN), for the Log-BEKK, Log-HE-BEKK, HE-BEKK and the constant coefficient $a$ of the sc-BEKK are illustrated in Figure 2. HE-BEKK and sc-BEKK have constant variance coefficients, but the former is specific to variances, while the latter is common to variances and covariances. The news effect is smaller for covariances than for variances, so the specific coefficient of HE-BEKK is 0.11 points larger than that of sc-BEKK. Models with time-varying coefficients show similar movements, with $a$ representing the minimum value of the variance coefficients (0.312 for Log-BEKK and 0.419 for Log-HE-BEKK). For both models, the increase in coefficients depends on the logistic function, but only on the variances for Log-HE-BEKK and the variances and covariances for Log-BEKK; consequently, the highest values of the Log-HE coefficients are higher than those of Log (0.487 versus 0.475) and move within a narrow range. Interestingly, the highest values of the Log-BEKK and Log-HE-BEKK variance coefficients correspond to the highest volatility peaks (top left panel): the Chinese stock market crash of August 2015; the peak of November 20, 2018, following the resurgence of US-China trade tensions and growing concerns about global growth and monetary tightening; the effects of the 2020 COVID pandemic; and Trump’s ”Liberation Day” tariffs of April 2025.

In terms of covariance coefficients, the Log-HE and HE models exhibit very similar behavior, depending on the same HE structure. The Log-BEKK coefficients exhibit similar movements, but a small variability (ranging between 0.28 and 0.34), while the sc-BEKK model overestimates this coefficient.

The greater flexibility of the Log-HE-BEKK model explains its superior performance in terms of the likelihood function, AIC and BIC, as shown in Table 2. A further comparison is conducted using loss functions; we consider a statistical loss, the Quasi-Likelihood (QL),⁸ defined as:

and an economic loss function, the Global Minimum Variance (GMV, Engle and Colacito, [19]):

with weights $\hat{w}_t = \sqrt{n}\hat{S}_t^{-1}j_n/(j_n^{\prime}\hat{S}_t^{-1}j_n)$. $T^*$ represents the number of periods for which losses are calculated (referring to in-sample and out-of-sample periods).

Based on these statistics, for both the in-sample and out-of-sample spans, we select the models with the best forecasting performance using the Model Confidence Set (MCS) approach of Hansen, Lunde, and Nason [26].⁹ Of course, this exercise selects the models with the best fit in the in-sample analysis and the models with the best prediction performance for the out-of-sample case. Hansen et al. [26] suggest a test significance level of 10% or 25% (the latter to make a more restricted selection).

The comparison is enriched by two alternative approaches, which are very popular in applications: the vecHAR model (Chiriac and Voev, [14]), which can be estimated with the ordinary least squares method, and the exponentially weighted moving average (EWMA) model (J.P.Morgan, [29]), with a decay parameter equal to 0.96.

The two loss functions show similar results for both in-sample and out-of-sample periods. The Log-HE model achieves the best performance in terms of QL, but the Log, HE, and Lin-HE versions also fall into the best set at the 10% significance level in the in-sample case, and Lin-HE in the out-of-sample case. It is worth noting that at the 25%

	EWMA	vecHAR	sc	Log	Lin	HE	Log-HE	Lin-HE
QL	65.977	67.667	64.446	64.031*	64.095	64.057*	64.023	64.057*
GMV	31.730	32.000	32.175	24.974	24.957	23.281	24.095	23.278
	EWMA	vecHAR	sc	Log	Lin	HE	Log-HE	Lin-HE
QL	80.564	83.324	80.391	80.316	80.377	80.300	80.287	80.299*
GMV	18.866	23.647	23.592	19.737	18.879	21.474	19.020	18.010

significance level, only Log-HE is in the set of superior models, for both the in-sample and out-of-sample exercises.

Considering the GMV loss function, the Lin-HE version significantly outperforms all other models. Interestingly, in terms of AIC and BIC, Lin-HE does not improve on the HE-BEKK model (see Table 2), but the changes to the variance coefficient it captures imply a significant improvement in GMV, especially in out-of-sample analysis.

The new models show better performance than the classic EWMA, vecHAR, and sc-BEKK. This is an expected result for in-sample analysis, but is non-trivial for out-of-sample; the fact that Hansen and Timmerman (2015) and Xie and Yu (2020) found that the simpler models performed better for out-of-sample analysis further emphasizes the performance of Log-HE-BEKK and Lin-HE-BEKK.

Final Remarks

The availability of accurate measures of variances and covariances, such as realized covariance matrices, together with advances in computational power, has stimulated growing interest in covariance matrix models. When the objective is to model high-dimensional matrices, the curse of dimensionality compounds the challenge of ensuring that the model produces positive definite matrices.

We propose a family of parsimonious yet flexible models for realized covariance matrices. By estimating a relatively small number of parameters, the proposed models accommodate time-varying coefficients while ensuring the positive definiteness of the estimated conditional covariance matrices. The Monte Carlo results provide evidence on the finite-sample behavior of the estimators and their performance as the sample size increases, although some parameters exhibit non-normal sampling distributions when the likelihood is relatively flat. The empirical application indicates that the proposed models can outperform commonly used benchmark specifications.

Several extensions are possible. First, we develop the models within the BEKK framework, but the proposed approach can be extended to other baseline specifications, such as DCC [18], multivariate HAR [4], and DRD (Oh and Patton, 2016) models.

Second, we adopt a scalar specification for $B_t$ in (2.1), motivated by previous empirical evidence. However, extending the model to a non-scalar $B_t$ is straightforward. In this case, the sufficient condition for uniform contraction of the homogeneous component of the associated time-varying linear recursion in Proposition 2 becomes

where $b_{ij,t}$ denotes the $(i,j)$-th element of $B_t$.

Third, we use the standard exponential function in the Log-BEKK model and the Hadamard exponential function in the HE-type models. Either function could be used in both classes of models. We adopt the specifications that perform best in the empirical application.

Fourth, we use lagged realized covariance matrices—or lagged correlation matrices in the HE specifications—to drive the time variation in $A_t$. More generally, this role could be assigned to any symmetric, time-varying positive semidefinite matrix $M_t$. Furthermore, if $B_t$ is parameterized in the same way as $A_t$, the two coefficient matrices need not be driven by the same matrix $M_t$.

The proposed specifications can be viewed as multivariate covariance recursions with state-dependent coefficients. They are related to the literature on random-coefficient multivariate ARCH/GARCH models (e.g. Fong, Li, and An [22]) and, more generally, to stochastic recurrence equations with random coefficient matrices [13].