Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
74,762 characters · 13 sections · 52 citation commands
Matrix GARCH Model: Inference and Application
\if11 { \affil[1]{\it Center for Statistical Science, Tsinghua University, Beijing 100084, China} \affil[2]{\it Department of Statistics and Data Science, School of Management, Fudan University, Shanghai, China} \affil[3]{\it Department of Statistics and Actuarial Science, University of Hong Kong, Hong Kong} \setcounter{Maxaffil}{0}
} \fi \if01 {
} \fi
\setcounter{Maxaffil}{0}
Keywords: matrix factor GARCH model; matrix GARCH model; matrix time series; portmanteau test; quasi maximum likelihood estimator.
The increasing complexity of time series data has led to the collection of data that extend beyond the conventional vectorial mode. When time series data are observed at the intersections of two classifications, they naturally form a matrix structure. Examples of such data include economic indicator data chen2021autoregressive, international trade flow data chen2022tradenetwork, global stock market data chen2023testing, and many others. To study the conditional mean of matrix time series, chen2021autoregressive propose a matrix autoregressive (MAR) model. The estimation of MAR model is tractable when the dimension of matrix time series is small, however, its implementation tends to be computationally infeasible when the dimension of matrix time series is large. To handle large dimensional matrix time series, Xiao2022reduced propose a reduced rank MAR model to achieve dimension reduction. Their model is distinct from matrix factor models, which aim to learn a much smaller dimensional latent matrix factor that is feasible for modelling or prediction. For a burgeoning literature on matrix factor models, one can refer to Wang2019Factor, chen2019constrained, chen2021statistical, YU2022projected, zhang2022modeling, and chang2021modelling.
To date, the aforementioned works focus solely on modelling conditional mean of matrix time series, without accounting for a prevalent fact that many economic and financial time series data exhibit conditional heteroskedasticity (engle1982autoregressive). As an illustrating example, Figure (ref) plots a time series list of tables recording the sector indices in various international stock markets. The returns of these sector indices naturally form an $m\times n$ matrix time series, where $m$ and $n$ are the numbers of sectors (rows) and markets (columns), respectively. For this return matrix time series, we find that the sample autocorrelations of each entry series are insignificant, suggesting the absence of dynamic structure in its conditional mean. However, we observe that the sample autocorrelations of each squared entry series are very significant. This is a clear indication of the conditional heteroskedasticity, a characteristic cannot be captured by matrix time series models that only account for the conditional mean.
Since the seminal works of engle1982autoregressive and bollerslev1986generalized, the generalized autoregressive conditional heteroskedasticity (GARCH) model and its variants have become benchmark/standard tools for studying the conditional heteroskedasticity of time series. Along this line, one can transform an $m\times n$ matrix time series into an $mn$-dimensional vector time series by the vectorization operator, and then investigate the conditional heteroskedasticity of this vector time series via multivariate GARCH-type models, such as the constant conditional correlation (CCC) model (bollerslev1990modelling), the BEKK model (engle1995multivariate), the dynamic conditional correlation (DCC) model (engle2002dynamic), and many others. See bauwens2006multivariate and francq2019garch for comprehensive surveys on multivariate GARCH-type models. However, due to such a vectorization operation, the multivariate GARCH-type models have two major drawbacks to study the conditional heteroskedasticity of matrix time series. First, they bring a damage to the original matrix data structure, making it difficult to informatively interpret the conditional row and column covariance matrices. Second, they need to estimate model parameters of order $O(m^2n^2)$, so the related estimation becomes computationally challenging even for the moderate values of $m$ and $n$.
In this paper, we first propose a new matrix GARCH model to tackle the aforementioned limitations. Our matrix GARCH model applies two different BEKK specifications to capture the conditional row and column covariance matrices, respectively, and at the same time, it adopts a univariate GARCH specification to depict the trace of conditional row or column covariance matrix. Formally, the matrix GARCH model is related to the existing BEKK and univariate GARCH models, and its trace specification is essential to identifying the conditional row and column covariance matrices. Compared with the aforementioned multivariate GARCH-type models, the matrix GARCH model offers the advantage of directly interpreting the conditional row and column covariance matrices, while also reducing computational complexity with fewer parameters of order $O(m^2+n^2)$.
Next, in order to provide statistical inference for the proposed matrix GARCH model, we develop a quasi maximum likelihood estimator (QMLE) of its parameters and a portmanteau test for its diagnostic checking. Under certain conditions, the asymptotics of the QMLE and portmanteau test are established. Further, we introduce a new matrix factor GARCH model to handle the situations when either $m$ or $n$ is large. Finally, we assess the finite-sample performance of the QMLE and portmanteau tests by simulation studies, and provide three applications on credit default swap prices, global stock sector indices, and future prices to demonstrate the advantage of our matrix GARCH and matrix factor GARCH models over existing multivariate GARCH-type models in volatility forecasting and portfolio allocation.
The remaining paper is organized as follows. Section (ref) introduces the specification of matrix GARCH model. Section (ref) studies the QMLE with its asymptotics for the model. Section (ref) proposes the portmanteau test for the model diagnostic checking. The matrix factor GARCH model is formed in Section (ref). Simulation results are reported in Section (ref). Three applications are presented in Section (ref). Concluding remarks are offered in Section (ref). The proofs of all theorems are deferred to the supplementary materials.
Throughout the paper, $\mathbb{R}$ is a real line, $\mathbf{I}_{p}$ is the $p\times p$ identity matrix, $\mathbf{1}_p$ is the $p$-dimensional vector of ones; “$\stackrel{p}{\to}$” and “$\stackrel{d}{\to}$” denote the convergence in probability and in distribution, respectively. For a square matrix $\mathbf{X}\in \mathbb{R}^{p\times p}$, $\text{tr}(\mathbf{X})$, $|\mathbf{X}|$, and $\rho(\mathbf{X})$ are its trace, determinant, and spectral radius, respectively. For a rectangular matrix $\mathbf{X}\in \mathbb{R}^{p\times q}$, $\mathbf{X}'$, $\text{vec}(\mathbf{X})$, $\|\mathbf{X}\|:=\sqrt{\mathrm{tr}(\mathbf{X}'\mathbf{X})}$, and $\|\mathbf{X}\|_{spec}:=\sqrt{\rho(\mathbf{X}'\mathbf{X})}$ are its transpose, vectorization, Frobenius norm, and spectral norm, respectively.
Consider a time series of length $T$, in which at each time point $t$, an $m \times n$ matrix $\mathbf{X}_t$ is observed. Let $\mathcal{F}_{t}=\sigma(\mathbf{X}_s; s\leq t)$ be the natural filtration containing all available information up to time $t$. To accommodate the matrix structure of $\mathbf{X}_t$, we study its dynamics of conditional heteroskedasticity via the following matrix time series model
where $\mathbf{U}_t\in\mathcal{F}_{t-1}$ and $\mathbf{V}_t\in\mathcal{F}_{t-1}$ are $m\times m$ and $n\times n$ positive definite random matrices, respectively, and $\{\mathbf{Z}_t\}_{t=1}^{T}$ is a sequence of independent and identically distributed (i.i.d.) $m\times n$ random matrix innovations satisfying $E(\mathbf{Z}_t)=\mathbf{0}$ and $E\left[\text{vec}(\mathbf{Z}_t) \text{vec}(\mathbf{Z}_t)'\right]=\mathbf{I}_{mn}$.
Clearly, model (ref) is motivated by the conventional multivariate GARCH models. The two elements, namely, $\mathbf{U}_t$ and $\mathbf{V}_t$, are key to studying the conditional row and column covariance matrices of $\mathbf{X}_t$. In fact, the conditional row covariance matrix of $\mathbf{X}_t$ in ((ref)) is
by noting that $E\left(\mathbf{Z}_t\mathbf{V}_t\mathbf{Z}_t'|\mathcal{F}_{t-1}\right)=\text{tr}(\mathbf{V}_t)\mathbf{I}_m$. Similarly, the conditional column covariance matrix of $\mathbf{X}_t$ in ((ref)) is
The results of (ref) and (ref) demonstrate that the conditional row and column covariances of $\mathbf{X}_t$ are entangled through $\mathbf{U}_t$ and $\mathbf{V}_t$. However, since $\text{tr}(\mathbf{V}_t)\mathbf{U}_t=\text{tr}(a\mathbf{V}_t)[\mathbf{U}_t/a]$ for any $a\neq0$, $\mathbf{U}_t$ and $\mathbf{V}_t$ can not be identified without a constraint on $\text{tr}(\mathbf{U}_t)$ or $\text{tr}(\mathbf{V}_t)$.
To tackle the identification issue of $\mathbf{U}_t$ and $\mathbf{V}_t$, we model $\mathbf{U}_t$ and $\mathbf{V}_t$ by
where $\mathbf{S}_{1,t}\in\mathcal{F}_{t-1}$ and $\mathbf{S}_{2,t}\in\mathcal{F}_{t-1}$ are $m\times m$ and $n\times n$ positive definite random matrices, respectively, $y_t\in\mathcal{F}_{t-1}$ is a positive random variable, and they satisfy
Here, $\mathbf{A}_0$, $\mathbf{A}_1$, and $\mathbf{A}_2$ are $m\times m$ parameter matrices, $\mathbf{B}_0$, $\mathbf{B}_1$, and $\mathbf{B}_2$ are $n\times n$ parameter matrices, $\mathbf{A}_0$ and $\mathbf{B}_0$ are lower triangular matrices with non-negative diagonal elements and
and the scale parameters $w > 0$, $\alpha \ge 0$, and $\beta \ge 0$ ensure the positivity of $y_t$.
The way we model $\mathbf{U}_t$ and $\mathbf{V}_t$ in ((ref))--((ref)) is new to the literature, and it is well suited to the matrix nature of $\mathbf{X}_t$. In ((ref))--((ref)), both $\mathbf{S}_{1,t}$ and $\mathbf{S}_{2,t}$ have the BEKK specifications (engle1995multivariate) to capture the dynamics of conditional row and column covariance matrices, respectively. However, due to the aforementioned identification issue, $\mathbf{U}_t$ and $\mathbf{V}_t$ are not directly related to $\mathbf{S}_{1,t}$ and $\mathbf{S}_{2,t}$ directly, but rather to their normalized versions $\mathbf{S}_{1,t}/\text{tr}(\mathbf{S}_{1,t})$ and $\mathbf{S}_{2,t}/\text{tr}(\mathbf{S}_{2,t})$. The normalization of $\mathbf{S}_{2,t}$ in ((ref)) guarantees $\text{tr}(\mathbf{V}_t)=1$, so that the identification problem can be solved. The normalization of $\mathbf{S}_{1,t}$ in ((ref)) ensures that $\text{tr}(\mathbf{U}_t)=y_t$, and together with the condition $\text{tr}(\mathbf{V}_t)=1$, we have
The above result implies that under ((ref))--((ref)), $y_t$ is indeed the trace of conditional row or column covariance matrix of $\mathbf{X}_t$, which motivates its univariate GARCH specification (bollerslev1986generalized) in ((ref)). Moreover, since the specifications of $\mathbf{U}_t$ and $\mathbf{V}_t$ are invariant to scale changes of $\mathbf{S}_{1,t}$ and $\mathbf{S}_{2,t}$, we add the constraints in ((ref)) to avoid this ambiguity.
In the sequel, we call model ((ref)) with the specifications of $\mathbf{U}_t$ and $\mathbf{V}_t$ in ((ref))--((ref)) the matrix GARCH model. By ((ref))--((ref)), we know that the conditional row and column covariance matrices of $\mathbf{X}_t$ under the matrix GARCH model are $\mathbf{U}_t$ and $y_t\mathbf{V}_t$, respectively. Intrinsically, the matrix GARCH model is related to univariate or multivariate GARCH-type models. For example, when $m=n=1$, we have $\mathbf{U}_{t} = y_t$, $\mathbf{V}_{t}= 1$, and $\mathbf{X}_t = y_t\mathbf{Z}_t$, so the matrix GARCH model reduces to the univariate GARCH model in bollerslev1986generalized. When $n=1$, $m>1$, and $y_t/\text{tr}(\mathbf{S}_{1,t}) = c$ for some constant $c>0$, we have $\mathbf{U}_t=c\mathbf{S}_{1,t}$, $\mathbf{V}_t=1$, and $\mathbf{X}_t=c\mathbf{S}_{1,t}\mathbf{Z}_t$, so the matrix GARCH model becomes the BEKK model in engle1995multivariate. Note that the condition $y_t/\text{tr}(\mathbf{S}_{1,t}) = c$ holds, when $\mathbf{A}_1 = a_1\mathbf{I}_m$, $\mathbf{A}_2 = a_2\mathbf{I}_m$, $w = c\|\mathbf{A}_0 \|^2$, $\alpha = c(a_1m)^2$, and $\beta = c(a_2m)^2$.
Using the vectorization operator, the matrix GARCH model can also be represented as
where $\mathbf{\Sigma}_t$ is the $mn \times mn$ conditional covariance matrix of $\text{vec}(\mathbf{X}_t)$ with the Kronecker product structure
However, from ((ref))--((ref)), we see that except for some special cases (as demonstrated above), the matrix GARCH model in general does not have the form of the BEKK model, for which $\mathbf{\Sigma}_t$ has the specification
where $\mathbf{C}_0$, $\mathbf{C}_1$, and $\mathbf{C}_2$ are $mn\times mn$ parameter matrices defined similarly as those in ((ref)). One may use the BEKK model in ((ref)) (or other multivariate GARCH-type models) to model $\text{vec}(\mathbf{X}_t)$ directly. However, this approach has two disadvantages. First, the matrix structure of $\mathbf{X}_t$ is brutally broken, so the fitting results from the BEKK model lack interpretations on conditional row or column covariance. Second, the number of parameters in the BEKK model has an order of $O(m^2n^2)$, whereas that in the matrix GARCH model has a much smaller order of $O(m^2+n^2)$, even for the moderate values of $m$ and $n$.
Needless to say, the matrix GARCH model can be easily extended to the case where $\mathbf{S}_{1,t}$ and $\mathbf{S}_{2,t}$ have higher order BEKK specifications, and $y_t$ has a higher order GARCH specification. With slight modifications, our proposed methodologies and their related technical treatments below are available for this extended higher order model, which is thus not investigated in details for ease of exposition.
This section studies the QMLE of our matrix GARCH model. Let $\theta = (\gamma', \delta', \zeta')' \in \Theta$ be a vector of unknown parameters with the true value $\theta_0$ for the matrix GARCH model, where $\gamma=(w, \alpha, \beta)'\in\Theta_{\gamma}\subset \mathbb{R}^{3}$, $\delta=(\text{vech}(\mathbf{A}_0)',\text{vec}(\mathbf{A}_1)',\text{vec}(\mathbf{A}_2)')'\in\Theta_{\delta}\subset \mathbb{R}^{p_{\delta}}$ with $p_{\delta}=2m^2+[m(m+1)]/2-1$, $\zeta=(\text{vech}(\mathbf{B}_0)', \text{vec}(\mathbf{B}_1)', \text{vec}(\mathbf{B}_2)')'\in\Theta_{\zeta}\subset \mathbb{R}^{p_{\zeta}}$ with $p_{\zeta}=2n^2+[n(n+1)]/2-1$, $\Theta=\Theta_{\gamma}\times\Theta_{\delta}\times\Theta_{\zeta}\subset\mathbb{R}^{p}$ is a compact parameter space with $p=3+p_{\delta}+p_{\zeta}$, and $\theta_0\in\Theta$. Here, due to the constraints in ((ref)), we exclude $\mathbf{A}_{0,11}$ and $\mathbf{B}_{0,11}$ in $\text{vech}(\mathbf{A}_0)$ and $\text{vech}(\mathbf{B}_0)$, respectively.
By assuming that $\mathbf{Z}_t$ follows the standard matrix normal distribution $\text{MN}(\mathbf{0}, \mathbf{I}_m, \mathbf{I}_n)$ (see its definition in Chapter 2 of gupta1999matrix), we have $\text{vec}(\mathbf{X}_t) | \mathcal{F}_{t-1} \sim \text{N}(\mathbf{0}, \mathbf{\Sigma}_t)$ with $\mathbf{\Sigma}_t = \mathbf{V}_t \otimes \mathbf{U}_t$. Then the quasi log-likelihood function (ignoring constants) of $\{\text{vec}(\mathbf{X}_t)\}_{t=1}^{T}$ is
where $$\mathbf{\Sigma}_t(\theta)= \frac{y_t(\gamma)}{\text{tr}(\mathbf{S}_{1,t}(\delta))\text{tr}(\mathbf{S}_{2,t}(\zeta))} \mathbf{S}_{2,t}(\zeta) \otimes \mathbf{S}_{1,t}(\delta)$$ with $\mathbf{S}_{1,t}(\delta)$, $\mathbf{S}_{2,t}(\zeta)$, and $y_t(\gamma)$ being calculated recursively by
To compute $L_T(\theta)$, the initial values $\{ \mathbf{X}_t: t \leq 0\}$ are necessary, which are unfortunately unobservable. Thus, we consider the computationally feasible version of $L_T(\theta)$ as follows
where $\widetilde{l}_t(\theta)$ is defined in the same way as $l_t(\theta)$ with $\mathbf{\Sigma}_t(\theta)$ being replaced by $\widetilde{\mathbf{\Sigma}}_t(\theta)$, and $\widetilde{\mathbf{\Sigma}}_t(\theta)$ is calculated in the same way as $\mathbf{\Sigma}_t(\theta)$ based on the given initial values $\mathbf{X}_0=\mathbf{0}$, $\mathbf{S}_{1,0}=\mathbf{0}$, $\mathbf{S}_{2,0}=\mathbf{0}$, and $y_{0}=0$.
Let $\widehat{\theta}$ be the minimizer of $\widetilde{L}_T(\theta)$, that is,
Although $\widehat{\theta}$ is defined by assuming $\mathbf{Z}_t\sim\text{MN}(\mathbf{0}, \mathbf{I}_m, \mathbf{I}_n)$, the consistency and asymptotic normality of $\widehat{\theta}$ below still hold for general distributions of $\mathbf{Z}_t$ with some moment constraints. Thus, we call $\widehat{\theta}$ the QMLE of $\theta_0$.
To study the asymptotic properties of $\widehat{\theta}$, some technical assumptions are needed.
Assumption (ref) is commonly used for statistical inference on nonlinear time series models. For multivariate GARCH-type models, the conditions for their stationarity and ergodicity are obtained by using the stochastic recurrence equation theory (Ling2003asymptotic,Hafner2009asymptotic,boussama2011stationarity). However, such an approach seems inapplicable for the matrix GARCH model, since its $\mathbf{\Sigma}_t$ does not have a tractable recursive form. Assumption (ref) is the standard identification condition originated from multivariate GARCH-type models; see comte2003asymptotic, Hafner2009asymptotic, and Zhou2021CBF. Assumption (ref) gives some technical conditions to bound $E\|\mathbf{\Sigma}_t^{-1}(\theta) \|$ and is used to deal with the initial value problem. Assumption (ref)(i) is similar to sufficient conditions for the covariance stationarity of the BEKK or univariate GARCH models, and Assumption (ref)(ii) is mild as it holds trivially when $s\rightarrow 1$.
Now, we are ready to establish the consistency and asymptotic normality of $\widehat{\theta}$.
In the two aforementioned theorems, the moment conditions of $\mathbf{X}_t$ are deemed sufficient. However, we are of the opinion that a more intricate analysis could potentially aid in reducing these requirements. Using $\widehat{\theta}$, we can estimate $\mathcal{C}_0$ and $\mathcal{C}_1$ in Theorem (ref) by
respectively. Under the conditions of Theorem (ref), it is not hard to see that $\widehat{\mathcal{C}}_0$ and $\widehat{\mathcal{C}}_1$ are consistent estimators of $\mathcal{C}_0$ and $\mathcal{C}_1$, respectively. Based on $\widehat{\mathcal{C}}_0$ and $\widehat{\mathcal{C}}_1$, the asymptotic standard errors of $\widehat{\theta}$ can be easily computed.
Diagnostic checking is a crucial step in time series analysis. For multivariate GARCH-type models, several diagnostic checking methods have been proposed in the literature; see, for example, ling1997diagnostic, tse2002residual, and escanciano2013automatic. In this section, we follow the idea of ling1997diagnostic and then propose a new portmanteau test for checking the adequacy of a fitted matrix GARCH model.
Let $\widehat{\mathbf{Z}}_t = \widehat{\mathbf{U}}_t^{-1/2} \mathbf{X}_t\widehat{\mathbf{V}}_t^{-1/2}$ be the residual of a fitted matrix GARCH model, where $\widehat{\mathbf{U}}_t$ and $\widehat{\mathbf{V}}_t$ are computed based on the QMLE $\widehat{\theta}$. Note that if the matrix GARCH model is correctly specified, $\{\text{vec}(\mathbf{Z}_t)\}_{t=1}^{T}$ is a sequence of i.i.d. random vectors with $E\left[\text{vec}(\mathbf{Z}_t) \text{vec}(\mathbf{Z}_t)'\right] = \mathbf{I}_{mn}$. In this case, it turns out that $\{\text{vec}(\mathbf{Z}_t)'\text{vec}(\mathbf{Z}_t)-mn\}_{t=1}^{T}$ is a sequence of uncorrelated random variables. Hence, if the matrix GARCH model is correctly specified, it is expected that the sample autocorrelation of $\{\text{vec}(\widehat{\mathbf{Z}}_t)'\text{vec}(\widehat{\mathbf{Z}}_t)-mn\}_{t=1}^{T}$ at lag $l(\neq 0)$, denoted as $\widehat{R}_l$, is close to zero, where
with $\widehat{\mathbf{\Sigma}}_t=\widehat{\mathbf{V}}_t\otimes\widehat{\mathbf{U}}_t$.
For a given integer $L \ge 1$, let $\widehat{R} = (\widehat{R}_1, \widehat{R}_2, ..., \widehat{R}_L)'\in\mathbb{R}^{L}$, $\mathbf{M} = (M_1, M_2, ..., M_L)'\in\mathbb{R}^{L\times p}$ and $\mathbf{N} = (N_1, N_2, ..., N_L)'\in\mathbb{R}^{L \times p}$, where
Moreover, let $\kappa = E\big( \sum_{i =1}^{m}\sum_{j=1}^{n}\mathbf{Z}_{t,ij}^2 \big)^2 - m^2n^2$ and $\boldsymbol{\eta}=(\eta_{11},\eta_{12},...,\eta_{mn})'\in\mathbb{R}^{mn}$ with $\eta_{ij} = E(\mathbf{Z}_{t,ij}^4) - 1$. To facilitate our portmanteau test, we need the limiting distribution of $\widehat{R}$ below.
Based on Theorem (ref), we construct the portmanteau test statistic at the lag $L$ as
where $\widehat{\mathbf{\Omega}}$ is the sample counterpart of $\mathbf{\Omega}$. If the matrix GARCH model is correctly specified, then $Q_T(L) \stackrel{d}{\rightarrow} \chi_{L}^2$ as $T \rightarrow \infty$. Thus, when the value of $Q_T(L)$ is larger than the $\alpha^*$-th upper-tailed critical value of $\chi_{L}^2$, the fitted matrix GARCH model is inadequate at the significance level $\alpha^*$. Otherwise, it is adequate. In practice, the choice of $L$ depends on the frequency of the series, and one can often choose $L$ to be $O(\log(T))$, which delivers 6 or 8 for a moderate value of $T$ (see tsay2008).
Considering that the number of parameters in the matrix GARCH model has an order of $O(m^2 + n^2)$, the estimation would be computationally challenging when the dimension $m$ or $n$ is large. To deal with the matrix data with large dimension, we adopt the idea of matrix factor model (Wang2019Factor) by assuming that $\mathbf{X}_t$ is driven by a $k_1 \times k_2$ latent common factor matrix $\mathbf{F}_t$, with $k_1$ and $k_2$ much smaller than $m$ and $n$. In particular, we assume
where $\mathbf{R}^*$ is an $m \times k_1$ row factor loading matrix exploiting the variations of $\mathbf{X}_t$ across the rows, $\mathbf{C}^*$ is an $n \times k_2$ column factor loading matrix reflecting the variations of $\mathbf{X}_t$ across the columns, $\mathbf{E}_t$ is an $m\times n$ idiosyncratic matrix white noise error with $E(\mathbf{E}_t)=0$ and $\mathbf{\Sigma}_e = E\left[\text{vec}(\mathbf{E}_t) \text{vec}(\mathbf{E}_t)' \right]<\infty$, and $\mathbf{E}_t$ is also uncorrelated with $\mathbf{F}_{t}^*$. For simplicity, we assume that $(\mathbf{R}^*)'\mathbf{R}^*=\mathbf{I}_{k_1}$ and $(\mathbf{C}^*)'\mathbf{C}^*=\mathbf{I}_{k_2}$.
For model ((ref)), YU2022projected propose the projection estimators $\widehat{\mathbf{R}}^*$ and $\widehat{\mathbf{C}}^*$, and under some regular conditions, they show that
for two orthogonal matrices $\mathbf{H}_{1}\in\mathcal{R}^{k_1\times k_1}$ and $\mathbf{H}_{2}\in\mathcal{R}^{k_2\times k_2}$. See also Wang2019Factor and chen2021statistical for other estimators of loading matrices. Due to the factor nature of model ((ref)), instead of considering the estimated loading matrices $\widehat{\mathbf{R}}^*$ and $\widehat{\mathbf{C}}^*$ directly, we follow chen2022factor to consider their varimax rotated versions (kaiser1958varimax) given by $\widehat{\mathbf{R}}^*\mathbf{G}_1$ and $\widehat{\mathbf{C}}^*\mathbf{G}_2$, where $\mathbf{G}_1$ and $\mathbf{G}_2$ are $k_1\times k_1$ and $k_2\times k_2$ orthogonal matrices, respectively. We focus on the varimax rotated loading matrices, because varimax rotation brings simpler structure so that each factor will have fewer variables with large loadings, leading to a better interpretation of factors.
The consistency results in ((ref)) and varimax rotation motivate us to reshape model ((ref)) as
where $\mathbf{R}=\mathbf{R}^*\mathbf{H}_{1}\mathbf{G}_1$, $\mathbf{C}=\mathbf{C}^*\mathbf{H}_{2}\mathbf{G}_2$, and $\mathbf{F}_t=\mathbf{G}_1'\mathbf{H}_{1}'\mathbf{F}_t^*\mathbf{H}_{2}'\mathbf{G}_2'$ that is assumed to follow the matrix GARCH model: $\mathbf{F}_t=\mathbf{U}_{f,t}^{1/2}\mathbf{Z}_t\mathbf{V}_{f,t}^{1/2}$. Here, $\mathbf{U}_{f,t}$ and $\mathbf{V}_{f,t}$ are defined in the same way as $\mathbf{U}_{t}$ and $\mathbf{V}_{t}$ in ((ref)) with $\mathbf{X}_t$ replaced by $\mathbf{F}_t$, and $\mathbf{Z}_t$ is independent of $\mathcal{F}_{t-1}$. We call model ((ref)) the matrix factor GARCH model, since $\mathbf{F}_t$ has the matrix GARCH structure.
Using the matrix factor GARCH model, we are able to study the high-dimensional conditional covariance matrix of $\text{vec}(\mathbf{X}_t)$, which is useful for portfolio selections based on many assets that have a matrix time series structure in nature (see our two applications in Section (ref) below). Specifically, we have $\text{vec}(\mathbf{X}_t) = (\mathbf{C}\otimes \mathbf{R})\text{vec}(\mathbf{F}_t) + \text{vec}(\mathbf{E}_t)$ under ((ref)), so the conditional covariance matrix of $\text{vec}(\mathbf{X}_t)$ is
where $\mathbf{\Sigma}_{f,t}=\mathbf{V}_{f,t}\otimes \mathbf{U}_{f,t}$. Let $\widehat{\mathbf{R}}^{**}=\widehat{\mathbf{R}}^*\mathbf{G}_1$, $\widehat{\mathbf{C}}^{**}=\widehat{\mathbf{C}}^*\mathbf{G}_2$, and $\widehat{\mathbf{F}}_t = (\widehat{\mathbf{R}}^{**})' \mathbf{X}_t \widehat{\mathbf{C}}^{**}$. In view of ((ref)), we can estimate $\mathbf{\Sigma}_{x,t}$ by
where $\widehat{\mathbf{V}}_{f,t}$ and $\widehat{\mathbf{U}}_{f,t}$ are estimators of $\mathbf{V}_{f,t}$ and $\mathbf{U}_{f,t}$ computed by using the QMLE based on $\{\widehat{\mathbf{F}}_t\}_{t=1}^{T}$, $\widehat{\mathbf{\Sigma}}_e = T^{-1}\sum_{t=1}^{T}\text{vec}(\widehat{\mathbf{E}}_t)\text{vec}(\widehat{\mathbf{E}}_t)'$ with $\widehat{\mathbf{E}}_t = \mathbf{X}_t - \widehat{\mathbf{R}}^{**} \widehat{\mathbf{F}}_t (\widehat{\mathbf{C}}^{**})'$, and the number of factors $k_1$ and $k_2$ is determined by the eigenvalue ratio method in YU2022projected.
In this subsection, we assess the finite sample performance of the QMLE $\widehat{\theta}$. We generate $1000$ replications of sample size $T=1000$ and $2000$ from the following matrix GARCH model
where
with the parameters $w =0.4, \alpha = 0.3, \beta = 0.6$,
and $\{\mathbf{Z}_t\}_{t=1}^{T}$ is an i.i.d. sequence with $\mathbf{Z}_t\sim \text{MN}(\mathbf{0}, \mathbf{I}_3, \mathbf{I}_3)$, $\text{SMT}_{15}$, and $\text{SMT}_{25}$. Here, $\text{SMT}_{\nu}$ is the standardized matrix $t$ distribution with degrees of freedom $\nu$, row covariance $\mathbf{I}_3$, and column covariance $\mathbf{I}_3$.
Based on 1000 replications, Tables (ref), (ref), and (ref) report the bias, sample root mean squared error (SE), and averaged asymptotic standard error (AE) of $\widehat{\theta}$ when $\mathbf{Z}_t$ follows $\text{MN}(\mathbf{0}, \mathbf{I}_3, \mathbf{I}_3)$, $\text{SMT}_{15}$, and $\text{SMT}_{25}$ distributions, respectively. From these tables, we find that (i) regardless of the distribution of $\mathbf{Z}_t$, $\widehat{\theta}$ has small biases, and its values of SE are close to those of AE; (ii) the values of bias and AE decrease as the sample size $T$ increases; (iii) as expected, the value of AE in the case of $\mathbf{Z}_t\sim \text{MN}(\mathbf{0}, \mathbf{I}_3, \mathbf{I}_3)$ is smaller than that in the case of $\mathbf{Z}_t\sim \text{SMT}_{15}$ or $\text{SMT}_{25}$. Overall, our proposed QMLE $\widehat{\theta}$ has a good finite-sample performance in all examined cases.
In this subsection, we exam the finite-sample performance of the portmanteau test $Q_T(L)$. First, we generate 1000 replications of sample size $T=1000$, $2000$, and $4000$ from the following matrix GARCH model
where
with the parameters $\mathbf{A}_0$, $\mathbf{A}_1$, $\mathbf{B}_0$, and $\mathbf{B}_1$ chosen as before and the remaining parameters
and $\{\mathbf{Z}_t\}_{t=1}^{T}$ is an i.i.d. sequence with $\mathbf{Z}_t\sim\text{MN}(\mathbf{0}, \mathbf{I}_m, \mathbf{I}_n)$. Here, we choose the values of $\delta_1$ and $\delta_2$ as follows: Case 1: $\delta_1 = 0.038d$ and $\delta_2 = 0$; Case 2: $\delta_1 = 0$ and $\delta_2 = 0.038d$,\\ where $d = 0, 1, ..., 10$. For each case, we take the null (or alternative) model with respect to $d=0$ (or $d\not=0$).
Next, we fit each replication by the null model, and then apply the portmanteau test $Q_T(L)$ to detect whether the fitted model is adequate. Based on 1000 replications, the empirical power of $Q_T(L)$ across the value of $d$ at the level $\alpha^* = 5\%$ is plotted in Figure (ref), where the lag $L = 2, 4, 6$, and $8$. From Figure (ref), we find that (i) all portmanteau tests have precise sizes; (ii) the power of all portmanteau tests becomes large as the value of $T$ or $d$ increases; (iii) $Q_T(2)$ is generally more powerful than other three portmanteau tests; (iv) all portmanteau tests are more powerful to detect the mis-specification of ARCH part in Case 1 than that of GARCH part in Case 2. In summary, the portmanteau test $Q_T(L)$ performs well especially when $T$ is large.
In this section, we demonstrate the usefulness of the matrix GARCH model and matrix factor GARCH model by three applications. Application 1 applies the matrix GARCH model to study the conditional row and column covariance matrices of a $3\times 3$ matrix time series of Credit Default Swap (CDS) returns. Applications 2 and 3 adopt the matrix factor GARCH model to investigate the portfolio allocations in global stock market and China future market, respectively.
In this application, we consider a $3\times 3$ matrix time series, consisting of the daily log-returns (in percentage) of CDS for three financial institutions: Deutsche Bank (DB), Bank of America (BAC), Barclays Bank (BARC) (with respect to columns) in three tenors of 3 years, 5 years, and 7 years (with respect to rows). The return data for each company and tenor range from $2015/01/01$ to $2018/10/31$, having $T = 1017$ observations in total. As the sample autocorrelations of each entry series are insignificant, we directly fit the demeaned sequence $\{\mathbf{X}_t\}_{t=1}^{1017}$ of this matrix time series (see Figure (ref) for its visualization) using the matrix GARCH model.
Table (ref) reports the QMLE and corresponding standard errors for the fitted matrix GARCH model. The portmanteau tests $Q_T(2)$, $Q_T(4)$, and $Q_T(6)$ have p-values of 0.476, 0.692, and 0.417, respectively, indicating the adequacy of the fitted model. Based on the fitted matrix GARCH model, we plot the estimated conditional row and column covariance matrices in Figures (ref) and (ref), respectively. As expected, we observe from these figures that (i) each entry of both conditional row and column covariance matrices exhibits the clustering phenomenon over time, and (ii) in most of times, each entry of conditional row and column covariance matrices has positive values, indicating that the returns of CDS are generally moving together. Meanwhile, we also make some unexpected but interesting observations from these figures. First, the conditional row covariance matrices generally are much more stable than the conditional column covariance matrices. Second, the conditional row variance of each tenor has several spikes simultaneously around years 2015 and 2016, and similar spikes also largely exist in the conditional row covariances between different tenors. However, this informative phenomenon cannot be found from the conditional column covariance matrices. These two findings indicate that the returns of CDS in the same tenor commonly do not usually exhibit large change for all banks, unless there is a systematic risk in the bank system that causes the CDS price in all tenors to uplift during a very short time period. This kind of systematic risk clearly cannot be captured solely by examining at the conditional covariance of the CDS returns in different tenors across banks.
Next, we apply the matrix GARCH model to forecast the volatility of all nine entry series $\mathbf{X}_{t,ij}$, where $i,j=1, 2, 3$. We choose the first $T_{train}=900$ matrix observations as the training data to estimate model parameters and then fix those estimated parameters to do volatility prediction for the testing data that contain the remaining $T_{test}=100$ matrix observations. In addition, we also compare with the forecasts based on the (first order) univariate GARCH model, BEKK model, column BEKK model, and row BEKK model. Specifically, the univariate GARCH model fits each series $\mathbf{X}_{t,ij}$ independently; the BEKK model fits the $9$-dimensional vector $\textrm{vec}(\mathbf{X}_t)$; the column BEKK model fits each $3$-dimensional column vector $(\mathbf{X}_{t,1j},\mathbf{X}_{t,2j},\mathbf{X}_{t,3j})'$ independently; and the row BEKK model fits each $3$-dimensional row vector $(\mathbf{X}_{t,i1},\mathbf{X}_{t,i2},\mathbf{X}_{t,i3})'$ independently. All the competing models are fitted by the QMLE method, except for the BEKK model that is estimated by the variance targeted QMLE method in Pedersen2014TarBEKK with the diagonal ARCH and GARCH parameter matrices.
After evaluating the forecasting result on the testing data, we summarize the values of mean squared error (MSE), mean absolute error (MAE), and quasi-likelihood (QLIKE) for all considered five models in Table (ref). To determine whether the matrix GARCH model provides significantly more accurate volatility predictions than its competitors, we also perform the DM test (diebold2002comparing). We find from the results of DM test that (i) under the criterions of MSE and MAE, the matrix GARCH model can give more accurate volatility predictions than each competing model, and (ii) under the criterion of QLIKE, the matrix GARCH model does not have a significant prediction advantage over the BEKK and column BEKK models, though it can significantly outperform the two remaining competing models. Overall, the matrix GARCH model has demonstrated better forecasting abilities than its competing models. This is probably because the matrix GARCH model makes use of the complete matrix structure of $\mathbf{X}_t$ while avoiding the inclusion of too many parameters.
In this application, we apply the matrix factor GARCH (MF-GARCH) model to perform global portfolio selections across international stock markets. Specifically, we consider a $10\times 10$ matrix time series, consisting of ten different daily sector indexes (with respect to rows) in ten different stock markets (with respect to columns) from $2017/01/03$ to $2019/11/18$. Here, the sector indexes are Financials, Energy, Industrials, Materials, Information Technology, Telecommunication Service, Utilities, Consumer Staples, Consumer Discretionary, and Healthcare, while the stock markets are SHSZ in China, HSCI in Hong Kong, KOSPI in Korea, SPX in America, SPTSX in Canada, AS in Australia, SBF in France, NZSE in New Zealand, OSEAX in Norway, and SDEURO in Europe. See a visualization of this daily matrix time series in Figure (ref). For this matrix time series, we analyze its daily simple returns (in percentage), and denote this return matrix time series (after demean) by $\{\mathbf{X}_t\}_{t=1}^{700}$ that exhibits no dynamic structure in conditional mean based on its insignificant sample autocorrelations of each entry series. Moreover, we split the entire time period into two parts: the training period containing the first $T_{train}=600$ timepoints and the testing period including the remaining $T_{test}=100$ timepoints. At each timepoint $t_0$ in the testing period, we fit the MF-GARCH model based on the latest $T_{train}$ observations $\{\mathbf{X}_t\}_{t=t_0-T_{train}}^{t_0-1}$, and then use this fitted model to predict the conditional covariance matrix of $\textrm{vec}(\mathbf{X}_{t_0})$ by $\widehat{\mathbf{\Sigma}}_{t_0}$, where $\widehat{\mathbf{\Sigma}}_{t_0}$ is computed as in ((ref)). Note that we choose $k_1=k_2=3$ in all fitted MF-GARCH models, and this choice of $k_1$ and $k_2$ is guided by the eigenvalue ratio method (YU2022projected) based on the training data.
Using $\widehat{\Sigma}_{t_0}$, we select the unconstrained minimum variance portfolio (MVP) from 100 sector indexes by choosing the weight vector as
for $m=n=10$. If the short sales are not preferred or allowed, we can also select the constrained MVP by choosing the weight vector as
where $\widehat{w}_{t_0}^{c}$ has no closed form and need be computed by numerical optimization methods. To evaluate the performance of the proposed unconstrained and constrained MVPs, we follow Engle2019Large to consider the out-of-sample averaged returns (AV), standard deviation of returns (SD), and information ratio (IR) defined by
respectively, where $\widehat{w}_{t_0}$ is $\widehat{w}_{t_0}^{u}$ (or $\widehat{w}_{t_0}^{c}$) with respect to unconstrained (or constrained) MVP, and $R_{t_0}=\mathrm{vec}(\mathbf{X}_{t_0})$. Certainly, many other GARCH-type models can also be used to predict the conditional covariance matrix of $\textrm{vec}(\mathbf{X}_{t_0})$, and the related unconstrained and constrained MVPs can be constructed similarly as done in ((ref))--((ref)). For comparison purpose, we report the results based on the generalized orthogonal GARCH (GO-GARCH) model (van2002go), DCC model (engle2002dynamic), Risk Metrics model (zumbach2007riskmetrics), and Equal Weights model. For the DCC model, its QMLE may not perform well when the dimension of $\textrm{vec}(\mathbf{X}_{t_0})$ is 100. Therefore, we further follow Engle2019Large to estimate the DCC model by using either the linear shrinkage (LS) method in ledoit2004well or nonlinear shrinkage (NLS) method in ledoit2012nonlinear, leading to the DCC-LS and DCC-NLS models for comparison.
Table (ref) reports the annualized AV, SD, and IR of the unconstrained and constrained MVPs proposed by the MF-GARCH model and its six competitors. From this table, we find that the unconstrained and constrained MVPs selected by the MF-GARCH model perform the best in terms of AV, SD, and IR. In particular, the values of IR for the unconstrained and constrained MVPs from the MF-GARCH model are 28% and 20% higher than those from the DCC-NLS model, which has the best performance among all four competing multivariate GARCH-type models. The advantage of the MF-GARCH model over its four multivariate GARCH-type competitors is most likely due to its matrix structure, which can largely alleviate the trouble of dimensionality and reduce the computational burden in model estimation. Moreover, we find that the Risk Metrics model performs the worst. This is expected since a simple weighting parameter used by the Risk Metrics model could not adequately capture the dynamics of $\textrm{vec}(\mathbf{X}_{t})$. Another observation is that the Equal Weights model has an unsatisfactory performance, indicating the necessity of weight selection from the dynamics of conditional covariance matrix of $\textrm{vec}(\mathbf{X}_{t})$.
Last, we plot the out-of-sample cumulative returns of unconstrained and constrained MVPs constructed by all considered models in Figure (ref). From this figure, we find that the unconstrained and constrained MVPs selected by the MF-GARCH model have the best cumulative returns over time. This finding is consistent to the one from Table (ref) that the largest value of AV for each MVP is achieved by the MF-GARCH model.
To further demonstrate the usefulness of the MF-GARCH model, we consider another $8\times 12$ matrix time series that is comprised of weekly closing prices of futures in $8$ varieties (with respect to rows) with $12$ different delivery months (with respect to columns) from $2009/08/19$ to $2022/09/26$, where the varieties include Aluminum, Cuprum, Zinc, Polyethylene, Deformed Steel Bar, Palm Oil, Polyvinyl Chloride, and Pure Terephthalic Acid, and the delivery months are from January to December. See a visualization of this weekly matrix time series in Figure (ref).
Below, we study the weekly simple returns (in percentage) of this matrix time series and denote the corresponding return matrix time series (after demean) by $\{\mathbf{X}_t\}_{t=1}^{637}$. Similar to the implementations in Application 2, we apply the MF-GARCH model and other six competing models to propose the unconstrained and constrained MVPs during the testing period, which contains the last $T_{test}=100$ timepoints. For the MF-GARCH model, we pick up $k_1=2$ and $k_2=3$ according to the eigenvalue ratio method (YU2022projected) based on the first $T_{train}=537$ observations.
Table (ref) reports the annualized AV, SD, and IR of the unconstrained and constrained MVPs proposed by all considered models. From this table, we find that the unconstrained and constrained MVPs selected by the MF-GARCH model have the best values of SD and IR, while those selected by the DCC-LS or DCC-NLS model have the best value of AV. This finding implies that the higher return of MVPs from the DCC-LS and DCC-NLS models is obtained at the price of higher risk, and the MF-GARCH model having the largest value of IR tends to make a better balance between return and risk. This advantage of MF-GARCH model can also be observed from Figure (ref), in which the out-of-sample cumulative returns of two MVPs from the MF-GARCH model seem to be more stable over time than those from other models.
In this paper, we propose a novel matrix GARCH model to study the conditional heteroskedasticity of matrix time series. This matrix GARCH model is comprised of three specifications, including two BEKK specifications (normalized by their own traces) to describe the conditional row and column covariance matrices, and one univariate GARCH specification to capture the trace of conditional row or column covariance matrix. We show that the univariate GARCH specification for the trace is necessary to identify the conditional row and column covariance matrices. Next, we propose the QMLE for the matrix GARCH model and construct the portmanteau test for the model diagnostic checking. Under certain conditions, we establish the asymptotics of the QMLE and portmanteau test. Moreover, we introduce a matrix factor GARCH model to deal with large dimensional matrix time series. Compared with existing multivariate GARCH-type models that can also study the conditional heteroskedasticity of matrix time series after vectorization, the matrix GARCH and matrix factor GARCH models not only have a clear interpretation on conditional row and column covariance matrices, but also significantly reduce the number of model parameters. These two advantages are extensively demonstrated by three real examples.
Our matrix GARCH model opens a new door for studying the conditional heteroskedasticity of complex structured time series. In line with the vast body of literature on multivariate GARCH-type models, one could consider the extension of our matrix-GARCH model to better capture the heavy-tailedness and asymmetry of matrix time series data. Another interesting future work is to investigate the conditional heteroskedasticity of tensor-variant time series (wang2021high,chen2022factor) by using a similar approach as our matrix GARCH model. These future works along with the current work, have the potential to generate valuable toolkits for practitioners to analyze complex time series.
\linespread{1.3}\selectfont