# CNN-based Realized Covariance Matrix Forecasting

Yanwen Fang1 Philip L. H. Yu1,2 Corresponding author. Email address:
<EMAIL_ADDRESS> Yaohua Tang1
1Department of Statistics and Actuarial Science, The University of Hong Kong
2Department of Mathematics and Information Technology, Education University of
Hong Kong

###### Abstract

It is well known that modeling and forecasting realized covariance matrices of
asset returns play a crucial role in the field of finance. The availability of
high frequency intraday data enables the modeling of the realized covariance
matrices directly. However, most of the models available in the literature
depend on strong structural assumptions and they often suffer from the curse
of dimensionality. We propose an end-to-end trainable model built on the CNN
and Convolutional LSTM (ConvLSTM) which does not require to make any
distributional or structural assumption but could handle high-dimensional
realized covariance matrices consistently. The proposed model focuses on local
structures and spatiotemporal correlations. It learns a nonlinear mapping that
connect the historical realized covariance matrices to the future one. Our
empirical studies on synthetic and real-world datasets demonstrate its
excellent forecasting ability compared with several advanced volatility
models.

Keywords: Deep Learning, Convolutional Neural Network, Convolutional LSTM,
Realized Covariance Matrix.

##  1 Introduction

Modeling and forecasting covariance matrices of financial asset returns have
long been an important problem in asset pricing, portfolio allocation and risk
management. The availability of high-frequency intraday financial data enables
us to estimate daily volatilities and co-volatilities of asset returns
directly, which leads to the so-called realized covariance (RCOV) matrices. In
the field of finance, modeling RCOV matrices has always been a hot research
topic.

In this modeling, the first issue is that the fitted RCOV matrices should be
guaranteed to be symmetric and positive definite. A natural choice for this
issue is the family of matrix-valued Wishart distributions which automatically
generates random positive definite matrices without imposing additional
constraints. One example is the conditional autoregressive Wishart (CAW) model
proposed by Golosnoy et al.,  [3]. Another approach is to transform the RCOV
matrices. Bauer and Vorkink,  [1] suggested applying matrix logarithm
transformation, and Wang,  [10] suggested square root transformation to the
RCOV matrices. High-dimensionality is the second issue. For $d$ assets, there
are $d(d+1)/2$ entries in the RCOV matrices, which means the number of
parameters needed to model RCOV matrices grows quickly with $d$.

To overcome the difficulty due to the high-dimensionality, structural
assumptions are needed in order to estimate the RCOV matrix consistently. Some
assumptions of the existing effective models for the RCOV matrices are so
restrictive that they cannot model the dynamic of the RCOV matrices properly.
Tao et al.,  [9] proposed the matrix factor analysis (MFA) model which
projects high dimensional RCOV matrices to a time series of low dimensional
matrices, and then builds a vector autoregressive (VAR) model to the
vectorized time series. Shen et al.,  [7] replaced the VAR model by a diagonal
CAW (DCAW) model with diagonal matrices. As the projected dimension is usually
low, say 3-10, both VAR and DCAW models may contain too few parameters to
determine the complex dependency structure of a high-dimensional RCOV matrix.

With the recent advancement in powerful computing resources (like GPUs [5])
and the easy access to an abundance of data (like ImageNet [2]), deep learning
models have received great attention. They are proved to be excellent in
solving hard problems, like speech recognition, object recognition, and
machine translation where the data are notably complex, and they have won
numerous contests in pattern recognition and machine learning, beating the
stat-of-art traditional methods. Convolutional neural network (CNN) and
recurrent neural network (RNN) are the most common types of neural networks,
which capture image data and sequential data well respectively. Long short-
term memory (LSTM) as a special RNN structure has been proven to perform
extremely well on temporal data [hochreiter1997long, gers2002learning,
sundermeyer2012lstm, gers2002applying, cao2019financial,
srivastava2015unsupervised, sutskever2014sequence]. If a sequence of RCOV
matrices can be viewed as a sequence of maps, the goal of our task is to give
precise prediction of future maps given the previously observed maps, which in
essence is a spatiotemporal sequence forecasting problem. Xingjian et al.,
[12] extended the traditional LSTM models to have convolutional structures,
called ConvLSTM. ConvLSTM performs well on spatiotemporal sequence forecasting
problem which provides some novel insights on our task.

In this paper, we take the local structure of RCOV matrices into consideration
and formulate RCOV matrices forecasting problem as a spatiotemporal sequence
forecasting problem. Base on this, we propose an end-to-end forecasting model
for this problem. The model has almost no assumption on the RCOV matrices and
has a great flexibility that the model could fit to both low-dimensional and
high-dimensional matrices. Unlike other statistical models, the parameter size
of our model does not grow dramatically as the data dimension grows. The
proposed CNN and ConvLSTM based model detects local structures from the
historical RCOV matrices and try to forecast the future RCOV matrix after
applying non-linear mapping to these patterns. We will see later that our
proposed model could obtain a better estimation accuracy in much shorter time
than the existing methods studied.

The rest of this paper is organized as follows. Section 2 gives a brief review
of recent development of time series modeling for RCOV matrices and some
preliminary knowledge about CNN and ConvLSTM. Our model structure is proposed
in Section 3. Simulation results on synthetic data and experimental results on
real data are reported in Section 4. Section 5 concludes with some final
remarks and discusses some future work.

##  2 Preliminaries

###  2.1 Realized Covariance Matrix Modeling

Suppose that there are $d$ assets and their log price process
$\bm{X}_{t}=\left\\{X_{1}(t),\dots X_{d}(t)\right\\}^{\top}$ obeys a
continuous diffusion model

$d\bm{X}(t)=\bm{\mu}_{t}dt+\bm{\sigma}_{t}d\bm{W}_{t},t\in[0,L],$ (1)

where $L$ is an integer, $\bm{W}_{t}$ is a $d$-dimensional standard Brownian
motion, $\bm{\mu}_{t}$ is a drift taking values in $\mathbb{R}^{d}$, and
$\bm{\sigma}_{t}$ is a $d\times d$ matrix. $\bm{\mu}_{t}$ and
$\bm{\sigma}_{t}$ are assumed to be continuous in $t$. Let a day be a unit
time. The integrated volatility matrix for the $t$-th day is defined as

$\bm{\Sigma}_{x}(t)=\int_{t-1}^{t}\bm{\sigma}_{s}\bm{\sigma}_{s}^{\top}ds,~{}~{}~{}~{}t=1,\dots,L.$
(2)

