Time series models for realized covariance matrices based on the matrix-F distribution

Time series models for realized covariance matrices based on the matrix-F distribution

Introduction

Modeling the multivariate volatility of many asset returns is crucial for asset pricing, portfolio selection, and risk management. After the seminal work of Barndorff-Nielsen and Shephard (2002, 2004) and Andersen et al. (2003), the realized covariance (RCOV) matrix, estimated from the intra-day high frequency return data, has been recognized as a better estimator than the daily squared returns for daily volatility. Consequently, increasing attention has been focused on the modeling and forecasting of these RCOVs; see, e.g., McAleer and Medeiros (2008), Hansen et al. (2012), Noureldin et al. (2012), Bollerslev et al. (2016), and many others.

Existing models for the RCOV matrices can be roughly categorized into two types: transformation-based models and likelihood-based models. Models in the first category capture the dynamics of the RCOV matrices in an indirect way via transformation. Bauer and Vorkink (2011) used a factor model for the vectorization of the log transformation of RCOV matrix; Chiriac and Voev (2011) applied a vector autoregressive fractionally integrated moving average process to model the Cholesky decomposition of RCOV matrix; Callot et al. (2017) transformed the RCOV matrix into a large vector by the $vech$ operator, and then fitted this transformed vector by a vector autoregressive model. In the first two models, the dimension of RCOV matrix has to be moderate (e.g., less than 6) for a feasible manipulation. In the third model, the dimension of RCOV matrix is allowed to be 30 in applications with the help of the LASSO method.

Models in the second category deals with RCOV matrices directly by assuming that the innovation, which drives the RCOV time series, has a specific matrix distribution to generate random positive definite matrices automatically without imposing additional constraints. This important feature results in positive-definite estimated RCOV matrices. Unlike scalar or vector distributions, so far only few matrix distributions have been found to have explicit forms. The primary choice for the innovation distribution is Wishart, leading to the Wishart autoregressive (WAR) model in Gouri\'{e}roux et al. (2009), the conditional autoregressive Wishart (CAW) model in Golosnoy et al. (2012), the mixture Wishart model in Jin and Maheu (2013, 2016), and the generalized CAW model in Yu et al. (2017) to name a few. The other choice for the innovation distribution is matrix-F, which was recently adopted by Opschoor et al. (2018). Generally speaking, matrix-F distribution is the generalization of the usual F distribution, while Wishart distribution is the generalization of the $\chi^2$ distribution (see, e.g., Konno (1991) and Opschoor et al. (2018) for more discussions). Therefore, matrix-F distribution could be more appropriate than Wishart distribution in capturing the heavy-tailed innovation, which is an important stylized fact in many applications (see, e.g., Bollerslev (1987), Fan et al. (2014), Zhu and Li (2015), and Oh and Patton (2017)). These likelihood models have at least three edges over the transformation-based models. First, the likelihood-based models will preserve the useful and important matrix structural information, which makes them more interpretable compared with transformation-based models. Second, the number of estimated parameters in the transformation-based models has order $O(n^4)$, while the one in the likelihood-based models has order $O(n^2)$, where $n$ is the dimension of the RCOV matrix. When $n$ is large, the likelihood-based models can bring more convenience and a less daunting task in computation. Third, the likelihood-based models make use of the likelihood function of the RCOV matrices, and hence their statistical inference methods could be easily provided.

This paper contributes to the literature from three aspects. First, we propose a new Conditional BEKK matrix-F (CBF) model to study the time-varying RCOV matrices. Our CBF model has matrix-F distributed innovations with two degrees of freedom parameters $\nu_1$ and $\nu_2$. When $\nu_2\to\infty$, our CBF model reduces to the CAW model (Golosnoy et al., 2012), which has Wishart distributed innovations. Hence, the degrees of freedom $\nu_2$ is designed to capture the heavy-tailedness of the RCOV. Since the RCOV is also well documented to have long memory phenomenon, we further introduce a special CBF model which has a similar conditional heterogeneous autoregressive (HAR) structure as in Corsi (2009). This special model is coined the CBF-HAR model. Although the CBF-HAR model is not formally a long memory model, it gives rise to persistence in the RCOV time series. Two real examples demonstrate that our CBF model (especially the CBF-HAR model) can have a significantly better forecasting performance than the corresponding CAW model, and hence a simple incorporation of $\nu_2$ to capture the heavy-tailed RCOV is necessary from a practical viewpoint.

Second, we provide a systematically statistical inference procedure for the CBF model. Specifically, we explore its stationarity conditions, establish the strong consistency and asymptotic normality of its maximum likelihood estimator (MLE), and investigate some new inner-product-based tests for model diagnostic checking. Moreover, the performance of our entire methodology is assessed by simulation studies. Compared to the existing BEKK-type multivariate time series models, our proofs of the entire inference procedure are much involved, since the CBF model is tailored for matrix time series. Particularly, our inner-product-based tests seem to be the first diagnostic checking tool for matrix time series models, and the related idea can be easily extended to other models.

Third, we construct two reduced CBF models --- the variance targeted (VT) CBF (VT-CBF) model and the factor CBF (F-CBF) model, to handle moderately large and high dimensional RCOV matrix respectively. For both reduced models, the asymptotic theory of the estimated parameters is derived. The dimension of the RCOV matrix is allowed to be a moderate but fixed number in the VT-CBF model, while it is allowed to grow with the sample size $T$ and the intra-day sample size in the F-CBF model. Therefore, this makes the prediction of large dimensional RCOV matrices feasible in many cases. The importance of both reduced models is illustrated by two real applications.

The remainder of this paper is organized as follows. Section 2 introduces the CBF model and studies its probabilistic properties. Section 3 investigates the asymptotics of the MLE. Section 4 presents inner-product-based tests to check the model adequacy. Two reduced CBF models and their related asymptotic theories are provided in Section 5. Some simulation studies are carried out in Section 6. Applications are given in Section 7. Section 8 concludes this paper. Proofs of all theorems are relegated to the supplementary material.