Due to the existence of microstructure noise, we cannot observe $X_{i}(t)$ in
reality. What we observe is $Y_{i}(t_{ij})$, a noisy version of $X_{i}(\cdot)$
at times $t_{ij}\in(t-1,t],j=1,\dots,n_{i},i=1,\dots,d$. Here, $n_{i}$ is the
total trading times and $t_{ij}$ is the $j$-th trading times of asset $i$
during a giving trading day $t$. The observations $Y_{i}(t_{ij})$ are allowed
to be non-synchronized, i.e. $t_{i_{1}j}\neq t_{i_{2}j}$ for any $i_{1}\neq
i_{2}$. In this article, we assume that

$Y_{i}(t_{ij})=X_{i}(t_{ij})+\varepsilon_{i}(t_{ij}),$ (3)

where $\epsilon_{i}(t_{ij})$ are i.i.d. microstructure noise with mean zero
and variance $\eta_{i}$. $\varepsilon_{i}(\cdot)$ and $X_{i}(\cdot)$ are
assumed to be independent with each other. With high-frequency data on the
$t$-th day, we can calculate Threshold Averaging Realized Volatility Matrix
(TARVM) estimator $\bm{\hat{\Sigma}}_{y}(t)$ as an estimator for
$\bm{\Sigma}_{x}(t)$, which is shown to be consistent for the integrated RCOV
matrix under certain conditions [9]. To reduce the effective number of entries
in $\bm{\Sigma}_{x}(t)$, Tao et al.,  [9] proposed a matrix-factor analysis
(MFA) model as follows:

$\bm{\Sigma}_{x}(t)=\bm{A}\bm{\Sigma}_{f}(t)\bm{A}^{\top}+\bm{\Sigma}_{0},~{}~{}~{}~{}~{}t=1,\dots,L,$
(4)

where $\bm{\Sigma}_{0}$ is a $d\times d$ positive definite constant matrix,
$\bm{\Sigma}_{f}(t)$ are $r\times r$ positive definite matrices and $\bm{A}$
is a $d\times r$ factor loading matrix. Here $r(<d)$ is usually a fixed small
integer. (4) assumes that the daily dynamical structure of the matrix process
$\bm{\Sigma}_{x}(t)$ is driven by that of a lower-dimensional latent process
$\bm{\Sigma}_{f}(t)$, where $\bm{\Sigma}_{0}$ represents the static part of
$\bm{\Sigma}_{x}(t)$.

Put

$\displaystyle\bar{\bm{\Sigma}}_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{L}\sum_{t=1}^{L}\bm{\Sigma}_{x}(t),$ (5)
$\displaystyle\bar{\bm{S}}_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{L}\sum_{t=1}^{L}\\{\bm{\Sigma}_{x}(t)-\bar{\bm{\Sigma}}_{x}\\}^{2},$
(6)

and as $\bm{\Sigma}_{x}(t)$ is unknown in practice, we use
$\hat{\bm{\Sigma}}_{y}(t)$ as a proxy. Let

$\displaystyle\bar{\bm{\Sigma}}_{y}$ $\displaystyle=$
$\displaystyle\frac{1}{L}\sum_{t=1}^{L}\hat{\bm{\Sigma}}_{y}(t),$ (7)
$\displaystyle\bar{\bm{S}}_{y}$ $\displaystyle=$
$\displaystyle\frac{1}{L}\sum_{t=1}^{L}\\{\hat{\bm{\Sigma}}_{y}(t)-\bar{\bm{\Sigma}}_{y}\\}^{2},$
(8)

where $\hat{\bm{\Sigma}}_{y}(t)$ is the RCOV matrix on day $t$ obtained by
method TARVM from high-frequency data. Let $\hat{\bm{A}}$ be the $r$
orthonormal eigenvectors of $\bar{\bm{S}}_{y}$, corresponding to the $r$
largest eigenvalues, as its columns. Tao et al.,  [9] proposed the estimated
factor volatilities as

$\hat{\bm{\Sigma}}_{f}(t)=\hat{\bm{A}}^{\top}\hat{\bm{\Sigma}}_{y}(t)\hat{\bm{A}},t=1,\dots,L,$
(9)

and the estimator for $\bm{\Sigma}_{0}$ as

$\hat{\bm{\Sigma}}_{0}=\bar{\bm{\Sigma}}_{y}-\hat{\bm{A}}\hat{\bm{A}}^{\top}\bar{\bm{\Sigma}}_{y}\hat{\bm{A}}\hat{\bm{A}}^{\top}.$
(10)

Then Tao et al.,  [9] built up the dynamical structure of $\bm{\Sigma}_{x}$ by
fitting a VAR model to the vectorized $\hat{\bm{\Sigma}}_{f}(t)$. For a
$r\times r$ matrix $\bm{\Sigma}$, let $vech(\bm{\Sigma})$ be the
$r(r+1)/2\times 1$ vector obtained by stacking together the lower triangle of
$\bm{\Sigma}$ (including the diagonal). Then the VAR model for
$\bm{\Sigma}_{f}(t)$ is of the form

$vech\\{\hat{\bm{\Sigma}}_{f}(t)\\}=\bm{\alpha}_{0}+\sum_{j=1}^{q}\bm{\alpha}_{j}vech\\{\hat{\bm{\Sigma}}_{f}(t-j)\\}+\bm{e}_{t},$
(11)

where $q\geq 1$ is an integer, $\bm{\alpha}_{0}$ is a vector,
$\bm{\alpha}_{1},\dots,\bm{\alpha}_{q}$ are square matrices, and $\bm{e}_{t}$
is a white noise process with zero mean and finite fourth moments.

The VAR model cannot guarantee the positive definiteness of the forecasted
RCOV matrix. To resolve this problem, Shen et al.,  [7] proposed a diagonal
Conditional Autoregressive Wishart (DCAW) model to $\hat{\bm{\Sigma}}_{f}(t)$.
Let $\mathcal{F}_{t-1}=\sigma(\bm{\Sigma}_{f}(s),s<t)$ be the past history of
the process at time $t$. Conditional on $\mathcal{F}_{t-1}$,
$\bm{\Sigma}_{f}(t)$ follows a central Wishart distribution

$\bm{\Sigma}_{f}(t)|\mathcal{F}_{t-1}\sim\mathcal{W}_{n}(\bm{\nu},\bm{S}_{f}(t)/\bm{\nu}),$
(12)

with $\bm{\nu}$ degrees of freedom and the scaling matrix $\bm{S}_{f}(t)$.
Moreover, the scaling matrix $\bm{S}_{f}(t)$ follows a linear BEKK recursion
of order $(p,q)$:

$\begin{split}\bm{S}_{f}(t)=&\bm{C}\bm{C}^{\top}+\sum_{i=1}^{p}\bm{B}_{i}\bm{S}_{f}(t-i)\bm{B}_{i}^{\top}\\\
&+\sum_{j=1}^{q}\bm{A}_{j}\bm{\Sigma}_{f}(t-j)\bm{A}_{j}^{\top},\end{split}$
(13)

where $\bm{A}_{j},\bm{B}_{i}$ and $\bm{C}$ are all $r\times r$ matrices of
coefficients. Different orders $(p,q)$ are used to make comparison among
models.

The CAW process depends on the parameters
$\\{\bm{\nu},\bm{C},(\bm{B}_{i})_{1\leq i\leq p},(\bm{A}_{j})_{1\leq j\leq
q}\\}$, so that the total number of parameters is equal to
$(p+q)r^{2}+\frac{r(r+1)}{2}+1=\mathcal{O}(r^{2})$ which still grows quick
with the number of factors $r$ and the order $p$ and $q$. The DCAW model
restricts the coefficient matrices $\bm{C},(\bm{B}_{i})_{1\leq i\leq
p},(\bm{A}_{j})_{1\leq j\leq q}$ to diagonal matrices. Therefore, the number
of parameters becomes $(p+q+1)r+1=\mathcal{O}(r)$. The estimation of the
parameters
$\bm{\theta}=(\bm{\nu},diag(\bm{C})^{\top},diag(\bm{B}_{i})^{\top}_{1\leq
i\leq p},diag(\bm{A}_{j})^{\top}_{1\leq j\leq q})^{\top}$ of the DCAW($p$,$q$)
model is carried out by maximizing the log-likelihood function using the
Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization procedure.

###  2.2 CNN and ConvLSTM

Convolutional Neural Network (CNN) is a well known variety of deep neural
network which has been widely used to extract local patterns from images. The
convolutional layer aims to learn feature representations of the inputs. It is
composed of several convolution kernels which are used to compute different
feature maps. Specifically, each neuron of a feature map is connected to a
neighborhood of neurons in the previous layer. The new feature map can be
obtained by first convolving the input with a learned kernel and then applying
an element-wise nonlinear activation function on the convolved results. The
activation function introduces nonlinearities to CNN, which are desirable for
multi-layer networks to detect nonlinear features.

Xingjian et al.,  [12] extended the fully connected LSTM (FC-LSTM) to
convolutional LSTM (ConvLSTM) which has convolutional structures in both the
input-to-state and state-to-state transitions for precipitation nowcasting.
Due to the inherent convolutional structure, the ConvLSTM layer is very
suitable for spatiotemporal data. We follow the formulation of ConvLSTM as in
Xingjian et al.,  [12], which use $\mathcal{X}_{1},...,\mathcal{X}_{t}$ as
inputs, $\mathcal{C}_{1},...,\mathcal{C}_{t}$ as cell outputs,
$\mathcal{H}_{1},...,\mathcal{H}_{t}$ as hidden states. $i_{t},f_{t},o_{t}$
are input gate, forget gate and output gate respectively. Then the ConvLSTM is
of the form:

$\displaystyle
i_{t}=\sigma(W_{xi}\otimes\mathcal{X}_{t}+W_{hi}\otimes\mathcal{H}_{t-1}+W_{ci}\circ\mathcal{C}_{t-1}+b_{i})$
(14) $\displaystyle
f_{t}=\sigma(W_{xf}\otimes\mathcal{X}_{t}+W_{hf}\otimes\mathcal{H}_{t-1}+W_{cf}\circ\mathcal{C}_{t-1}+b_{f})$
$\displaystyle\mathcal{C}_{t}=f_{t}\circ\mathcal{C}_{t-1}+i_{t}\circ
tanh(W_{xc}\otimes\mathcal{X}_{t}+W_{hc}\otimes\mathcal{H}_{t-1}+b_{c})$
$\displaystyle
o_{t}=\sigma(W_{xo}\otimes\mathcal{X}_{t}+W_{ho}\otimes\mathcal{H}_{t-1}+W_{co}\circ\mathcal{C}_{t}+b_{o})$
$\displaystyle\mathcal{H}_{t}=o_{t}\circ tanh(\mathcal{C}_{t})$

where $\otimes$ denotes the convolution operation, $\circ$ denotes the
pointwise product. The ConvLSTM captures the local neighbors of the inputs and
past states to determine the future state of a certain cell on a spatial grid
as shown in (14).

##  3 Proposed Model

For simplicity, we use $\bm{\Sigma}_{t}$ to denote $\hat{\bm{\Sigma}}_{y}(t)$
which is viewed as a ground true RCOV matrix. Consider a stochastic, positive
definite RCOV matrix $\bm{\Sigma}_{t}=(\Sigma_{{ij},t})$, with dimension
$d\times d$ at time $t(t=1,\dots,L)$. We now present our proposed model in the
following aspects.

###  3.1 Data Preprocessing

The RCOV matrices are symmetric and positive definite, but the outputs of the
end-to-end network are directly treated as the forecasting RCOV matrices,
which can not guarantee the symmetry and positive definiteness. Therefore,
before we train the network on the training set, data preprocessing methods
must be applied to transform and compress the data.

To make sure the predicted RCOV matrices are positive definite, one method is
to use Cholesky decomposition to transform the original RCOV matrices and then
train the decomposed triangular matrices. The Cholesky decomposition for a
$\mathbf{\Sigma}_{t}$ is of the form

$\mathbf{\Sigma}_{t}=\bm{C}_{t}\bm{C}_{t}^{\top},$ (15)

where $\bm{C}_{t}$ is a lower triangular matrix with real and positive
diagonal entries, and $\bm{C}_{t}^{\top}$ denotes the transpose of
$\bm{C}_{t}$. After we get a predicted $\hat{\bm{C}_{t}}$, we could obtain a
positive definite matrix
$\hat{\bm{\Sigma}}_{t}=\hat{\bm{C}}_{t}\hat{\bm{C}}_{t}^{\top}$. Note that
Cholesky decomposition helps to decrease the range of the data. Such an
operation reduces the large covariance values while boosts the small values,
which makes the data suiting the neural networks better.