Some notations are used throughout the paper. $I_{n}$ is the identity matrix of order $n$, and $\otimes$ represents the Kronecker product. For an $n\times n$ matrix $A$, $tr(A)$ is its trace, $A'$ is its transpose, $\left|A\right|$ is its determinant, $\rho(A)$ is its biggest eigenvalue, $\left\Vert A\right\Vert =\sqrt{tr(A'A)}$ is its Euclidean (or Frobenius) norm, $\left\Vert A\right\Vert _{spec}=\sqrt{\rho(A'A)}$ is its spectral norm, $vec(A)$ is a vector obtained by stacking all the columns of $A$, $vech(A)$ is a vector obtained by stacking all columns of the lower triagular part of $A$, and $A^{\otimes2}=A\otimes A$.

Model and Properties

Model specification

Let $Y_{t}^*$ be the integrated volatility matrix of $n$ asset returns $X_{t}$ at time $t=1,...,T$. After the seminal work of Barndorff-Nielsen and Shephard (2002, 2004) and Andersen et al. (2003), the $n\times n$ positive definite realized covariance (RCOV) matrix $Y_{t}$ calculated from the high-frequency return data of $X_{t}$ has been widely applied to estimate $Y_{t}^*$ in the literature; see, e.g., Barndorff-Nielsen et al. (2011), Lunde et al. (2016), A\"{i}t-Sahalia and Xiu (2017), Kim et al. (2018) and references therein. Moreover, $Y_{t}$ is often viewed as a precise estimate for the conditional variances and covariances of these $n$ low-frequency asset returns $X_{t}$, and hence how to predict $Y_{t}$ by some dynamic models is important in practice. Motivated by this, a new dynamic model for $Y_t$ is proposed in the current paper.

Let $\mathcal{G}_{t}=\sigma(Y_{s}; s\leq t)$ be a filtration up to time $t$. We assume that

where $\{\Delta_{t}\}_{t=1}^{T}$ is a sequence of independent and identically distributed (i.i.d.) $n\times n$ positive definite random innovation matrices with $E(\Delta_{t}|\mathcal{G}_{t-1})=I_{n}$, each $\Delta_{t}$ follows the matrix-F distribution $F(\nu,\frac{\nu_{2}-n-1}{\nu_{1}}I_{n})$, and the density of $F(\nu,\Sigma)$ is

where $\nu=\left(\nu_{1},\nu_{2}\right)'$ with degrees of freedom $\nu_{1}>n+1$ and $\nu_{2}>n+1$, $\Sigma$ is an $n\times n$ positive definite matrix, and

moreover, $\Sigma_{t}^{1/2}\in\mathcal{G}_{t-1}$ is the square root of the $n\times n$ positive definite matrix $\Sigma_{t}$, which has a BEKK-type dynamic structure (see Engle and Kroner, 1995):

where $\Omega$, $A_{ki}$, $B_{kj}$ are all $n\times n$ real matrices, the integers $P,Q,K$ are known as the orders of the model, and $\Omega$ as well as the initial states $\Sigma_{0},\Sigma_{-1},...,\Sigma_{-Q+1}$ are all positive definite. Under model ((ref)),

with $E(Y_{t}|\mathcal{G}_{t-1})=\Sigma_{t}$, that is, the conditional distribution of $Y_{t}$ is matrix-F with a BEKK-type mean structure. In this sense, we call model ((ref)) the Conditional BEKK matrix-F (CBF) model.

The CBF model is related to the CAW model in Golosnoy et al. (2012), in which $\Delta_{t}$ follows the Wishart distribution. To see it clearly, we follow Konno (1991) and Leung and Lo (1996) to re-write $Y_{t}$ in model ((ref)) as

where $L_{t}\sim\mbox{Wishart}(\nu_{1},I_{n})$ and $R_{t}\sim\mbox{Wishart}(\nu_{2},I_{n})$ are independent. As $\lim \limits_{\nu_{2}\to\infty}\nu_{2}^{-1}R_{t}$ $=I_{n}$ in probability, the identity ((ref)) implies that when $\nu_{2}\rightarrow\infty$, $Y_{t}|\mathcal{G}_{t-1}\sim \mbox{Wishart}(\nu_{1}, $ $\nu_{1}^{-1}\Sigma_{t})$, which is exactly the CAW model. Therefore, compared to the CAW model, the degrees of freedom $\nu_{2}$ in the CBF model accommodates the heavy-tailed RCOV, meaning that each $Y_{t,ij}$ from $Y_t$ satisfying ((ref)) could have a heavier tail than that from $Y_t$ satisfying $Y_{t}|\mathcal{G}_{t-1}\sim \mbox{Wishart}(\nu_{1}, \nu_{1}^{-1}\Sigma_{t})$ (see, e.g., Opschoor et al. (2018) for more discussions and examples). Clearly, the identity ((ref)) also guarantees $Y_{t}$ to be symmetric and positive definite, and it can be used to generate $Y_{t}$ by using Wishart random variables.

Besides the heavy-tailedness, long memory is another well documented feature for the RCOV, and it has been taken into account by many RCOV models, including the heterogeneous autoregressive (HAR) model in Corsi (2009) as a benchmark. Although the HAR model does not formally belong to the class of long memory models, it is able to reproduce the persistence of RCOV observed in many empirical data. Inspired by the HAR model, we consider a special CBF model, which has the following specification for $\Sigma_t$:

where $Y_{t-1,d}=Y_{t-1}$, $Y_{t-1,w}=(1/5)\sum_{i=1}^{5}Y_{t-i}$, and $Y_{t-1,m}=(1/22)\sum_{i=1}^{22}Y_{t-i}$ are the daily, weekly, and monthly averages of RCOV matrices, respectively. In this case, we label model ((ref)) as the CBF-HAR model, since we put “HAR dynamics” on $\Sigma_t$. Clearly, the CBF-HAR model is simply a constrained CBF model with $P=22$, $K=3$ and $Q=0$. Figure (ref) plots the sample autocorrelation functions (ACFs) up to lag 100 of one simulated data from the CBF-HAR model with $\nu=(20,10)$ and

From this figure, we can find that all entries of $Y_{t}$ exhibit long memory phenomenon as expected.

Note that when $K=1$, sufficient identifiability conditions of model ((ref)) are that the main diagonal elements of $\Omega$ and the first diagonal element of each $A_{1i}$, $B_{1j}$ are positive; when $K>1$, some sufficient identifiability conditions of model ((ref)) can be found in Engle and Kroner (1995). For simplicity, we assume subsequently that model ((ref)) is identifiable.

Of course, the BEKK specification in model ((ref)) is not the only way to describe the dynamics of $\Sigma_{t}$. The multivariate ARCH-type models such as the VEC model in Bollerslev et al. (1988), the component model in Engle and Lee (1999), the dynamic conditional correlation model in Engle (2002) and many others can also be adopted to model $\Sigma_{t}$. Using these models together with the matrix-F distribution to fit and predict the RCOV matrices could be a promising direction for future study.

Stationarity

Stationarity is an important issue for most RCOV models, but so far it has been rarely studied. Denote $M=max(P,Q)$. For $i=1,2,\cdots,M$, let $$A_{i}^{*}=\sum_{k=1}^{K}A_{ik}^{\otimes2}\,\,\,\mbox{ and }\,\,\,B_{i}^{*}=\sum_{k=1}^{K}B_{ik}^{\otimes2},$$ where $A_{ik}=0$ for $i>P$ and $B_{ik}=0$ for $i>Q$. A sufficient condition for the stationarity of the CBF model is given below, and it works for other general distributions of $\Delta_{t}$.

As a special case, the results in Theorem (ref) hold for the CAW model, in which $\Delta_{t}$ follows the Wishart distribution. Under conditions (H1) and (H2), condition (H3) is necessary and sufficient for the strict stationarity of $Y_{t}$ with $E\|Y_{t}\|<\infty$. However, the necessary and sufficient condition for the higher moments of $Y_t$ is still unclear at this stage. Let $K_{n^2}$ be the $n^{2}\times n^{2}$ permutation matrix such that $K_{n^2}vec(A)=vec(A')$ for any $n\times n$ matrix $A$. If $E\|Y_{t}\|^2<\infty$, it is not hard to see

(i) $\overline{y}:=E\left(vec(Y_{t})\right)=\Big[I_{n^{2}}-\sum_{i=1}^{M}\left(A_{i}^{*}+B_{i}^{*}\right)\Big]^{-1}vec(\Omega)$;

(ii) $vec\left[E\left(vec(Y_{t})vec(Y_{t})'\right)\right]=\left(\Pi+I_{n^{4}}\right) \left(I_{n^{4}}-\sum_{i=1}^{\infty}\Phi_{i}^{\otimes2}\Pi\right)^{-1}vec(\overline{y})\otimes vec(\overline{y})$,

where $\Pi=\left[s_{1}(\nu)-1\right]I_{n^{4}}+\left[s_{2}(\nu)I_{n^{2}}\otimes(I_{n^{2}}+K_{n^2})\right]\left[I_{n}\otimes K_{n^2}\otimes I_{n}\right]$ with $$ s_{1}(\nu)=\frac{(\nu_{2}-n-1)[\nu_{1}(\nu_{2}-n-2)+2]}{\nu_{1}(\nu_{2}-n)(\nu_{2}-n-3)},\,\,\, s_{2}(\nu)=\frac{(\nu_{2}-n-1)(\nu_{1}+\nu_{2}-n-1)}{\nu_{1}(\nu_{2}-n)(\nu_{2}-n-3)}, $$ and $\Phi_{0}=I_{n^{2}}$, $\Phi_{i}=-B_{i}^{*}+\sum_{j=1}^{i}\big(A_{j}^{*}+B_{j}^{*}\big)\Phi_{i-j}$ for $i>0$. Result (ii) above clearly indicates that the parameters $\nu_1$ and $\nu_2$ have impact on the second moment of $Y_t$ in a non-linear way. Although a closed form of third moment of $Y_t$ is absent, similar impact from $\nu_1$ and $\nu_2$ is expected for the third moment of $Y_t$ and hence the asymptotic distribution of the proposed estimator (see Theorem (ref) below).

Maximum Likelihood Estimation

Let $\theta=(\gamma',\nu')'\in\Theta$ be the unknown parameter of model ((ref)) with the true value $\theta_{0}=(\gamma_{0}',\nu_{0}')'$, where $\Theta=\Theta_{\gamma}\times\Theta_{\nu}$ is the parametric space with $\Theta_{\gamma}\subset\mathbb{R}^{\tau_{1}}$ and $\Theta_{\nu}\subset\mathbb{R}^{2}$, $\gamma=(w', u')'$, $w=vech(\Omega)$, $u=(vec(A_{11})',...,vec(A_{KP})', vec(B_{11})',...,vec(B_{KQ})')$, and $\tau_{1}=\frac{1}{2}n+[(P+Q)K+\frac{1}{2}]n^2$. Below, we assume that $\Theta_{\gamma}$ and $\Theta_{\nu}$ are compact and $\theta_{0}$ is an interior point of $\Theta$.

Given the observations $\{Y_{t}\}_{t=1}^{T}$ and the initial values $\{Y_{t}\}_{t\leq 0}$, the negative log-likelihood function based on ((ref)) is

where

with $C(\nu)=-\log \Lambda(\nu)$ and $\Sigma_{t}(\gamma)$ calculated recursively by

Clearly, $\Sigma_{t}(\gamma_{0})=\Sigma_{t}$.

As the initial values $\{Y_{t}\}_{t\leq 0}$ are not observable, we shall modify $L(\theta)$ as

where $\widehat{l}_{t}(\theta)$ is defined in the same way as $l_{t}(\theta)$ with $\Sigma_{t}(\gamma)$ being replaced by $\widehat{\Sigma}_{t}(\gamma)$, and $\widehat{\Sigma}_{t}(\gamma)$ is calculated in the same way as $\Sigma_{t}(\gamma)$ based on a sequence of given constant matrices $h:=\{Y_{0},\cdots,Y_{-M+1},\Sigma_{0},...,\Sigma_{-M+1}\}$. The minimizer, $\widehat{\theta}=(\widehat{\gamma}',\widehat{\nu}')'$, of $\widehat{L}(\theta)$ on $\Theta$ is called the maximum likelihood estimator (MLE) of $\theta_{0}$. That is,

To study the asymptotic properties of $\widehat{\theta}$, we need two assumptions below.

Assumption (ref) is standard, and Assumption (ref) which is in line with Comte and Lieberman (2003) and Hafner and Preminger (2009) is the identification condition. The following two theorems give the consistency and asymptotic normality of $\widehat{\theta}$, respectively.

Based on the observations $\{Y_{t}\}_{t=1}^{T}$ and a sequence of given constant matrices $h$, we can use the analytic expression of $\partial^2 l_{t}(\theta) /( \partial\theta\partial\theta')$ (see Appendix D in the supplementary material) to estimate $\mathcal{O}$ by its sample counterpart. As the univariate ARCH-type models, the coefficients in the main diagonal line of $\Omega$ are positive to ensure the positive definite of $\Sigma_t$. Hence, the classical $t$ or Wald test, which is constructed by the estimate of $\mathcal{O}$, can not be used to detect whether their values are zeros; see Li et al. (2018) for more discussions on this context.

Model Diagnostic Checking

Diagnostic tests are crucial for model checking in multivariate time series analysis; see, e.g., Li and McLeod (1981), Ling and Li (1997), Tse (2002) and many others. However, no attempt has been made for the stationary matrix time series. In this section, we propose some new inner-product-based tests to check the adequacy of model ((ref)).

Let $\mathfrak{Z}_{t}(\gamma) = vec\big(\Sigma_{t}^{-1/2}(\gamma)Y_{t}\Sigma_{t}^{-1/2}(\gamma)-I_{n}\big)$ be the vectorized residual for a given $\gamma$, and $ \mathbf{b}_{t,j}(\gamma)=\mathfrak{Z}_{t}'(\gamma)\mathfrak{Z}_{t-j}(\gamma)$ be the inner product of two vectorized residuals at lag $j$. Then, we stack $\mathbf{b}_{t,j}(\gamma)$ up to lag $l$ to construct $\mathcal{V}_{l}(\gamma)$, where $$ \mathcal{V}_{l}(\gamma)=\frac{1}{T}\sum_{t=l+1}^{T}\big( \mathbf{b}_{t,1}\left(\gamma\right), \mathbf{b}_{t,2}\left(\gamma\right), \cdots, \mathbf{b}_{t,l}\left(\gamma\right)\big)', $$ and $l\geq1$ is a given integer. Our testing idea is motivated by the fact that if model ((ref)) is adequate, $\mathfrak{Z}_{t}(\gamma_{0})$ is a sequence of i.i.d. random vectors with mean zero, and hence the value of $\mathcal{V}_{l}(\widehat{\gamma})$ is expected to be close to zero. To implement our test, we need study the asymptotic property of $\mathcal{V}_{l}(\widehat{\gamma})$ in the following theorem.

Based on Theorem (ref), we construct the inner-product-based test statistic

to detect the adequacy of model ((ref)), where $\widehat{\mathbf{V}}$ is the sample counterpart of $\mathbf{V}$. If $\Pi(l)$ is larger than the upper-tailed critical value of $\chi^{2}(l)$, the fitted model ((ref)) is not adequate at a given significance level. Otherwise, it could be deemed as adequate.

Note that if we consider a test based on $\{\mathfrak{Z}_{t}(\widehat{\gamma})\}$ directly, the resulting limiting distribution shall still be chi-squared, but its degrees of freedom increases fast with the dimension $n$. To avoid this dilemma, we use the inner product of the residuals to construct our test $\Pi(l)$, whose limiting distribution is independent of $n$. This new idea is different from the portmanteau test in Ling and Li (1997) in which the test statistic is constructed based on the auto-correlations of the transformed scale residuals, while our test $\Pi(l)$ is based on the auto-covariances of the original vectorized residuals. Clearly, our idea can be easily extended to the framework in Ling and Li (1997). Meanwhile, our inner-product-based test $\Pi(l)$ takes the auto-covariances of all entries of $\mathfrak{Z}_{t}(\widehat{\gamma})$ into account, while the idea of regression-based test in Tse (2002) only considers one entry of $\mathfrak{Z}_{t}(\widehat{\gamma})$ at a time. In view of this, we prefer to use the proposed inner-product idea for testing purpose.

The Reduced CBF Models

As the number of parameters in the CBF model is $O(n^2)$, the estimation of the CBF model could be very computationally demanding when $n$ is large. This section introduces two reduced CBF models, which are feasible in fitting RCOV matrices with a large $n$.

The VT-CBF model

This subsection proposes a reduced CBF model by using the variance target (VT) technique in Engle and Mezrich (1996). The idea of VT is to re-parameterize the drift matrix $\Omega$ by using the theoretical mean of $Y_{t}$, so that the estimation of $\Omega$ is excluded in the implementation of the maximum likelihood estimation. Other related studies on the VT time series models can be found in Francq et al. (2011) and Pedersen and Rahbek (2014).

To define our reduced model, we assume that $Y_{t}$ is strictly stationary with a finite mean $S=E(Y_{t})$. By taking expectation on both sides of ((ref)), we have

due to the fact that $S=E(Y_{t})=E(\Sigma_{t})$. With the help of ((ref)), model ((ref)) becomes

where all notations are inherited from model ((ref)), except that

We call model ((ref)) the VT-CBF model. Clearly, this reduced model shares the same probabilistic properties as the full CBF model. Although the VT-CBF model has the same amount of parameters as the full CBF model, its two-step estimator given below is computationally easier than the MLE for the full CBF model.

To present this two-step estimator, we let $\theta_{v}=(\delta',\nu')'\in \Theta_{v}$ be the unknown parameters of model ((ref)) and its true value be $\theta_{v0}=(\delta_{0}',\nu_{0}')'$, where $\Theta_{v}=\Theta_{\delta}\times\Theta_{\nu}$ is the parametric space with $\Theta_{\delta}=\Theta_{s}\times\Theta_{u}\subset\mathbb{R}^{\tau_{2}}$, $\tau_{2}=[(P+Q)K+1]n^2$ and $\Theta_{\nu}\subset\mathbb{R}^{2}$. Let $\delta=(s',u')'$ with $s= vec(S)$, $\Theta_s\in\mathbb{R}^{n^2}$ and $\Theta_u\in\mathbb{R}^{[(P+Q)Kn^2]}$. As before, we assume that $\Theta_{\delta}$ and $\Theta_{\nu}$ are compact and $\theta_{v0}$ is an interior point of $\Theta_{v}$.

In the first step, we estimate $s$ by $\widehat{s}_{v}$, where $\widehat{s}_{v}=vec\left(\overline{Y_t}\right):=vec\big(\frac{1}{T}\sum_{t=1}^{T}Y_{t}\big)$. In the second step, we estimate the remaining parameters $\zeta=(u',\nu')'$ by the constrained MLE based on the following modified log-likelihood function:

where

and $\widehat{\Sigma}_{vt}(\delta)$ is calculated recursively by

based on a sequence of given constant matrices $h$. Clearly, $\widehat{L}_{v}(\theta_{v})$ is analogous to $\widehat{L}(\theta)$ in ((ref)), and it is the modification of the following log-likelihood function:

where $l_{vt}(\theta_{v})$ is defined in the same way as $\widehat{l}_{vt}(\theta_{v})$ with $\widehat{\Sigma}_{vt}(\delta)$ being replaced by $\Sigma_{vt}(\delta)$, and $\Sigma_{vt}(\delta)$ is calculated recursively by

based on the observations $\{Y_{t}\}_{t=1}^{T}$ and the initial values $\{Y_{t}\}_{t\leq 0}$. The minimizer, $\widehat{\zeta}_{v}=(\widehat{u}_{v}',\widehat{\nu}_{v}')'$, of $\widehat{L}_{v}(\widehat{s}_{v},\zeta)$ on $\Theta_{u}\times\Theta_{\nu}$ is the constrained MLE of $(u_{0}', \nu_{0}')'$. That is,

Now, we call $\widehat{\theta}_{v}=(\widehat{s}_{v}',\widehat{\zeta}_{v}')'$ the two-step estimator of $\theta_{v}$ in model ((ref)). Let $\Psi(u)=\big(I_{n^{2}}-\sum_{i=1}^{M}A_{i}^{*}-\sum_{i=1}^{M}B_{i}^{*}\big)^{-1}\big(I_{n^{2}}-\sum_{i=1}^{M}B_{i}^{*}\big)$ and $w_{t}(\theta_{v})=\Big(

\Big)$. The following two theorems give the consistency and asymptotic normality of $\widehat{\theta}_{v}$, respectively.

As before, we can use the sample counterpart of the analytic expressions of $\partial l_{vt}(\theta_{v})/\partial\theta_{v}$ and $\partial^2 l_{vt}(\theta_{v})/\partial\theta_{v}\partial\theta_{v}'$ to estimate $\mathcal{O}_{v}$. Although the VT-CBF model can be estimated by the aforementioned two-step estimation procedure, it still has to handle a large number of estimated parameters with order $O(n^{2})$ caused by the parameter matrices $A_{ki}$ and $B_{kj}$. To make a more parsimonious VT-CBF model, we can further impose some restrictions on $A_{ki}$ and $B_{kj}$. McCurdy and Stengos (1992) and Engle and Kroner (1995) have suggested to use diagonal volatility models, which not only avoid over-parameterization, but also reflect the fact that the variances and the covariances rely more on its own past than the history of other variances or covariances. Motivated by this, we can assume that all $A_{ki}$ and $B_{kj}$ have a diagonal structure, leading to a diagonal VT-CBF model. Clearly, the number of estimated parameters in the diagonal VT-CBF model has order $O(n)$, which is feasible to be handled for a moderate large but fixed $n$.

Next, similar to $\Pi(l)$ in ((ref)), we can construct the inner-product-based test statistics to check the adequacy of model ((ref)) based on the two-step estimator $\widehat{\theta}_{v}$. Let $\delta_{0}=(s_{0}',u_{0}')'$, $\widehat{\delta}_{v}=(\widehat{s}_{v}',\widehat{u}_{v}')'$, $\mathfrak{Z}_{vt}(\delta) = vec(\Sigma_{vt}^{-1/2}(\delta)Y_{t}\Sigma_{vt}^{-1/2}(\delta)-I_{n})$ be the residual vector for a given $\delta$, $ \mathbf{b}_{vt,j}(\delta)=\mathfrak{Z}_{vt}'(\delta)\mathfrak{Z}_{vt-j}(\delta)$ be the inner product of the residuals at lag $j$, and $$ \mathcal{V}_{vl}(\delta)=\frac{1}{T}\sum_{t=l+1}^{T}\big( \mathbf{b}_{vt,1}\left(\delta\right), \mathbf{b}_{vt,2}\left(\delta\right), \cdots, \mathbf{b}_{vt,l}\left(\delta\right)\big)'. $$ The asymptotic property of $\mathcal{V}_{vl}(\widehat{\delta}_{v})$ is given in the following theorem.

By the preceding theorem, we can adopt the test statistic

to detect the adequacy of model ((ref)), where $\widehat{\mathbf{V}}_{v}$ is the sample counterpart of $\mathbf{V}_{v}$. If $\Pi_{v}(l)$ is larger than the upper-tailed critical value of $\chi^{2}(l)$ at a given significance level, the fitted model ((ref)) is not adequate. Otherwise, it is adequate.

The factor CBF model

In modern data analysis, the dimension $n$ could be growing with the sample size $T$ in many cases, and this makes the CBF (or VT-CBF) models computationally infeasible. Also, the dimension $n$ may be proportional to $m$ (the average intra-day sample size across all assets and all days), and then the methods to calculate $Y_{t}$ used for the fixed $n$ deliver an inconsistent estimator of $Y_{t}^{*}$; see, e.g., Wang and Zou (2010) and Tao et al. (2011) for surveys. To overcome this difficulty, we use the thresholding average realized volatility matrix (TARVM) estimator in Tao et al. (2011) to calculate $Y_{t}$. The TARVM is built based on the ARVM (Wang and Zou, 2010), which is constructed by taking the average of the constructed realized volatility matrices according to different predetermined sampling frequencies. The TARVM further thresholds the elements in each estimated RCOV matrix from the ARVM method, so that certain sparsity structure is retained and the resulting estimator is consistent for large $n$, which can be growing with (or even larger than) $T$. For more recent works in this direction, we refer to A\"{i}t-Sahalia and Xiu (2017), Kim et al. (2018), and the references therein.

Since the dimension of $Y_{t}$ could be very large, it seems hard to study the dynamics of $Y_{t}$ without imposing some specific structure. Here, we adopt the factor model proposed by Tao et al. (2011) by assuming that

where $Y_{ft}^{*}$ is an $r\times r$ positive definite factor covariance matrix with $r$ being a fixed integer (much smaller than $n$), $Y_{0}^{*}$ is an $n\times n$ positive definite constant matrix, and $F$ is an $n\times r$ factor loading matrix normalized by the constraint $F'F=I_{r}$. In model ((ref)), the dynamic structure of $Y_{t}^{*}$ is driven by that of a lower-dimensional latent process $Y_{ft}^{*}$, while $Y_{0}^{*}$ represents the static part of $Y_{t}^{*}$.

In ((ref)), we shall highlight that only the column space of $F$ can be identified, and $F$ is not identified even if $F'F=I_r$ is imposed. This is because $Y_{t}^{*}$ is unchanged when $F$ and $Y_{ft}^*$ are replaced by $F_{\dag}=FR$ and $Y_{ft,\dag}^{*}=R^{-1}{Y}_{ft}^*R^{-1'}$, respectively, where $R$ is any $r\times r$ matrix satisfying $R'R=I_r$.

Define

and

Then, we estimate $Y_{ft}^{*}$, $Y_{0}^{*}$ and $F$ by

respectively, where $\widehat{f}_{1},\cdots,\widehat{f}_{r}$ are the eigenvectors of $\overline{S}$ corresponding to its $r$ largest eigenvalues. As suggested by Lam and Yao (2012) and Ahn and Horenstein (2013), we may select $r$ such that the $r$ largest ratios of adjacent eigenvalues are significantly larger.

In order to study the asymptotics of the proposed estimators, we introduce the following technical assumptions.

Assumptions (ref)-(ref) are sufficient to prove the consistency of $\widehat{Y}_{ft}$. For TARVM, we can take $A(n,m,T)=\pi(n)[e_{m}(n^2T)^{1/\beta}]^{1-\delta_*}\log T$ and $B(T)=\log T$ with $e_m=m^{-1/6}$ so that $A(n,m,T)B^{5}(T)=o(1)$ for large $\beta$; see Tao et al. (2011). Note that Assumptions (ref)-(ref) do not rule out the case that $n$ is larger than $T$, as long as $n^2T$ grows more slowly than $m^{\beta/6}$. For other estimators, the rate $A(n,m,T)$ may be improved; see Tao et al. (2013) for more discussions.

The above theorem indicates that $\widehat{Y}_{ft}$ is a consistent estimator of $Y_{ft}$ rather than $Y_{ft}^{*}$. Next, we assume that $Y_{ft}$ satisfies the CBF model, that is,

with $E(Y_{ft}|\mathcal{G}_{t-1})=\Sigma_{ft}$, where $\Sigma_{ft}$ is defined in the same way as $\Sigma_{t}$ in ((ref)) with $Y_{t}$ replaced by $Y_{ft}$, and the remaining notations and set-ups inherent from model ((ref)). We call models ((ref)) and ((ref)) the factor CBF (F-CBF) model. Particularly, if $\Sigma_{ft}$ has the HAR dynamical structure as in ((ref)), the resulting model is called the factor CBF-HAR (F-CBF-HAR) model. Based on this model, we have $Y_{t}^{*}=F(Y_{ft}-F'Y_{0}^{*}F)F'+Y_{0}^{*}$. Since $Y_{t}\approx Y_{t}^{*}$, it implies that we can study the large dimensional matrix $Y_{t}$ by using an $r\times r$ low-dimensional matrix $Y_{ft}$.

As $Y_{ft}$ is not observable, we should estimate model ((ref)) based on $\widehat{Y}_{ft}$, and hence we consider a feasible MLE of $\theta_{0}$ in model ((ref)) given by

where $\widehat{L}_{f}(\theta)$ is defined in the same way as $\widehat{L}(\theta)$ in ((ref)) with $Y_{t}$ and $\widehat{\Sigma}_{t}(\gamma)$ replaced by $\widehat{Y}_{ft}$ and $\widehat{\Sigma}_{ft}(\gamma)$, respectively. The following theorem shows that $\widehat{\theta}_{1f}$ is consistent with the ideal MLE $\widehat{\theta}_{2f}$ based on $Y_{ft}$, where

and $L_{f}(\theta)$ is defined in the same way as $L(\theta)$ in ((ref)) with $Y_{t}$ and $\Sigma_{t}(\gamma)$ replaced by $Y_{ft}$ and $\Sigma_{ft}(\gamma)$, respectively.

Since the dimension of $Y_{ft}$ is $r$ (much smaller than $n$), the calculation of $\widehat{\theta}_{1f}$ is computationally feasible. In order to further reduce the number of parameters in model ((ref)), we can also assume that $Y_{ft}$ follows a VT-CBF model. This leads to the F-VT-CBF model, which includes the F-VT-CBF-HAR model as a special case. For this F-VT-CBF model, we consider its feasible two-step estimator $\widehat{\theta}_{1fv}=(\widehat{s}_{1fv}',\widehat{\zeta}_{1fv}')'$, where

and $\widehat{L}_{fv}(\theta_{v})$ is defined in the same way as $\widehat{L}_{v}(\theta_{v})$ in ((ref)) with $Y_{t}$ and $\widehat{\Sigma}_{vt}(\delta)$ replaced by $\widehat{Y}_{ft}$ and $\widehat{\Sigma}_{fvt}(\delta)$, respectively. Similar to Theorem (ref), $\widehat{\theta}_{1fv}$ is consistent with the ideal two-step estimator $\widehat{\theta}_{2fv}=(\widehat{s}_{2fv}',\widehat{\zeta}_{2fv}')'$ based on $Y_{ft}$, where

and $L_{fv}(\theta_{v})$ is defined in the same way as $L(\theta_{v})$ in ((ref)) with $Y_{t}$ and $\Sigma_{t}(\delta)$ replaced by $Y_{ft}$ and $\Sigma_{fvt}(\delta)$, respectively.

Particularly, if $Y_{ft}$ follows a diagonal VT-CBF model, the number of estimated parameters in model ((ref)) is $O(r)$, which is easy to calculate in practice. In view of model ((ref)) and the fact that $F'F=I_{r}$, we can predict $Y_{t}$ by either $\widehat{F}\widehat{\Sigma}_{ft}(\widehat{\gamma}_{1f})\widehat{F}'+\widehat{Y}_{0}^{*}$ based on $\widehat{\theta}_{1f}$ or $\widehat{F}\widehat{\Sigma}_{fvt}(\widehat{\delta}_{1fv})\widehat{F}'+\widehat{Y}_{0}^{*}$ based on $\widehat{\theta}_{1fv}$, where $\widehat{\delta}_{1fv}=(\widehat{s}_{1fv}',\widehat{u}_{1fv}')'$.

Simulation

In this section, we first assess the performance of the MLE $\widehat{\theta}$ and the two-step estimator $\widehat{\theta}_{v}$ in the finite sample. We generate 1000 replications of sample size $T=1000$ and $2000$ from the following model:

where

$\{\Delta_{t}\}$ is a sequence of independent $F\big(\nu_0,\frac{\nu_{20}-n-1}{\nu_{10}}I_{n}\big)$ distributed random matrices with $n=3$, and $\nu_0=(10, 8), (15, 10)$ or $(20, 10)$. For each repetition, we calculate $\widehat{\theta}$, $\widehat{\theta}_{v}$, and their related asymptotic standard deviations. For $\widehat{\theta}_{v}$, we report the results related to $\Omega$ instead of $S$, and hence the asymptotic standard deviation of the estimated parameters in $\Omega$ is absent in this case.

Table (ref) reports the sample bias, the sample standard deviation (SD) and the average asymptotic standard deviation (AD) of $\widehat{\theta}$ and $\widehat{\theta}_{v}$. From this table, we can see that the biases of both estimators are small comparing to the magnitude of the parameters, and they become smaller as the sample size $T$ increases. This assures the accuracy of both estimators. Furthermore, we find that the SDs are generally close to the ADs for both estimators, and all of the SDs and ADs become smaller as $T$ increases from 1000 to 2000. In terms of ADs or SDs, $\widehat{\theta}$ is generally more efficient than $\widehat{\theta}_{v}$, although this efficiency advantage is weak for many parameters. However, the estimation time for $\widehat{\theta}_{v}$ is almost 70% of the time for $\widehat{\theta}$, and this computation advantage can be more significant when $n$ increases.

Next, we examine the performance of the inner-product-based tests $\Pi(l)$ and $\Pi_{v}(l)$ in the finite sample. We generate 1000 replications of sample size $T=1000$ and $2000$ from the following model:

where the values of $\Omega_0$, $A_{10}$ and $B_{10}$ are chosen as in ((ref)), $A_{20}=\mbox{diag}\{\lambda,\lambda,\lambda\}$ is a diagonal matrix with $\lambda=0, 0.05, 0.1, 0.15, 0.2$, and $\{\Delta_{t}\}$ is a sequence of independent $F\big(\nu_0,\frac{\nu_{20}-n-1}{\nu_{10}}I_{n}\big)$ distributed random matrices with $n=3$ and $\nu_0=(10, 8)$. We fit each replication by the CBF model with $(K, P, Q)=(1, 1, 1)$, and use $\Pi(l)$ and $\Pi_{v}(l)$ to check whether the fitted model is adequate. Here, we set the significance level $\alpha=0.05$ and $l=2,3,4,5,6$. The empirical sizes and powers of both tests are reported in Table (ref), and their sizes correspond to the results for the case of $\lambda=0$. From Table (ref), we can find that both $\Pi(l)$ and $\Pi_{v}(l)$ always have accurate sizes, although they are slightly oversized for small $T$. For the power of both tests, it is generally as expected. First, all the power becomes larger as $T$ increases. Second, both tests become more powerful as $\lambda$ becomes larger. Third, the power of $\Pi(l)$ and $\Pi_{v}(l)$ is comparable, but the former need more computational time. Note that when $\nu_0=(15, 10)$ and $(20, 10)$, the testing results are similar to these for $\nu_0=(10, 8)$, and hence they are not reported for saving space.

Overall, both estimators $\widehat{\theta}$ and $\widehat{\theta}_{v}$ and both tests $\Pi(l)$ and $\Pi_{v}(l)$ have a good performance especially when the sample size $T$ gets larger. When the dimension of $Y_{t}$ is small, our simulation results show that $\widehat{\theta}_{v}$ is only slightly less efficient than $\widehat{\theta}$, and $\Pi_{v}(l)$ is generally as powerful as $\Pi(l)$. When the dimension of $Y_{t}$ is large, $\widehat{\theta}_{v}$ and $\Pi_{v}(l)$ can enjoy a faster computation speed than $\widehat{\theta}$ and $\Pi(l)$, respectively. Based on these grounds, we would recommend using $\widehat{\theta}_{v}$ and $\Pi_{v}(l)$ in practice.

Applications

In this section, we consider two applications on the U.S. stock market. Application 1 studies the low dimensional RCOV matrix series calculated by composite realized kernels (CRK) in Lunde, Shephard and Sheppard (2016). Application 2 studies the high dimensional RCOV series calculated by TARVM estimator in Tao et al. (2011).

Application 1

In this application, we revisit the RCOV matrix data of Hewlett-Packard Development Company, L.P. (HPQ), International Business Machines Corporation (IBM) and Microsoft Corporation (MSFT) in Lunde, Shephard and Sheppard (2016). This data set, denoted by $\{Y_t\}_{t=1}^{1474}$, ranges from January 2006 to December 2011 with 1474 observations in total. Here, two flash crashes are flagged in 6 May, 2010 and 9 August, 2011 and replaced by an average of the nearest five preceding and following matrices.

Figure (ref) plots the diagonal and off-diagonal components of $\{Y_{t}\}_{t=1}^{1474}$, exhibiting that $Y_t$ has a clear clustering feature. Meanwhile, Figure (ref) plots their sample autocorrelation functions (ACFs), which show the significant temporal dependence of $Y_t$. Based on these facts, we first fit $\{Y_{t}\}_{t=1}^{1474}$ by a diagonal VT-CBF model with $(P,Q,K)=(3,1,1)$, where the order $K$ is taken as one for ease of model identification, and the orders $P$ and $Q$ are selected by the Bayesian information criterion (BIC). Specifically, this diagonal VT-CBF model is estimated using the two-step estimation procedure, and the corresponding estimates are give in Table (ref). Second, since the sample ACFs of each component in Figure (ref) decay slowly, we also fit $\{Y_{t}\}_{t=1}^{1474}$ by a diagonal VT-CBF-HAR model, and the related estimation results are also listed in Table (ref). From this table, we find that the estimates of the degrees of freedom (especially for $\nu_2$) in both fitted models are close to each other, and both estimates of $\nu_2$ are small indicating the heavy-tailedness of the examined data. For the estimates of the mean parameter matrix $S$, its standard errors based on the VT-CBF model are smaller than those based on the VT-CBF-HAR model. For other estimates of parameter matrices, the estimated diagonal components in each parameter matrix seem to have close values, meaning that the examined three stocks possibly have similar temporal structures. This similarity can also be seen from the values of persistence of each stock in Table (ref), where the persistence of stock $s$ is defined by $\sum_{i=1}^{P}A_{1i,ss}^2+\sum_{j=1}^{Q}B_{1j,ss}^2$ for the VT-CBF model and $A_{(d),ss}^2+A_{(w),ss}^2+A_{(m),ss}^2$ for the VT-CBF-HAR model. After estimation, we then apply our test statistics $\Pi_v(l)$ to both fitted models, and the results summarized in Table (ref) imply that both fitted models are adequate at the 5% level.

Next, we consider the forecasting performance of our proposed diagonal VT-CBF and VT-CBF-HAR models. Specifically, we compute the 1-step, 5-step and 10-step predictions of the RCOV matrices, based on a rolling window procedure with window size equal to $T_0=800$. That is, for $T_0\leq t\leq T-t_0$, we fit models based on $T_0$ observations $\{Y_s\}_{s=t-T_0+1}^{t}$, and forecast $\widehat{Y}_{t+t_0}$ with $t_0=1,5,10$ and calculate the forecasting error by $\widehat{Y}_{t+t_0}-Y_{t+t_0}$. To examine the importance of $\nu_2$ in the CBF models, we also apply the diagonal VT-CAW and VT-CAW-HAR models to do prediction for the purpose of comparison. The diagonal VT-CAW and VT-CAW-HAR models are defined in the same way as the diagonal VT-CAW and VT-CAW-HAR models, except that the matrix-F distribution for $\Delta_t$ in the latter two models is replaced by the Wishart distribution. Besides the CAW-type models, we further include a diagonal VAR-HAR model for comparison, where this VAR model uses an HAR structure with the diagonal autoregressive parameter matrices to fit $y_{t}=vech(Y_{t})$.

Table (ref) gives the average of forecasting errors in Frobenius and spectral norms for all models. From this table, we can find that regardless of the prediction horizon, the diagonal VT-CBF-HAR model always has the smallest forecasting error in both norms. Moreover, we apply the DM test (Diebold and Mariano, 1995) to examine whether the diagonal VT-CBF-HAR model has a significant forecasting accuracy over other four competing models. The corresponding testing results are given in Table (ref), and they show that the VT-CBF-HAR model is significantly better than its four competing models in terms of 5-step and 10-step forecasts. For 1-step forecasts, the VT-CBF-HAR and VT-CBF model models have comparable forecasting accuracy, and the VT-CBF-HAR model is significantly better than the remaining three models at level 10%. Note that the VAR-HAR model always performs worst in all examined cases, and this is probably because the VAR-HAR model brutally disentangles the matrix-structure of the RCOV matrices, which may have some intrinsic and useful value for forecasts.

Application 2

In this section, we consider intra-day data of 112 stocks from four major sectors constituting S&P 500 index: 31 stocks from financial sector, 31 stocks from industrial sector, 25 stocks from health care sector, and 25 stocks from consumer discretionary sector, see Table (ref). All intra-day price data are downloaded from Wharton Research Data Services (WRDS) database, and they are taken from 1 July, 2009 to 30 December, 2016, including 1890 non-missing dates of trading data in total.

Based on 100 times log of the price data, the daily RCOV matrices $\{Y_{t}\}_{t=1}^{1890}$ are calculated by the TARVM method in Tao et al. (2011) for each sector.

For each sector, since the dimension of the RCOV matrix is large, we fit the RCOV matrix data by the diagonal F-VT-CBF and F-VT-CBF-HAR models. To do this, we first look for the value of $r$ in model ((ref)) by plotting the ratios $\{\frac{\lambda_{i}}{\lambda_{i+1}}\}$ for each sector in Fig\,(ref), where $\{\lambda_{i}\}$ are the eigenvalues of $\bar{S}$ in descending order. From Fig\,(ref), we can choose $r=3$ for financial sector, $r=2$ for industrial sector, $r=2$ for health care sector, and $r=1$ for consumer discretionary sector. To get more information, we also plot the ratios $\{\frac{\lambda_{i}}{\lambda_{i+1}}\}$ for all four pooled sectors in Fig\,(ref), from which $r=3$ is suggested. This implies that all 112 stocks considered may be driven by 3 latent factors, but among which only two may affect the industrial and health care sectors, and only one may affect the consumer discretionary sector. Hence, it is more reasonable to study the RCOV matrix data across sectors rather than together.

Next, we estimate the diagonal F-VT-CBF and F-VT-CBF-HAR models and choose the orders by a similar procedure as in Application 1, and the related results are reported in Table (ref). From this table, we can find that except for the mean parameter matrix, the diagonal components of other parameter matrices seem to have different values, meaning that each component of $Y_{ft}$ has a different dynamical structure. Moreover, the values of persistence for $Y_{ft,ss}$ show clear differences across four sectors, with the largest persistence in financial sector and the smallest persistence in health care sector. This finding indicates that the effect of past stock returns to its current volatility decays very slowly in the financial sector, while it behaves oppositely in the health care sector.

In the end, we examine the forecasting performance of our F-CBF models. As in Application 1, five different diagonal factor models (see Table (ref)) are considered to forecast $Y_t$, based on a rolling window procedure with window size equal to $1000$. Their forecasting performance is evaluated by the average of forecasting errors in Frobenius and spectral norms as well as the results of the related DM test in Table (ref). From this table, we can see that except for the health care sector, the diagonal F-VT-CBF-HAR model always has the smallest forecasting error and the diagonal F-VAR-HAR model has the largest forecasting error. For 1-step forecasts in the health care sector, the diagonal F-VT-CAW-HAR has slightly smaller forecasting error compared with the diagonal F-VT-CBF-HAR model. In view of the results of DM test, the diagonal F-VT-CBF-HAR model has a significantly better performance than the other four competing models in terms of 5-step and 10-step forecasts, but this advantage is slightly weak in terms of 1-step forecasts, for which the diagonal F-VT-CBF and F-VT-CAW-HAR models have similar performance in the industrial sector, and the diagonal F-VT-CAW-HAR and F-VAR-HAR models have comparative performance in the health care sector.

Concluding Remarks

This paper proposes a new CBF model to study the dynamics of the RCOV matrix. For this CBF model, we explore its stationarity and moment properties, establish the asymptotics of its maximum likelihood estimator, and investigate the inner-product-based tests for its model checking. Hence, a systematic inferential tool of this CBF model is available for empirical researchers. In order to deal with large dimensional RCOV matrices, we also construct two reduced CBF models: the VT-CBF model and the F-CBF model. For both reduced models, the asymptotic theory of the estimated parameters is derived. Compared with the CAW model with Wishart innovations, the CBF model with matrix-F innovations is more able in capturing the heavy-tailed RCOV. This advantage is demonstrated by two real examples on U.S. stock markets. As motivated by Chiriac and Voev (2011), one obvious future work is to introduce the fractional integration structure into our CBF models. Another interesting potential future work could extend the idea of using the matrix-F innovation in a number of ways resulting in a large family of models, which shall be important to study the positive definite dynamics.