Recall that Wang,  [10] used the square root transformation to the RCOV
matrix. Consider the spectral decomposition of the RCOV matrix
$\bm{\Sigma}_{t}=\bm{Q}_{t}\bm{\Lambda}_{t}\bm{Q}_{t}^{\top}$, where
$\bm{\Lambda}_{t}=diag(\lambda_{t,1},\dots,\lambda_{t,p})$ is a diagonal
matrix of the eigenvalues of $\bm{\Sigma}_{t}$, and $\bm{Q}_{t}$ is an
orthonormal matrix consisting of eigenvectors of $\bm{\Sigma}_{t}$. Assume
that $\lambda_{t,1}\geq\lambda_{t,2}\geq\dots\geq\lambda_{t,p}>0$. In our
datasets, we observed the same phenomenon as in Tao et al.,  [8] that the
first several $\lambda_{t,i}$ are much bigger than the others. As the range is
quite large, we found that taking square root of the eigenvalues could reduce
the range. Denote $\bm{M}_{t}$ as
$\bm{M}_{t}=diag(\lambda_{t,1}^{1/2},\dots,\lambda_{t,p}^{1/2})$, and denote
$\bm{O}_{t}=\bm{Q}_{t}\bm{M}_{t}\bm{Q}_{t}^{\top}$, then we have

$\bm{O}_{t}=\bm{\Sigma}_{t}^{1/2}.$ (16)

Both Cholesky decomposition and square root transformation could help maintain
the positive definiteness of the predicted RCOV matrices and reduce the range
of data. These two transformations can be applied together, that is, we may
apply square root transformation to the RCOV matrices and then apply Cholesky
decomposition from the square root transformed matrices.

###  3.2 Encoding-Generating Structure

We name our proposed model Covariance Matrix Convolutional LSTM (CM-ConvLSTM)
model. The proposed model has several appealing properties. First, its
structure is intentionally designed with simplicity in mind, and yet provides
superior accuracy compared with several traditional statistical methods.
Second, our model could handle higher dimensional RCOV matrix but traditional
methods could not. Our model could deal with high-dimensional and low-
dimensional RCOV matrices in a consistent way. Third, with moderate numbers of
filters and layers, our method achieves remarkably fast speed for practical
on-line usage. Our model is easy to implement with the excellent open deep
learning packages.

Given past $m$ observations, a 3D matrix
$\tilde{\bm{\Sigma}}_{t}=[\bm{\Sigma}_{t};\bm{\Sigma}_{t-1};\dots;\bm{\Sigma}_{t-m+1}]$
with dimension $m\times d\times d$, the corresponding transformed 3D matrix
$\tilde{\bm{L}}_{t}=[\bm{L}_{t};\bm{L}_{t-1};\dots;\bm{L}_{t-m+1}]$ can be
obtained using the data preprocessing methods mentioned above. The input of
the network is the transformed matrix $\tilde{\bm{L}}_{t}$, and our model is
an end-to-end model. Therefore, the output of the network is $\bm{L}_{t+1}$.
It is easy to obtain the target prediction $\bm{\Sigma}_{t+1}$ after applying
the same transformations in the data preprocessing stage in the reverse order.
We wish to learn a mapping $F$, which conceptually consists of the following
operations.

####  3.2.1 Encoding

The network consists of two types of encoding layers: ConvLSTM encoding layer
and CNN encoding layer. We first use a ConvLSTM layer to operate on the inputs
which can extract the local information of the inputs and past states at each
timestamp. This layer aims to learn the complex spatiotemporal patterns.

Then we apply one or more CNN layers which work as the local structure
representation and non-linear mapping layer which is equivalent to convolving
the input feature maps by a set of filters. These layers aim to further
capture the local information and conduct non-linear mapping on the features
learnt in the ConvLSTM layer. We apply the Leaky Rectified Linear Unit (LReLU)
as the activation function on the filter responses after each ConvLSTM layer
or CNN layer. Zero-padding is applied to make sure the output feature maps are
the same size as the input transformed RCOV matrices. It is possible to add
more CNN layers to increase the non-linearity according to the complexity of
the data. But this can increase the complexity of the model, and thus demands
more training time.

####  3.2.2 Generating

In the last stage, the predicted overlapping broader local patches are
averaged to produce the final transformed RCOV matrix. The averaging can be
considered as a pre-defined filter on a set of feature maps. Motivated by
this, we define a convolutional layer to generate the transformed RCOV matrix
$\bm{L}_{t+1}$ and so that obtain the final prediction RCOV matrix
$\bm{\Sigma}_{t+1}$. In this layer, instead of using non-linear LReLU
activation function, we use linear mapping here to act like regression on the
learned feature maps in the previous layer. By stacking multiple ConvLSTM
layers and CNN layers, we can build the end-to-end network with the encoding-
generating structure.

An overview of the network is depicted in Figure 1. This shows how the above
operations form a ConvLSTM and CNN based network. Unlike other CNN
architectures, our model has no pooling layers or fully-connected layers. Zero
padding is needed before a convolution operation in the model to ensure that
the ouputs of the network have the same size as the inputs.

Figure 1: Diagram of the proposed encoding-generating model for RCOV matrices
forecasting.

###  3.3 Training

Learning the end-to-end mapping function $F$ requires the estimation of
network parameters. This is achieved through minimizing the loss between the
predicted and the corresponding ground truth matrices. It is worth noting that
deep learning models do not preclude the usage of other kinds of loss
functions, if only the loss functions are derivable. If a better perceptually
motivated metric is given during training, it is flexible for the network to
adapt to that metric.

We were troubled by “outliers” in the RCOV matrices even after data
preprocessing. The network would tend to struggle between the majority of
small values and small proportion of extreme big values that the performances
stuck and are hard to improve. A popular robust loss function called Huber
Loss is chosen instead of L1-norm loss and L2-norm loss to relieve this
situation. Huber Loss is resistant to outliers in the data, while it maintains
sensitive to the small values. And this sensitivity enables this loss function
to perform well on our real-world datasets. Mathematically, the Huber loss
function is defined as

$\ell_{\delta}(r_{i})=\left\\{\begin{array}[]{ll}\frac{1}{2}(r_{i})^{2},&|r_{i}|\leq\delta\\\
\\\ \delta(|r_{i}|-\frac{1}{2}\delta),&otherwise\end{array}\right.$ (17)

where $r_{i}$ is the error of the prediction of the $i$-th observation. Figure
2 distinguishes Huber loss from L1-norm loss and L2-norm loss in terms of
$r_{i}$. When the absolute difference between the observed and predicted
values is larger than a prespecified value $\delta$, the Huber loss is reduced
to the L1-norm loss in order to prevent the squared differences with large
magnitude which may dominate the whole loss.

Figure 2: Huber loss (black solid line), L1-norm loss (red dotted line) and
L2-norm loss (blue dashed line) as a function of $r_{i}$.

As our data is a time series of matrices, the form of Huber loss function in
(17) has to be adapted to the 2D observed and predicted matrices. Given a set
of pairs of $\\{\bm{\hat{Y_{i}}},\bm{Y_{i}}\\}$, where $\bm{\hat{Y_{i}}}$ and
$\bm{Y_{i}}$ are the predicted matrix and the corresponding ground true matrix
of the $i^{th}$ observation, the Huber loss function becomes:

$\ell_{t}(\Theta)=\left\\{\begin{array}[]{ll}\frac{1}{2}||\bm{\hat{Y_{i}}}-\bm{Y_{i}}||^{2},&||\bm{\hat{Y_{i}}}-\bm{Y_{i}}||\leq\delta\\\
\\\
\delta(||\bm{\hat{Y_{i}}}-\bm{Y_{i}}||-\frac{1}{2}\delta),&otherwise\end{array}\right.$
(18)

where $||\bm{\hat{Y_{i}}}-\bm{Y_{i}}||$ and
$||\bm{\hat{Y_{i}}}-\bm{Y_{i}}||^{2}$ are the L1-norm loss and L2-norm loss of
the cells in these two matrices, respectively. Gradient descent method based
on Adam optimization algorithm and L1 regularization are used in the training
process.

##  4 Experimental Results

In this section, we study the forecasting performances of our model on
synthetic data and real data. The root mean squared error (RMSE) and mean
absolute error (MAE) are adopted as the evaluation metrics. Given a pair of
predicted and ground true matrices, the RMSE takes square root after the sum
of the squared errors of the corresponding cells, while the MAE is the sum of
the absolute errors.

###  4.1 Experiments on Synthetic Data

####  4.1.1 Simulations

To explore the superiority of our model, we first conduct simulation studies
in this section. Firstly, we generate a dataset based on the Wishart
distribution in (12) and DCAW model in (13). Secondly, we change the
underlying Wishart distribution to matrix-F distribtion [13] and generate
another dataset. That is,

$\bm{\Sigma}_{f}(t)|\mathcal{F}_{t-1}\sim
F(\bm{\nu}_{1},\bm{\nu}_{2},\frac{\bm{\nu}_{2}-n-1}{\bm{\nu}_{1}}\bm{S}_{f}(t)),$
(19)

with $E(\bm{\Sigma}_{f}(t)|\mathcal{F}_{t-1}=\bm{S}_{f}(t))$ and
$\bm{S}_{f}(t)$ is defined in (13). The conditional distribution of
$\hat{\bm{\Sigma}}_{f}(t)$ is matrix-F with a BEKK mean structure. Following
Konno,  [4] and Leung and Lo,  [6], $\bm{\Sigma}_{f}(t)$ in the above equation
could be written as

$\begin{split}\bm{\Sigma}_{f}(t)=&\left(\frac{\bm{\nu}_{2}-n-1}{\bm{\nu}_{1}}\right)\times\bm{S}_{f}(t)^{1/2}L(t)^{1/2}R(t)^{-1}\\\
&\times L(t)^{1/2}\bm{S}_{f}(t)^{1/2}\end{split}$ (20)

where $L(t)\sim Wishart(\bm{\nu}_{1},\bm{I}_{n})$, $R(t)\sim
Wishart(\bm{\nu}_{2},\bm{I}_{n})$ are independent and $\bm{I}_{n}$ is the
$n\times n$ identity matrix.

In both sets of simulations, the loading matrix $\bm{A}$ is formed by
eigenvector corresponding to the largest eigenvalues of $\bar{\bm{S}}_{y}$ in
(8) obtained from the S&P100 dataset (see Section 4.2.1). We set $r=3$ and
$(p,q)$ as $(2,2)$, while the matrices of coefficients are set as

$\bm{A}_{1}=\begin{bmatrix}0.2&0&0\\\ 0&0.4&0\\\
0&0&0.5\end{bmatrix},\bm{A}_{2}=\begin{bmatrix}0.3&0&0\\\ 0&0.5&0\\\
0&0&0.2\end{bmatrix}$ $\bm{B}_{1}=\begin{bmatrix}0.2&0&0\\\ 0&0.5&0\\\
0&0&0.4\end{bmatrix},\bm{B}_{2}=\begin{bmatrix}0.3&0&0\\\ 0&0.5&0\\\
0&0&0.2\end{bmatrix}$ $\bm{C}=\begin{bmatrix}0.5&0.2&0.3\\\ 0.2&0.5&0.25\\\
0.3&0.25&0.5\end{bmatrix}$

The initial values of $\bm{S}_{f}(t)$ are generated by the first two days’
data in S&P100 dataset. In the first simulation study, we set $\bm{\nu}=5$. In
the second simulation study, we set $\bm{\nu}_{1}=10,\bm{\nu}_{2}=8$. Each
simulation has $L=5000$ samples. The simulation from Wishart distribution can
be discribed in the following steps:

  1. 1.

Assign values to parameters,

  2. 2.

For $t=1$, set the initial value to $\bm{S}_{f}(0)$ and $\bm{S}_{f}(1)$,
sample $\bm{\Sigma}_{f}(0)$ and $\bm{\Sigma}_{f}(1)$ by (12); otherwise sample
$\bm{\Sigma}_{f}(t-1)$ based on $\bm{S}_{f}(t-1)$,

  3. 3.

Update $\bm{S}_{f}(t)$ from the BEKK recursion in (13),

  4. 4.

Set $t=t+1$ and repeat steps 2 and 3.

The simulation from F distribution would change step 2 by sampling $L(t)$ and
$R(t)$ from two independent Wishart distributions first and then obtain
$\bm{\Sigma}_{f}(t)$ by (20).

####  4.1.2 Implementation Details and Results

We conduct 10 simulations for each set of distribution and split 70%, 10% and
20% of the 5000 samples as training, validation and testing sets,
respectively. The parameters of networks are optimized by Adam optimization
algorithm with weight decay of 1e-5, $\beta_{1}$ of 0.9, $\beta_{2}$ of 0.999,
the initial learning rate of 0.001 and mini-batch sizes of 128. The $\lambda$
of L1 regularization is set as 0.005 and the $\delta$ of huber loss is set as
300. As all the samples are drawn from distributions, the simulated data do
not include ”outliers”. In this case, the $\delta$ of huber loss tends to be
larger (L2-norm loss). Table 1 illustrates the standard three-layer structure
of CM-ConvLSTM used in the simulations.

The evaluations in terms of RMSE of these 10 simulated data are displayed in
Table 2. Whether under Wishart distribution samples or F distribution samples,
our proposed CM-ConvLSTM model consistently outperforms DCAW model, by about
0.85 and 0.66 in terms of average RMSE over these two distribution samples,
respectively. These observations support our expectation that our model has
less dependence on the underlying distribution and is more robust. In the next
section, the superiority of our model will be further demonstrated on three
real-world datasets. In fact, the real data is too complex to be modeled by a
simple Wishart distribution. Our model, however, makes no distributional
assumption.

Table 1: ConvLSTM architecture for simulations with lag length 20.
Layer | Operation | Size-in | Size-out | Kernel size | Number of kernels
---|---|---|---|---|---
1 | ConvLSTM | $20\times 1\times 60\times 60$ | $8\times 60\times 60$ | $3\times 3$ | 8
leaky-relu | $8\times 60\times 60$ | $8\times 60\times 60$ | - | -
2 | Conv | $8\times 60\times 60$ | $16\times 60\times 60$ | $3\times 3$ | 16
leaky-relu | $16\times 60\times 60$ | $16\times 60\times 60$ | - | -
3 | Conv | $16\times 60\times 60$ | $1\times 60\times 60$ | $1\times 1$ | 1
linear | $1\times 60\times 60$ | $1\times 60\times 60$ | - | -
Table 2: Performance comparison of average RMSE for two set of simulation
studies, where W and F denote Wishart distribution and matrix-F distribtion
respectively.
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
DCAW (W) | 1.102 | 1.635 | 1.101 | 1.148 | 1.096 | 1.083 | 1.452 | 1.062 | 1.265 | 1.274
CM-ConvLSTM (W) | 0.359 | 0.360 | 0.362 | 0.395 | 0.371 | 0.335 | 0.351 | 0.374 | 0.357 | 0.411
difference (W) | 0.743 | 1.275 | 0.739 | 0.753 | 0.725 | 0.748 | 1.101 | 0.688 | 0.908 | 0.863
DCAW (F) | 1.473 | 1.350 | 1.575 | 1.534 | 1.559 | 1.770 | 1.602 | 1.887 | 1.620 | 1.595
CM-ConvLSTM (F) | 0.877 | 0.709 | 0.881 | 0.841 | 0.930 | 1.150 | 0.779 | 1.164 | 1.041 | 1.018
difference (F) | 0.596 | 0.642 | 0.694 | 0.693 | 0.629 | 0.621 | 0.822 | 0.723 | 0.579 | 0.577

###  4.2 Experiments on Real Data

We further show the outstanding forecasting performances of our model by
comparing the performances of some traditional time series models and some
well-known deep learning models on three real-world datasets. All the deep
learning models are trained on GPU Tesla V100 SXM2 with about 32GB memory.

####  4.2.1 Data Description

We consider three datasets: (i) all constituent stocks of the DJIA index, (ii)
all constituent stocks of the S&P100 index, (iii) all constituent stocks of
the S&P500 index. All three datasets are downloaded from the NYSE TAQ database
of WRDS, and daily RCOV matrices are calculated using the ARVM method [11].
The trading records are from 9:30 am to 4:00 pm each day, with observations
before 10:00 am deleted to avoid opening effects. The sampling frequency is
set to five minutes. Stocks with less than 100 daily trading records are also
deleted.

The first dataset comprises the 30 constituent stocks of the DJIA index from
2007 to the end of 2013. Only 25 stocks have full intra-day price data from
January 18, 2007 to December 31, 2013 (1752 trading days). RCOV matrices
($\bm{\Sigma}_{t},t=1,\dots,1752$) with dimension equals to 25 are obtained.
The second dataset contains intra-day data on 60 constituent stocks ($d=60$)
of the S&P100 index, which have full records during the period from September
10, 2003 to December 30, 2016($L=3345$). The third dataset contains a
selection of S&P500 constituent stocks that continuously traded over the
period from September 10, 2003 to December 30, 2016. The full dataset contains
244 stocks ($d=244$) and 3345 observations ($L=3345$).

The realized variances vary a lot, in some days, the variances (diagonal
values) of some stocks jump to another magnitude which distract the neural
network learning. Figure 3 reveals high volatilities around the 2008 subprime
mortgage crisis and during the flash crash on May 6,2010 for DJIA dataset. The
same phenomenon exists in the other two datasets.

Figure 3: Realized variances of the DJIA constituents from 01/18/2007 -
12/31/2013

The dimensions of the three datasets increase from 25 to 60, and finally 244,
aimed to demonstrate the power of our model in handling high dimensional RCOV
matrices. All datasets are split into three parts: the testing set contains
the last $252$ days’ RCOV matrices, the validation set contains the second
last $252$ days’ RCOV matrices, the remaining matrices all go to the training
set.

####  4.2.2 Models for Comparison

For the sake of model comparison, besides the proposed model, we also have
four time series models and three CNN-based deep learning models. The time
series models are moving average (MA), exponential moving average (EMA), MFA-
VAR and MFA-DCAW. For deep learning models, some well-known networks,
including CNN, ResNet and DenseNet are considered to apply to the datasets.

####  4.2.3 Implementation Details and Ablation Study

In the following, we conduct ablation experiments to gain a better
understanding of the impact of using different data preprocessing methods and
loss functions. Due to the limited time and computing resource, we only
conduct the comprehensive comparisons on the DJIA dataset instead of all three
datasets. The preprocessing methods include no preprocessing, Choleskey
decomposition, square root transformation and doing Choleskey decomposition
after square root transformation. The loss functions include L1-norm loss,
L2-norm loss, and Huber loss. We systematically tested different network
architectures of different combinations of filter size, the number of filters
and the number of layers on DJIA dataset. For S&P100 and S&P500 datasets, less
combinations are tested based on the experience of DJIA dataset.

Tables 3 and 4 show the average RMSE and MAE for each preprocessing method and
loss function on DJIA dataset, respectively. We can find that using Cholesky
decomposition after square root transformation of eigenvalues gives the best
performance. The Huber loss performs the best among all the loss functions.
Similar findings can be observed on S&P100 and S&P500 datasets.

Table 3: Evaluation of the impact of data preprocessing methods on DJIA
validation set.
Preprocessing | RMSE | MAE
---|---|---
$\bm{\Sigma}_{t}$ | 3.974 | 76.885
$\bm{\Sigma}_{t}\rightarrow\bm{L^{1}_{t}}\mathtt{(Cholesky~{}Decomposition)}$ | 3.774 | 72.613
$\bm{\Sigma}_{t}\rightarrow\bm{L^{2}_{t}}\mathtt{(Square~{}Root)}$ | 3.792 | 72.709
$\bm{\Sigma}_{t}\rightarrow\bm{L^{2}_{t}}\rightarrow\bm{L^{1}_{t}}$ | 3.760 | 72.139
Table 4: Evaluation of the impact of loss functions on DJIA validation set.
Loss Function | RMSE | MAE
---|---|---
L1-norm loss | 3.799 | 72.812
L2-norm loss | 3.896 | 75.437
Huber loss | 3.780 | 72.512

The parameters of networks are optimized by Adam optimization algorithm with
weight decay of 1e-5, $\beta_{1}$ of 0.9 and $\beta_{2}$ of 0.999. The initial
learning rates are 0.001, 5e-4, 5e-4 and mini-batch sizes are 128, 128, 64 for
DJIA, S&P100 and S&P500, respectively. We use L1 regularization in the
training process and we set the $\lambda$ as 0.005, 0.001, 0.0005 and the
$\delta$ of huber loss as 800, 2400, 12000 for these three datasets
respectively.

The architecture of the ConvLSTM model is very simple and is determined by the
validation RMSE and MAE. For DJIA dataset, the configuration is a standard
three-layer network, the same as in Figure 1, with the input lag length of 20,
as shown in Table 5. We find that adding more CNN encoding layers can not
further improve the model performance, then we stop at three layers for this
dataset. For S&P100 and S&P500 datasets, the networks have 1 or 2 more CNN
encoding layers respectively, and the architecture of S&P100 is explicated in
Table 6. The state-to-state kernels in the ConvLSTM encoding layer are of size
$5\times 5$ with larger receptive field, showing the spatiotemporal
correlations in the RCOVs. As time advances, the later states have the much
larger reception fields.

Table 5: ConvLSTM architecture for DJIA with lag length 20. “Conv” refers the
layer of convolution.
Layer | Operation | Size-in | Size-out | Kernel size | Number of kernels
---|---|---|---|---|---
1 | ConvLSTM | $20\times 1\times 25\times 25$ | $4\times 25\times 25$ | $3\times 3$ | 4
leaky-relu | $4\times 25\times 25$ | $4\times 25\times 25$ | - | -
2 | Conv | $4\times 25\times 25$ | $8\times 25\times 25$ | $3\times 3$ | 8
leaky-relu | $8\times 25\times 25$ | $8\times 25\times 25$ | - | -
3 | Conv | $8\times 25\times 25$ | $1\times 25\times 25$ | $1\times 1$ | 1
linear | $1\times 25\times 25$ | $1\times 25\times 25$ | - | -
Table 6: ConvLSTM architecture for S&P100 with lag length 20.
Layer | Operation | Size-in | Size-out | Kernel size | Number of kernels
---|---|---|---|---|---
1 | ConvLSTM | $20\times 1\times 60\times 60$ | $16\times 60\times 60$ | $5\times 5$ | 16
leaky-relu | $16\times 60\times 60$ | $16\times 60\times 60$ | - | -
2 | Conv | $16\times 60\times 60$ | $16\times 60\times 60$ | $3\times 3$ | 16
leaky-relu | $16\times 60\times 60$ | $16\times 60\times 60$ | - | -
3 | Conv | $16\times 60\times 60$ | $32\times 60\times 60$ | $3\times 3$ | 32
leaky-relu | $32\times 60\times 60$ | $32\times 60\times 60$ | - | -
4 | Conv | $32\times 60\times 60$ | $1\times 60\times 60$ | $5\times 5$ | 1
linear | $1\times 60\times 60$ | $1\times 60\times 60$ | - | -

Three well-known networks including CNN, ResNet and DenseNet are used to
compare with our proposed model. We first illustrate the structure of CNN
model. In our experiments, we find that the kernel size in the first layer
tends to be larger ($7\times 7$ or $9\times 9$), verifying that the existence
of local structure of RCOVs have positive effects on forecasting performances.
While the following layers are mainly non-linear mapping of the representation
of the features learnt before with $1\times 1$ or $3\times 3$ kernel size. L1
regularization is used to determine the number of kernels. For DJIA, the model
only has two layers and the numbers of kernels in the first and second layer
are 2 and 1, respectively. For S&P100 and S&P500 datasets, the architectures
are very similar, which have four and five layers, respectively. Based on the
structure of CNN, short connections are added to build the corresponding
ResNet and DenseNet model. To avoid the notable increase of the network
parameters of DenseNet for S&P100 and S&P500 datasets, we employ 1 kernel with
size $1\times 1$ to the feature maps before they are connected to each layer.

####  4.2.4 Forecasting Performances

Table 7 shows the RMSEs of all the mentioned models on DJIA, S&P100 and S&P500
testing sets. The best values of lag length are listed in the brackets next to
the RMSEs of MA and EMA models in the table, which are selected according to
the performances on validation sets. The values listed in the brackets of the
MFA-VAR models are the orders $q$ in (11) and the dimension of factor matrix
in (4). For example, (3, 1) means VAR(3) based on $r=1$ factor matrix is
chosen by the BIC criterion on training sets recommended by Shen et al.,  [7].
Similarly, for MFA-DCAW model, (1, 1, 1) means the order $(p,q)$ of DCAW in
(13) is (1, 1) based on $r=1$ factor matrix, which is chosen by the BIC
criterion on training sets. The architectures of the deep learning models for
the three datasets are selected with the smallest RMSE on validation sets, and
the number of parameters is listed next to the RMSE.

Table 7: Forecasting performances of different models on DJIA, S&P100 and
S&P500 testing sets.

Dataset | DJIA | S&P100 | S&P500
---|---|---|---
Model | Parameters | RMSE | Parameters | RMSE | Parameters | RMSE
Time Series Model | MA | (7) | 3.851 | (3) | 17.430 | (5) | 110.689
EMA | (7) | 3.734 | (5) | 16.277 | (6) | 107.586
MFA-VAR | (3, 1) | 9.542 | (3, 1) | 26.173 | (3, 1) | 136.627
MFA-DCAW | (1, 1, 1) | 7.758 | (1, 1, 1) | 25.852 | (1, 1, 1) | 134.200
Deep Learning Model | CNN | 4052 | 3.397 | 41376 | 15.054 | 61456 | 105.186
ResNet | 4078 | 3.407 | 41712 | 15.577 | 62032 | 106.593
DenseNet | 4329 | 3.401 | 42276 | 15.073 | 62526 | 106.735
ConvLSTM | 1016 | 3.394 | 34912 | 15.016 | 44128 | 104.924

Among these models, our proposed model achieved the best RMSE on all three
datasets. Our experiments show that the deep learning models perform
consistently better than the time series models, and our ConvLSTM model has
the minimum number of network parameters in these four end-to-end deep
learning models. Figure 4 compares different models by using RMSE and
correlation over time. Here, we show differnt part of time of the 252-day
testing data. According to the daily RMSE in the 252-day testing data, our
model performs the best among these days on all three datasets except for some
days with extremely large volatilities. These “outliers” come from some
positive efforts from the US Government to stimulate the economy after the
crisis. We randomly select two stocks for each dataset, which are AXP and BA
for DJIA, ABT and ACN for S&P100, AVP and HAL for S&P500, and then calculate
the correlation over time. From Figure 4, we can observe that our ConvLSTM
model catch the trend of the correlation over time.

Figure 4: Comparison of different models based on RMSE and correlation over
time on DJIA ((a) and (b)), S&P100 ((c) and (d)) and S&P500 datasets ((e) and
(f))(Best viewed in color).

For S&P100 and S&P500 datasets, the histograms of the distances of feature
maps and observed RCOV matrices of each layer are shown in Figure 5, where the
distance of each feature map and observed RCOV matrix is defined by the MAE.
Each histogram is ranged in the same scale for comparison. We find that the
distances are smaller layer by layer, and the output from the last layer shows
the smallest distance. Although the distances of the third and fourth layer of
S&P500 dataset look similar, but we find that the distances of 60% of the
covariances in the RCOV matrix are actually smaller after the third layer.

Figure 5: The histograms of distances of feature maps and ground true RCOV
matrices over four and five layers for S&P100 dataset and S&P500 dataset,
respectively.

##  5 Conclusion and Future Work

In this paper, we have successfully built an end-to-end deep learning model to
tackle the RCOV matrices forecasting problem. To our knowledge, our study is
the first attempt to apply ConvLSTM and CNN to model the RCOV matrices. The
proposed model could fit RCOV matrices well and empirically show its
outstanding performances on simulations and three real-world datasets of low
to high dimensions.

For future work, the best method of ordering the stocks with a reasonable
complexity needs further exploration. In this paper, several methods were used
to order the stocks in the RCOV matrix. The first method is to sort the stocks
by their average daily variances. The second method is to take the first PCA
component of the average RCOV matrices and use its coefficients to sort the
stocks. The third method is to use a block modeling method in social network
clustering that tries to find an ordering of stocks so that the clusters of
stocks become apparent as ’blocks’. We tested all the proposed methods of
ordering the stocks and some random orders. Based on the results, we found
that changing the ordering of stocks would affect the performances slightly
($\pm 0.4$ in MAE). This may be explained by the market globalization and most
of the stocks are weakly correlated (correlations ranged from 0.27 to 0.59)
and using filters on a local block of RCOV matrices could be sufficient in
modeling and forecasting RCOV matrices.

## References

  * Bauer and Vorkink, [2011] Bauer, G. H. and Vorkink, K. (2011).  Forecasting multivariate realized stock market volatility.  Journal of Econometrics, 160(1):93–101.
  * Deng et al., [2009] Deng, J., Dong, W., Socher, R., Li, L.-J., Li, K., and Fei-Fei, L. (2009).  Imagenet: A large-scale hierarchical image database.  In Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on, pages 248–255. IEEE.
  * Golosnoy et al., [2012] Golosnoy, V., Gribisch, B., and Liesenfeld, R. (2012).  The conditional autoregressive wishart model for multivariate stock market volatility.  Journal of Econometrics, 167(1):211–223.
  * Konno, [1991] Konno, Y. (1991).  A note on estimating eigenvalues of scale matrix of the multivariate f-distribution.  Annals of the Institute of Statistical Mathematics, 43(1):157–165.
  * Krizhevsky et al., [2012] Krizhevsky, A., Sutskever, I., and Hinton, G. E. (2012).  Imagenet classification with deep convolutional neural networks.  In Advances in Neural Information Processing Systems, pages 1097–1105.
  * Leung and Lo, [1996] Leung, P. L. and Lo, M. (1996).  An identity for the noncentral multivariate f distribution with application.  Statistica Sinica, pages 419–431.
  * Shen et al., [2018] Shen, K., Yao, J., and Li, W. K. (2018).  Forecasting high-dimensional realized volatility matrices using a factor model.  Quantitative Finance, Advance online publication(0):1–9.
  * Tao et al., [2013] Tao, M., Wang, Y., and Chen, X. (2013).  Fast convergence rates in estimating large volatility matrices using high-frequency financial data.  Econometric Theory, 29(04):838–856.
  * Tao et al., [2011] Tao, M., Wang, Y., Yao, Q., and Zou, J. (2011).  Large volatility matrix inference via combining low-frequency and high-frequency approaches.  Journal of the American Statistical Association, 106(495):1025–1040.
  * Wang, [2016] Wang, X. (2016).  Some topics in analyzing high-dimensional covariance matrices.  HKU Theses Online (HKUTO).
  * Wang et al., [2010] Wang, Y., Zou, J., et al. (2010).  Vast volatility matrix estimation for high-frequency financial data.  The Annals of Statistics, 38(2):943–978.
  * Xingjian et al., [2015] Xingjian, S., Chen, Z., Wang, H., Yeung, D.-Y., Wong, W.-K., and Woo, W.-c. (2015).  Convolutional lstm network: A machine learning approach for precipitation nowcasting.  In Advances in neural information processing systems, pages 802–810.
  * Zhou et al., [2019] Zhou, J., Jiang, F., Zhu, K., and Li, W. K. (2019).  Time series models for realized covariance matrices based on the matrix-F distribution.  Papers 1903.12077, arXiv.org.