# Next Generation Models for Portfolio Risk Management: An Approach Using
Financial Big Data

Kwangmin Junga, Donggyu Kimb and Seunghyeon Yub111corresponding author. Tel:
+82 2 958 3448.
E-mail addresses<EMAIL_ADDRESS>(K. Jung<EMAIL_ADDRESS>(D.
Kim<EMAIL_ADDRESS>(S. Yu)
aDepartment of Industrial and Management Engineering, POSTECH, Pohang, Korea
bCollege of Business, Korea Advanced Institute of Science and Technology
(KAIST), Seoul, Korea

###### Abstract

This paper proposes a dynamic process of portfolio risk measurement to address
potential information loss. The proposed model takes advantage of financial
big data to incorporate out-of-target-portfolio information that may be missed
when one considers the Value at Risk (VaR) measures only from certain assets
of the portfolio. We investigate how the curse of dimensionality can be
overcome in the use of financial big data and discuss where and when benefits
occur from a large number of assets. In this regard, the proposed approach is
the first to suggest the use of financial big data to improve the accuracy of
risk analysis. We compare the proposed model with benchmark approaches and
empirically show that the use of financial big data improves small portfolio
risk analysis. Our findings are useful for portfolio managers and financial
regulators, who may seek for an innovation to improve the accuracy of
portfolio risk estimation.

JEL classification:  C13, C32, C55, C58.

Key words and phrases:  Value at Risk, financial big data, blessing of
dimensionality, principal component analysis, multivariate GARCH, factor
model.

##  1 Introduction

The widespread failure of portfolio risk management models observed in the
2008 financial crisis raises the question of how risk managers can analyze the
risks of investment portfolios while taking into account systematic risks
caused by a systemic crisis. Extant literature has developed a range of risk
measurement methods, for example, Value at Risk (VaR), expected shortfall
(ES), entropic Value at Risk (EVaR), and superhedging price (Acerbi and
Tasche, , 2002; Ahmadi-Javid, , 2012; Bensaid et al., , 1992; Jorion, , 2000).
VaR, in particular, has been widely used to estimate the level of risk by
measuring the amount of loss for investments during a given period under
normal market conditions (Jorion, , 2000). Especially, the literature
documents the development of risk models with VaR to accommodate non-normal
market conditions, for example, asymmetric, heavy-tailed stock returns, and
volatility clustering (see, e.g., Huisman et al.,  (1998); Natarajan et al.,
(2008); Wu and Shieh,  (2007); Zoia et al.,  (2018) for non-normal stock
returns; Billio and Pelizzon,  (2000); Broda and Paolella,  (2009); Engle and
Manganelli,  (2004); Fantazzini,  (2008); McNeil and Frey,  (2000) for
volatility clustering).

However, the 2008 financial crisis demonstrates that the development of risk
models with dynamics for non-normal conditions does not guarantee financial
firms to be successful in managing extreme systematic risks. One of the
supporting arguments can be that risk managers still consider only stock data
in their target portfolio (Fan and Gu, , 2003; Hendricks, , 1996; Patton et
al., , 2019). They tend to overlook the risks of the assets that are not
included in the portfolio; hence, the portfolio size for risk analysis can be
relatively small. This so-called small portfolio risk analysis may not
sufficiently capture portfolio risk dynamics, particularly systematic risks
inherent in the market.

In this regard, this study highlights a clear message that the use of
financial big data can help risk managers take into consideration latent risks
that cannot be captured in the portfolios of interest. The contributions of
this paper are threefold. First, we propose a dynamic risk model to resolve
the above-mentioned problem by capturing systematic risk dynamics. At the core
of the proposed model is the use of financial big data to analyze portfolio
risks. The use of financial big data can help risk managers take into account
a universe of assets larger than the size of their portfolios. This usefulness
can lead to a higher accuracy of risk prediction, particularly, when it comes
to capturing the systematic risk inherent in the markets. Consequentially,
financial regulators can benefit from a high accuracy of risk prediction for
market players to maintain the stability of financial markets and thus reduce
the probability of a systemic risk occurrence.

One may be concerned about the fact that the use of financial big data can
lead multivariate VaR methods to be computationally intensive in parameter
estimation, which is called the curse of dimensionality.222For instance, when
the BEKK model (Engle and Kroner, , 1995) is conducted on $p$-dimensional log
return data, the number of parameters increases with the $p^{2}$ order. As the
dimension $p$ increases, the parameter estimation becomes computationally
demanding due to the exploding of computation time, which is the NP hard
problem. Even if it is possible to estimate the parameters, they are not
consistent estimators (Bickel and Levina, , 2008; Engle et al., , 2017). Our
second contribution is to address this point in that the proposed model
overcomes the curse of dimensionality with a two-step estimation. The first
step is to use the principal orthogonal complement thresholding (POET)
procedure (Fan et al., , 2013) to estimate the latent factor and idiosyncratic
component of the approximate latent factor model (Ait-Sahalia and Xiu, , 2017;
Fan and Kim, , 2018; Fan et al., , 2019; Kim et al., , 2018; Li et al., ,
2018). The second step is a maximum quasi-likelihood estimation for the
Generalized AutoRegressive Conditional Heteroskedasticity (GARCH) parameter
with the latent factor estimator from the first step. The estimated GARCH
parameters are used to predict one-step ahead portfolio volatility, with which
the one-step ahead VaR is predicted under the parametric and non-parametric
$\sigma$-based approaches (Brooks and Persand, , 2003; Giot and Laurent, ,
2004; Hull and White, , 1998).

Lastly, our contribution is the incorporation of financial big data featuring
the blessing of dimensionality (Donoho, , 2000; Fan and Kim, , 2019; Fan et
al., , 2013; Li et al., , 2018) by discussing the asymptotic properties of
$\sigma$-based approaches. Upon the risk measurement procedure with integrated
econometric approaches, our study is the first to suggest how a big data
approach can help financial risk management. The proposed model can be
particularly helpful for financial firms and regulators, who may be concerned
about the stability of financial markets and question how big data can be
optimally used for an improved prediction of risks.

The rest of the paper is organized as follows. In Section 2, we introduce the
concept of VaR and develop a dynamic risk model. In Section 3, we propose a
maximum quasi-likelihood estimation method and study its asymptotic behavior.
Section 4 shows how the large volatility matrix can be predicted, how the
proposed method can be applied to measure VaR, and how the curse of
dimensionality can be solved. In Section 5, we implement the Monte Carlo
simulation to check the finite sample performance of the proposed methods and
apply the proposed VaR measure to the empirical data. We conclude in Section
6.333All technical proofs are provided in the online supplement (see Appendix
A).

##  2 Model setup

Consider a market consisting of $p$ stocks. Let $\mathbf{y}_{t}$ denote the
$\mathbb{R}^{p}$-valued random vector of the assets’ log returns at time $t$.
Then, the portfolio is characterized by a vector of asset weights
$\mathbf{w}\in\mathbb{R}^{p}$, for example,
$r_{t}=\mathbf{w}^{\top}\mathbf{y}_{t}$. We assume that the log returns admit
the factor model (Ait-Sahalia and Xiu, , 2017; Bai, , 2003; Carhart, , 1997;
Fan et al., , 2013, 2019; Fama and French, , 1992, 1993, 2015) as follows:

$\displaystyle\mathbf{y}_{t}=\boldsymbol{\mu}_{t}+\mathbf{V}\mathbf{f}_{t}+\mathbf{u}_{t},$
(2.1)

where $\boldsymbol{\mu}_{t}$ is a mean vector, $\mathbf{f}_{t}$ is an
$r$-dimensional factor, $\mathbf{V}$ is a $p\times r$ factor loading matrix,
and $\mathbf{u}_{t}$ is an idiosyncratic component. We assume that the factor
$\mathbf{f}_{t}$ and idiosyncratic component $\mathbf{u}_{t}$ are independent,
and without loss of generality, we assume that the mean of $\mathbf{f}_{t}$
and $\mathbf{u}_{t}$ are zero. Usually, the number of market factors is much
smaller than the number of assets, so we assume that $r$ is finite. Moreover,
in this paper, we consider the latent factor model (Ait-Sahalia and Xiu, ,
2017; Fan et al., , 2013, 2019; Li et al., , 2018), so $\mathbf{f}_{t}$ and
$\mathbf{u}_{t}$ are not observable. In this section, under the latent factor
model, we propose a VaR estimation procedure.

###  2.1 $\sigma$-Based VaR method

One of the most popular risk measures is VaR, which estimates how much a
portfolio of assets might lose under normal market conditions. Specifically,
VaR at level $\alpha$ is defined as the $\alpha$-percentile of a portfolio log
return distribution as follows:

$\text{VaR}_{\alpha,t}=-\inf{\big{\\{}x:\mathbb{P}_{t-1}\\{r_{t}\leq
x\\}>\alpha\big{\\}}},$ (2.2)

where $\mathbb{P}_{t-1}$ denotes the conditional distribution of the portfolio
log returns $r_{t}$ given a filtration representing information to time $t-1$.
When the portfolio log return distribution follows a location-scale family,
such as normal distribution and t-distribution, we can obtain VaR from the
portfolio conditional mean $\mu_{t}$ and conditional volatility $\sigma_{t}$
given a filtration information to time $t-1$ as follows:

$\text{VaR}_{\alpha,t}=-\mu_{t}-c_{\alpha}\sigma_{t},$ (2.3)

where $c_{\alpha}$ is an $\alpha$-quantile value of the distributions of the
standardized $r_{t}$. This method is called a $\sigma$-based approach (Jorion,
, 2000). The performance of the $\sigma$-based approach depends on three
components: the conditional mean $\mu_{t}$, conditional volatility
$\sigma_{t}$, and quantile $c_{\alpha}$. Quantile $c_{\alpha}$ can be
determined using the parametric and non-parametric methods. The parametric
method involves the assumption of the parametric distribution for the
standardized log return $(r_{t}-\mu_{t})/\sigma_{t}$, such as normal
distribution or student-t distribution, and $c_{\alpha}$ is determined based
on the $\alpha$-quantile of each distribution. The non-parametric method uses
a sample quantile of the historical standardized portfolio log returns
$(r_{t}-\mu_{t})/\sigma_{t}$. Detailed implementations are described in
Section 4.2. However, several empirical studies indicate that the portfolio
conditional mean $\mu_{t}$ does not have strong time series patterns (Cajueiro
and Tabak, , 2004; Fama, , 1998; Malkiel, , 2003; Narayan and Smyth, , 2004),
which makes it difficult to predict the mean. In this paper, we simply assume
that the mean process is constant over time in order to reduce the model
misspecification error. Meanwhile, in the stock market, we observe volatility
time series structures, such as volatility clustering (Mandelbrot, , 1963).
From this point of view, it is crucial to determine the dynamic structure of
the conditional volatility $\sigma_{t}$. Thus, in this paper, we focus on how
the dynamic volatility structure can be incorporated in VaR. We note that the
mean structure would have a significant effect on VaR; hence, a parametric
structure of the mean process should be developed. We leave this for future
study.

###  2.2 Small portfolio risk analysis

When analyzing portfolio volatility, we usually investigate only the stock log
return data that belong to the portfolio. For example, let
$\mathbf{y}_{t}^{s}$ be the $s$-dimensional log return vector whose stocks
belong to the portfolio. We also denote the conditional co-volatility matrix
of $\mathbf{y}_{t}^{s}$ by $\boldsymbol{\Sigma}_{t}^{s}$ and the
$s$-dimensional portfolio weight vector by $\mathbf{w}^{s}$. Then, the
portfolio conditional volatility can be obtained from
$\boldsymbol{\Sigma}_{t}^{s}$ as follows:

$\displaystyle\sigma_{t}=\sqrt{\mathbf{w}^{s\top}\boldsymbol{\Sigma}_{t}^{s}\mathbf{w}^{s}}.$

To account for the market dynamics of the portfolio return, researchers have
introduced several dynamic models, such as the BEKK (Engle and Kroner, ,
1995), CCC (Bollerslev, , 1990), and portfolio univariate GARCH (Bollerslev, ,
1986). See also Engle,  (2002), Kawakatsu,  (2006), Paolella et al.,  (2021),
Paolella and Polak,  (2015), and Van der Weide,  (2002). The BEKK model with
the conditional mean vector $\boldsymbol{\mu}_{t}^{s}$ has the following
structure:

$\displaystyle\mathbf{y}_{t}^{s}$
$\displaystyle=\boldsymbol{\mu}_{t}^{s}+\left(\boldsymbol{\Sigma}_{t}^{s}\right)^{\frac{1}{2}}\boldsymbol{\varepsilon}_{t}^{s},$
$\displaystyle\boldsymbol{\Sigma}_{t}^{s}$
$\displaystyle=\mathbf{C}\mathbf{C}^{\top}+\mathbf{A}^{\top}(\mathbf{y}_{t-1}^{s}-\boldsymbol{\mu}_{t-1}^{s})(\mathbf{y}_{t-1}^{s}-\boldsymbol{\mu}_{t-1}^{s})^{\top}\mathbf{A}+\mathbf{B}^{\top}\boldsymbol{\Sigma}_{t-1}^{s}\mathbf{B},$
$\displaystyle\varepsilon_{it}^{s}$ $\displaystyle\sim i.i.d.\;F(0,1),$

where $\mathbf{C}$ is an $s\times s$ lower triangular matrix,
$\left(\boldsymbol{\Sigma}_{t}^{s}\right)^{\frac{1}{2}}$ is the square root of
the matrix $\boldsymbol{\Sigma}_{t}^{s}$, $\mathbf{A}$ and $\mathbf{B}$ are
$s\times s$ matrices, and $F(0,1)$ is some arbitrary distribution with mean
zero and variance one. As in the above BEKK model, the dynamic volatility
models have the autoregressive form of the historical squared returns.
Empirical studies shows that the volatility is heterogeneous, and these models
can explain the market dynamics (Bollerslev, , 1986, 1990; Engle, , 2002;
Engle and Kroner, , 1995; Kawakatsu, , 2006; Van der Weide, , 2002). Moreover,
the dynamic volatility model helps to improve the measurement of VaR (Engle
and Manganelli, , 2004; Giot and Laurent, , 2004; Kuester et al., , 2006).
That is, the performance of measuring VaR can be improved by identifying the
market dynamics. In financial markets, we often observe that much of the
market dynamic structure comes from the systematic risk, which is related to
the entire stock market; hence, the systematic risk information is contained
in the entire stock. From this point of view, even if we consider only the
small number of assets in the portfolio, to account for the market dynamics,
we need to investigate the whole stock market and extract the systematic risk
information from the stock returns. One of the naive extensions is to use a
multivariate volatility model, such as the BEKK and CCC models, for the whole
stock market. However, when applying these volatility models to the
incorporation of a large number of assets, we face the curse of dimensionality
(Bickel and Levina, , 2008; Pakel et al., , 2020; Engle et al., , 2017).
Specifically, let $p$ be the total number of assets, and then the number of
parameters quadratically increases with $p$. Hence, with a large $p$, the
estimated model is statistically inconsistent. Furthermore, in practice, the
parameter optimization takes exponential computation time, which is known as
the non-deterministic polynomial-time (NP) hard problem. Thus, applying the
usual multivariate volatility model to a large number of assets is not
feasible. In this paper, we discuss how the curse of dimensionality issue can
be handled while accounting for the market dynamics.

###  2.3 High-dimensional dynamic volatility model

In this section, we propose a dynamic factor model under the latent factor
structure in (2.1). The factor model implies that the risk of assets stems
from factor and idiosyncratic components, which are called systematic and
idiosyncratic risks, respectively. An idiosyncratic risk is related to the
firm-specific risk; hence, it does not affect the entire market. Moreover, it
can be mitigated by the portfolio diversification. By contrast, a systematic
risk arises from common market factor, such as interest rate, inflation, oil
price, etc. Because systematic risk affects the entire market, we cannot
mitigate this risk. From this point of view, we impose a dynamic structure on
the factor to account for the market risk. For the idiosyncratic risk, we
employ some martingale assumption, for example, the idiosyncratic volatility
is constant over time. Then, the conditional expected co-volatility matrix of
$\mathbf{y}_{t}$ given the current available information $\mathcal{F}_{t}$ is
expressed as follows:

$\displaystyle\mathbb{E}\left[(\mathbf{y}_{t+1}-\boldsymbol{\mu})(\mathbf{y}_{t+1}-\boldsymbol{\mu})^{\top}\middle|\mathcal{F}_{t}\right]=\boldsymbol{\Sigma}_{t+1}=\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}\mathbf{V}^{\top}+\boldsymbol{\Sigma}_{u}\text{
a.s.},$ (2.4)

where
$\boldsymbol{\Sigma}_{f,t+1}=\mathbb{E}\left[\mathbf{f}_{t+1}\mathbf{f}_{t+1}^{\top}\middle|\mathcal{F}_{t}\right]$,
and $\boldsymbol{\Sigma}_{u}$ represents the idiosyncratic volatility matrix.
Under the framework (2.4), the market dynamics can be explained by the factor
volatility dynamics.

The latent factor model has an identification problem in
$(\mathbf{V},\mathbf{f}_{t})$. For example, the pair
$(\mathbf{V},\mathbf{f}_{t})$ is not distinguishable from
$(\mathbf{V}\mathbf{H}^{\top},\mathbf{H}\mathbf{f}_{t})$ for any orthonormal
matrix $\mathbf{H}$. To uniquely define the latent factor model, we often
impose the following identifiability condition (Bai and Li, , 2012; Bai and
Ng, , 2013; Fan et al., , 2013; Kim and Fan, , 2019; Li et al., , 2018):

$\displaystyle\mathbf{V}^{\top}\mathbf{V}=p\mathbf{I}_{r}\quad\text{ and
}\quad\boldsymbol{\Sigma}_{f,t}\text{ is diagonal}.$ (2.5)

The identification condition (2.5) implies that the scaled factor loading
matrix $p^{-1/2}\mathbf{V}$ and elements of $p\boldsymbol{\Sigma}_{f,t}$ are
the eigenmatrix and eigenvalues of the conditional variance of the factor
part. Thus, the market dynamics are explained by the dynamic structure of the
eigenvalues, whereas its factor loading matrix is constant over time.

To account for the factor dynamics, we employ the famous GARCH model (Bauwens
et al., , 2006; Engle and Kroner, , 1995; Engle, , 2002, 2016; Van der Weide,
, 2002). The conditional expected volatility of the factor is modeled by the
following GARCH structure:

$\displaystyle\mathbf{f}_{t}=\left(\sqrt{h_{1t}(\boldsymbol{\theta})}\epsilon_{1t},\ldots,\sqrt{h_{rt}(\boldsymbol{\theta})}\epsilon_{rt}\right),$
$\displaystyle\mathbf{h}_{t}(\boldsymbol{\theta})=\boldsymbol{\omega}+\mathbf{A}\mathbf{f}_{t-1}^{2}+\mathbf{B}\mathbf{h}_{t-1}(\boldsymbol{\theta}),$
(2.6)

where $\epsilon_{it}$’s are i.i.d. random variables with mean zero and unit
variance; $\mathbf{h}_{t}(\boldsymbol{\theta})$ is an $r$-dimensional vector
of the conditional expected factor volatility $\boldsymbol{\Sigma}_{f,t}$,
that is, $\mathbf{h}_{t}=\big{(}h_{1t}(\boldsymbol{\theta}),\ldots,\allowbreak
h_{rt}(\boldsymbol{\theta})\big{)}=\operatorname{diag}(\boldsymbol{\Sigma}_{f,t})$,
where $\operatorname{diag}(\mathbf{X})$ is a vector whose elements are
diagonal entries of $\mathbf{X}$, $\mathbf{f}_{t}^{2}$ is an element-wise
squared vector of the factor log return $\mathbf{f}_{t}$, and
$\boldsymbol{\theta}=(\boldsymbol{\omega},\mathrm{vec}(\mathbf{A}),\mathrm{vec}(\mathbf{B}))$
is the GARCH parameters. The GARCH model in (2.3) indicates that the
conditional expected volatility of the factor is the autoregressive form of
the historical squared factor log returns. Thus, we expect that this model can
capture the stylized market features, such as the volatility clustering and
heavy tail. The empirical study in Section 5.2 also supports this.

The idiosyncratic component comes from the firm-specific risk, so their risk
is not strongly connected. Empirical studies show that there are some local
factors that affect few other idiosyncratic components (Ait-Sahalia and Xiu, ,
2017; Boivin and Ng, , 2006). In light of these, we allow a weak relationship
between idiosyncratic components so that the idiosyncratic volatility
$\boldsymbol{\Sigma}_{u}=\left((\Sigma_{u,ij})_{i,j=1,\ldots,p}\right)$
satisfies the following sparse condition:

$\displaystyle\max_{i\leq
p}{\sum_{j=1}^{p}{|\Sigma_{u,ij}|^{q}}(\Sigma_{u,ii}\Sigma_{u,jj})^{(1-q)/2}}$
$\displaystyle\leq s_{p},$ (2.7)

where $q\in[0,1)$, and the sparsity measure $s_{p}$ diverges slowly with the
dimension $p$, for example, $\log{p}$. We define $0^{0}=0$. This sparsity
condition is widely employed in the large convariance matrix inferences (Ait-
Sahalia and Xiu, , 2017; Bickel and Levina, , 2008; Cai and Liu, , 2011; Fan
et al., , 2013, 2019; Kim and Fan, , 2019). When the idiosyncratic volatility
satisfies the exact sparsity, that is, $q=0$, the sparsity condition indicates
that each asset has at most $s_{p}$ non-zero idiosyncratic correlations with
other assets.

Under the volatility structure (2.4) and (2.3), the VaR of portfolios in (2.3)
can be calculated as follows:

$\displaystyle\text{VaR}_{\alpha,t}=-\mathbf{w}^{\top}\boldsymbol{\mu}-c_{\alpha}\sqrt{\mathbf{w}^{\top}(\mathbf{V}\boldsymbol{\Sigma}_{f,t}\mathbf{V}^{\top}+\boldsymbol{\Sigma}_{u})\mathbf{w}}.$

To evaluate the above VaR value, we need to estimate the unobserved factor
components and the idiosyncratic volatility. However, to incorporate the
entire market asset information, we consider a large number of assets and
consequently run into the curse of dimensionality problem. In the following
section, we discuss how to overcome this issue and introduce an estimation
procedure to incorporate the financial big data for risk analysis.

##  3 Factor dynamics estimation

First, we define the notations. We denote $\|\cdot\|,\|\cdot\|_{F}$, and
$\|\cdot\|_{\max}$ by the matrix spectral norm, Frobenius norm, and max norm,
respectively. We use $O_{p}$ as a big-O in probability,
$\lambda_{k}(\mathbf{A})$ as the $k^{th}$ largest eigenvalue of the square
matrix $\mathbf{A}$, and $\mathrm{vec}(\mathbf{A})$ as the vectorization of
$\mathbf{A}$. We also denote the true parameters by
$\boldsymbol{\theta}_{0}=(\boldsymbol{\omega}_{0},\mathrm{vec}(\mathbf{A}_{0}),\mathrm{vec}(\mathbf{B}_{0}))$.

###  3.1 Latent factor estimation

To estimate the model parameter
$\boldsymbol{\theta}_{0}=(\boldsymbol{\omega}_{0},\mathrm{vec}(\mathbf{A}_{0}),\mathrm{vec}(\mathbf{B}_{0}))$,
we first need to estimate the latent factor components $\mathbf{f}_{t}$ and
$\mathbf{V}$. We recall that the conditional volatility matrix is

$\displaystyle\boldsymbol{\Sigma}_{t+1}=\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}\mathbf{V}^{\top}+\boldsymbol{\Sigma}_{u},$

and the mean conditional volatility matrix is

$\mkern 1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu=\mathbf{V}\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu_{f}\mathbf{V}^{T}+\boldsymbol{\Sigma}_{u},$

where $\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu_{f}=T^{-1}\sum_{t=1}^{T}\boldsymbol{\Sigma}_{f,t}$. Under the
identification and sparsity conditions (2.4) and (2.7), the mean conditional
volatility matrix has the low-rank plus sparse structure that is widely used
in the high-dimensional factor analysis (Ait-Sahalia and Xiu, , 2017; Fan et
al., , 2013, 2019; Kim and Fan, , 2019). With the low-rank plus sparse
structure, Fan et al.,  (2013) introduced the POET estimation procedure to
estimate the latent factor volatility and idiosyncratic volatility matrices.
To harness the POET procedure, we need a good proxy of the mean conditional
volatility matrix. We use the sample covariance matrix $\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu}=T^{-1}\sum_{t=1}^{T}(\mathbf{y}_{t}-\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)(\mathbf{y}_{t}-\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern 1.5mu)^{\top}$ with
sample mean $\mkern 1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu=T^{-1}\sum_{t=1}^{T}{\mathbf{y}_{t}}$, and under some mild condition,
the martingale convergence theorem implies that the sample covariance matrix
converges to $\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern 1.5mu$.
Thus, we apply the POET procedure with the sample covariance matrix to
estimate the latent factor components. Specifically, the eigenvalue
decomposition admits

$\displaystyle\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu}=\sum_{i=1}^{p}{\widehat{\lambda}_{i}\widehat{\mathbf{q}}_{i}{\widehat{\mathbf{q}}_{i}}^{\top}},$
(3.1)

where $\widehat{\lambda}_{k}$ and $\widehat{\mathbf{q}}_{k}$ are the $k^{th}$
largest eigenvalues and eigenvectors of the sample covariance matrix
$\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern 1.5mu}$,
respectively. Then, the factor loading matrix estimator $\widehat{\mathbf{V}}$
is $\sqrt{p}(\widehat{\mathbf{q}}_{1},\ldots,\widehat{\mathbf{q}}_{r})$, and
the mean factor volatility matrix estimator $\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern 1.5mu}_{f}$
is
$p^{-1}\mathrm{Diag}\big{(}(\widehat{\lambda}_{1},\ldots,\widehat{\lambda}_{r})\big{)}$,
where $\mathrm{Diag}(\mathbf{x})$ is a diagonal matrix whose diagonal entries
are $\mathbf{x}$. Note that $\widehat{\mathbf{V}}$ has a multiplier $\sqrt{p}$
so that $\widehat{\mathbf{V}}$ satisfies the condition
$\mathbf{V}^{\top}\mathbf{V}=p\mathbf{I}_{r}$. Then, the latent factors can be
obtained as follows:

$\displaystyle\widehat{\mathbf{f}}_{t}=\frac{1}{p}\widehat{\mathbf{V}}^{\top}(\mathbf{y}_{t}-\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern 1.5mu).$ (3.2)

Theorem S2.1 in the online supplement shows that the convergence rate of
$\widehat{\mathbf{f}}_{t}^{2}$ is $1/\sqrt{T}+\sqrt{s_{p}/p}$. We note that
when the factor is not observable and the number of assets is finite, the
latent factors are impossible to estimate. By contrast, Theorem S2.1 shows the
blessing of financial big data in the latent factor estimation. Specifically,
every asset contains the common market information. Thus, as the number of
assets increases, the information for the latent factors also increases.
Therefore, more stock samples exhibit a clearer latent signal, and we can
estimate the latent factors consistently with the convergence rate
$\sqrt{s_{p}/p}$.

###  3.2 Maximum quasi-likelihood estimation

In this section, we propose a model parameter estimation procedure.
Specifically, we adopt the following quasi-maximum likelihood estimation
(QMLE) method:

$\displaystyle\min_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}{\sum_{t=1}^{T}{\sum_{i=1}^{r}{\left(\log{h_{it}(\boldsymbol{\theta})}+h_{it}^{-1}(\boldsymbol{\theta})f_{it}^{2}\right)}}},$
(3.3)

where $\boldsymbol{\Theta}$ is a compact parameter space. Then, the well-
developed asymptotic theorems for the QMLE provide the consistency (Comte and
Lieberman, , 2003). However, the factors are not observable; hence, to
evaluate the quasi-likelihood function, we use the non-parametric factor
estimator $\widehat{\mathbf{f}_{t}^{2}}$. For example, the conditional co-
volatility $\mathbf{h}_{t}(\boldsymbol{\theta})$ is estimated by

$\displaystyle\widehat{\mathbf{h}}_{t}(\boldsymbol{\theta})=\boldsymbol{\omega}+\mathbf{A}\widehat{\mathbf{f}}_{t-1}^{2}+\mathbf{B}\widehat{\mathbf{h}}_{t-1}(\boldsymbol{\theta}),$
(3.4)

and we use
$\widehat{\mathbf{h}}_{1}(\boldsymbol{\theta})=(\mathbf{I}_{r}-\mathbf{A}-\mathbf{B})^{-1}\boldsymbol{\omega}$
as the initial value. Note that
$\widehat{\mathbf{h}}_{1}(\boldsymbol{\theta})$ is the unconditional
volatility of the factor. Then, we calculate the quasi-likelihood function
with the conditional co-volatility estimator
$\widehat{\mathbf{h}}_{t}(\boldsymbol{\theta})$ and obtain the maximum quasi-
likelihood estimator $\widehat{\boldsymbol{\theta}}$ as follows:

$\displaystyle\widehat{\boldsymbol{\theta}}=\arg\min_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}{\sum_{t=1}^{T}{\sum_{i=1}^{r}{\left(\log{\widehat{h}_{it}(\boldsymbol{\theta})}+\widehat{h}_{it}^{-1}(\boldsymbol{\theta})\widehat{f}_{it}^{2}\right)}}}.$
(3.5)

Theorem S3.1 in the online supplement shows the consistency of
$\widehat{\boldsymbol{\theta}}$.

##  4 Large volatility matrix estimation and VaR forecast

###  4.1 Estimation of the one-step ahead large volatility matrix

In the previous section, we find the blessing of dimensionality in the latent
factor estimation. However, when estimating large volatility matrices, we
still suffer from the curse of dimensionality. For example, sample volatility
matrix estimators are inconsistent when both the number of assets and sample
size go to infinity (Bickel and Levina, , 2008; Cai and Liu, , 2011). To
overcome the curse of dimensionality, we impose the sparse structure (2.7) on
the idiosyncratic volatility matrix. As discussed in Section 2.3, in the stock
market, the co-movement of stocks can be explained by the common factor, and
the remaining idiosyncratic co-volatilities are weakly correlated. Thus, the
sparsity condition is realistic. To estimate the sparse idiosyncratic
volatility matrix, we employ the POET procedure (Fan et al., , 2013) as
follows. First, we estimate the input idiosyncratic volatility matrix
estimator by using the non-pervasive eigen-components as follows:

$\widehat{\boldsymbol{\Sigma}}_{u}=\sum_{i=r+1}^{p}{\widehat{\lambda}_{i}\widehat{\mathbf{q}}_{i}{\widehat{\mathbf{q}}_{i}}^{\top}},$

where the eigenvalues $\widehat{\lambda}_{i}$’s and eigenvectors
$\widehat{\mathbf{q}}_{i}$’s are defined in (3.1). Then, we apply the
thresholding method to the input idiosyncratic volatility matrix estimator as
follows:

$\displaystyle[{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})]_{ij}=\begin{cases}\widehat{\Sigma}_{u,ii},&\text{
if }i=j\\\
s_{ij}(\widehat{\Sigma}_{u,ij})\mathbbm{1}_{\big{(}|\widehat{\Sigma}_{u,ij}|\geq\tau_{T}\sqrt{\widehat{\Sigma}_{u,ii}\widehat{\Sigma}_{u,jj}}\big{)}},&\text{
if }i\neq j\end{cases},$ (4.1)

where $\mathbbm{1}_{(\cdot)}$ is an indicator function, and $s_{ij}(\cdot)$ is
a shrinkage function satisfying
$|s_{ij}(x)-x|\leq\tau_{T}\sqrt{\widehat{\Sigma}_{u,ii}\widehat{\Sigma}_{u,jj}}$.
The examples of the shrink function $s_{ij}(x)$ are the soft thresholding
function
$s_{ij}(x)=x-sign(x)\tau_{T}\sqrt{\widehat{\Sigma}_{u,ii}\widehat{\Sigma}_{u,jj}}$
and the hard thresholding function $s_{ij}(x)=x$. The thresholding level
$\tau_{T}$ will be given in Theorem 4.1. The working principle of the
thresholding method is that the co-volatility is zero if the estimated
correlation is weak. This makes the estimated idiosyncratic volatility matrix
estimator sparse, so the estimated idiosyncratic volatility satisfies the
sparse condition.

With the estimated factor volatility $\widehat{\boldsymbol{\Sigma}}_{f,t+1}$
and idiosyncratic volatility matrix ${\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})$, we estimate the conditional volatility
matrix as follows:

$\displaystyle\widehat{\boldsymbol{\Sigma}}_{t+1}=\widehat{\mathbf{V}}\widehat{\boldsymbol{\Sigma}}_{f,t+1}\widehat{\mathbf{V}}^{\top}+{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u}).$ (4.2)

We call the conditional volatility matrix estimator
$\widehat{\boldsymbol{\Sigma}}_{t+1}$ the P-GARCH estimator.

The following theorem shows the asymptotic behaviors for the P-GARCH
estimator.

######  Theorem 4.1.

Under the model in Section 2.3, suppose that Assumptions 1–3 in the online
supplement hold. Take the thresholding level as
$\tau_{T}=C_{\tau}\big{(}\sqrt{\log{p}/T}+\sqrt{s_{p}/p}\big{)}$ for some
positive constant $C_{\tau}$. Then, we have

$\displaystyle\left\|{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})-\boldsymbol{\Sigma}_{u}\right\|_{\max}=O_{p}\left(\tau_{T}\right),\quad\left\|{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})-\boldsymbol{\Sigma}_{u}\right\|=O_{p}\left(s_{p}\tau_{T}^{1-q}\right),$
$\displaystyle\left\|\widehat{\boldsymbol{\Sigma}}_{t+1}-\boldsymbol{\Sigma}_{t+1}\right\|_{\max}=O_{p}\left(\tau_{T}\right),\quad\left\|\widehat{\boldsymbol{\Sigma}}_{t+1}-\boldsymbol{\Sigma}_{t+1}\right\|_{\boldsymbol{\Sigma}_{t+1}}=O_{p}\left(\frac{\sqrt{p}}{T}+s_{p}\tau_{T}^{1-q}\right),$

where the relative Frobenius norm
$\|\mathbf{G}_{1}\|_{\mathbf{G}_{2}}^{2}=p^{-1}\|\mathbf{G}_{2}^{-\frac{1}{2}}\mathbf{G}_{1}\mathbf{G}_{2}^{-\frac{1}{2}}\|_{F}^{2}$
for any given $p\times p$ matrices $\mathbf{G}_{1}$ and $\mathbf{G}_{2}$.

######  Remark 1.

Theorem 4.1 shows that the P-GARCH is the consistent estimator as long as
$p=o(T^{2})$.

With the realistic low-rank plus sparse structure, we can enjoy the blessing
of dimensionality for estimating the factor volatility matrix. Moreover, using
the regularization method, as shown in Theorem 4.1, we can overcome the curse
of dimensionality.

###  4.2 One-step ahead VaR forecast from financial big data

In this section, we discuss how the VaR value can be measured with the P-GARCH
estimator $\widehat{\boldsymbol{\Sigma}}_{t+1}$ in (4.2). Using the plug-in
method, we estimate the one-step ahead VaR as follows:

$\displaystyle\widehat{\text{VaR}}_{\alpha,t+1}=-\mathbf{w}^{\top}\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern 1.5mu-
c_{\alpha}\sqrt{\mathbf{w}^{\top}\widehat{\boldsymbol{\Sigma}}_{t+1}\mathbf{w}},$
(4.3)

where $c_{\alpha}$ is an $\alpha$-quantile, and $\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern 1.5mu$ is a sample
mean vector. To evaluate the VaR value, we need to determine the
$\alpha$-quantile value. To do this, we assume that the standardized portfolio
log returns are i.i.d. Then, when the standardized portfolio log returns
follow the standard normal distribution, $c_{\alpha}$ is the $\alpha$-quantile
of the standard normal $z_{\alpha}$. When the standardized portfolio log
returns follow the multivariate t-distribution with the degrees of freedom
$\nu$, then $c_{\alpha}=t_{\nu,\alpha}\sqrt{(\nu-2)/\nu}$, where
$t_{\nu,\alpha}$ is an $\alpha$-quantile of the t-distribution with the
degrees of freedom $\nu$ (Glasserman et al., , 2002). We call them the
parametric $\sigma$-based VaR estimator. The performance of the parametric
$\sigma$-based VaR estimator depends on the distribution assumption. On the
other hand, to obtain distribution robust VaR estimators, we use the non-
parametric sample quantile method, which we call the non-parametric
$\sigma$-based VaR estimator. For example, $c_{\alpha}$ is set to be
$\lceil\alpha T\rceil$-th smallest value of
$\\{(\mathbf{w}^{\top}(\mathbf{y}_{t}-\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)/(\mathbf{w}^{\top}\widehat{\boldsymbol{\Sigma}}_{t}\mathbf{w})^{1/2}\\}_{t=1}^{T}$,
where $\lceil\cdot\rceil$ is a ceiling function.

Then, the following theorem shows the convergence rates of VaR.

######  Theorem 4.2.

Under the model in Section 2.3, suppose that Assumptions 1–(a) in the online
supplement hold. Then, for an arbitrary portfolio weight $\mathbf{w}$ with the
gross exposure constraint $\|\mathbf{w}\|_{1}\leq C$, the parametric
$\sigma$-based VaR estimator has

$\displaystyle\left|\widehat{\text{VaR}}_{\alpha,t+1}-\text{VaR}_{\alpha,t+1}\right|=O_{p}\left(\sqrt{\frac{\log{p}}{T}}+\sqrt{\frac{s_{p}}{p}}\right).$

Moreover, the non-parametric $\sigma$-based VaR estimator is

$\displaystyle\left|\widehat{\text{VaR}}_{\alpha,t+1}-\text{VaR}_{\alpha,t+1}\right|=O_{p}\left(\sqrt{\log{T}}\left(\sqrt{\frac{\log{p}}{T}}+\sqrt{\frac{s_{p}}{p}}\right)\right).$

In this paper, we focus on investigating the effect of the VaR estimator with
the financial big data for a small portfolio. When we do not incorporate
financial big data, that is, $p$ is finite, the absolute VaR error is not
consistent for all parametric and non-parametric estimators. This is because
the term $\sqrt{s_{p}/p}$ does not converge and so the latent factor
estimation error is dominant. Therefore, Theorem 4.2 supports that
incorporating financial big data leads to the blessing in the VaR forecast by
capturing common factor dynamics.

######  Remark 2.

In this paper, we consider the parametric and non-parametric quantile
structures. However, imposing the parametric structure, such as normal and
Student $t$ distributions, on the standardized returns does not explain the
skewness of the asset returns. Furthermore, the empirical quantile approach
can have large errors in the extreme tail and does not work well when we
calculate other risk meausures, such as expected shortfall. Thus, it is
important and interesting to develop a quantile estimation procedure, which is
able to account for stylized features of the asset returns, such as skewness
and heavy-tailedness. We leave this for future study.

##  5 Numerical studies

###  5.1 Simulation study

In this section, we conducted Monte Carlo simulations to check the finite
sample performances of the proposed P-GARCH and corresponding VaR model. The
data-generating process is analogous to the models (2.1) and (2.4). The number
of factors was chosen to be $r=3$; the mean of $\mathbf{y}_{t}$ is set to be
zero $\boldsymbol{\mu}=\mathbf{0}$; and the factor loading matrix $\mathbf{V}$
was randomly sampled from the first $r$ right singular vector of random
matrix, which has the elements $\text{i.i.d.}\,\text{Unif}(0,1)$. The latent
factors $\mathbf{f}_{t}$ were generated by the multivariate normal
distribution with the conditional expected volatility
$\boldsymbol{\Sigma}_{f,t}$ as follows:

$\displaystyle\boldsymbol{\Sigma}_{f,t}=\mathrm{Diag}(\mathbf{h}_{t}(\boldsymbol{\theta}_{0}))\quad\text{
and
}\quad\mathbf{h}_{t}(\boldsymbol{\theta}_{0})=\boldsymbol{\omega}_{0}+\mathbf{A}_{0}\mathbf{f}_{t-1}^{2}+\mathbf{B}_{0}\mathbf{h}_{t-1}(\boldsymbol{\theta}_{0}),$

where
$\mathbf{h}_{1}(\boldsymbol{\theta}_{0})=\left(\mathbf{I}_{r}-\mathbf{A}_{0}-\mathbf{B}_{0}\right)^{-1}\boldsymbol{\omega}_{0}$,

$\displaystyle\boldsymbol{\omega}_{0}=\begin{pmatrix}0.003\\\ 0.002\\\
0.001\\\ \end{pmatrix},\qquad\mathbf{A}_{0}=\begin{pmatrix}0.2&0.3&0.4\\\
0.15&0.12&0.2\\\
0.1&0.1&0.1\end{pmatrix},\qquad\mathbf{B}_{0}=\begin{pmatrix}0.2&0.1&0.1\\\
0.2&0.05&0.07\\\ 0.1&0&0.05\end{pmatrix}.$

The idiosyncratic risk $\mathbf{u}_{t}$ was generated from the multivariate
normal distribution with the following sparse co-volatility:

$\displaystyle\Sigma_{u,ij}=0.01\times 0.5^{|i-j|}.$ (5.1)

We varied $p$ from 20 to 500 and $T$ from 500 to 10,000, and employed the QMLE
procedure in Section 3.2 to estimate the GARCH parameters. We repeated the
entire procedure 500 times.

Table 1: Mean absolute error (MAE) of the QMLE estimate
$\widehat{\boldsymbol{\theta}}$ with 500 replications.

$p$ | $T$ | MAE$\times 10^{2}$
---|---|---
$\omega_{1}$ | $\omega_{2}$ | $\omega_{3}$ | $A_{11}$ | $A_{12}$ | $A_{13}$ | $B_{11}$ | $B_{12}$ | $B_{13}$
20 | 500 | 0.138 | 0.088 | 0.068 | 5.706 | 9.483 | 12.637 | 14.279 | 18.768 | 21.801
2000 | 0.060 | 0.047 | 0.042 | 2.816 | 4.343 | 6.830 | 11.338 | 13.448 | 15.169
4000 | 0.046 | 0.040 | 0.034 | 2.033 | 3.359 | 5.621 | 9.738 | 10.703 | 14.345
6000 | 0.039 | 0.032 | 0.031 | 1.736 | 2.887 | 5.436 | 8.871 | 9.450 | 12.538
10000 | 0.030 | 0.031 | 0.029 | 1.333 | 2.825 | 4.947 | 8.279 | 8.912 | 12.923
100 | 500 | 0.124 | 0.079 | 0.060 | 5.816 | 8.651 | 12.264 | 13.739 | 16.993 | 19.134
2000 | 0.056 | 0.040 | 0.028 | 2.594 | 4.306 | 5.933 | 10.340 | 12.012 | 13.912
4000 | 0.040 | 0.027 | 0.020 | 1.846 | 2.925 | 3.978 | 9.869 | 10.098 | 12.730
6000 | 0.029 | 0.021 | 0.015 | 1.550 | 2.263 | 3.147 | 8.539 | 8.335 | 11.924
10000 | 0.023 | 0.018 | 0.012 | 1.047 | 1.874 | 2.606 | 7.993 | 7.731 | 10.932
500 | 500 | 0.118 | 0.079 | 0.055 | 5.788 | 8.539 | 11.931 | 14.059 | 16.716 | 18.482
2000 | 0.058 | 0.039 | 0.027 | 2.734 | 4.208 | 5.619 | 10.161 | 11.007 | 13.436
4000 | 0.039 | 0.027 | 0.019 | 1.786 | 2.916 | 3.948 | 8.975 | 9.451 | 12.525
6000 | 0.031 | 0.023 | 0.015 | 1.546 | 2.448 | 3.571 | 8.615 | 8.533 | 11.790
10000 | 0.024 | 0.018 | 0.012 | 1.057 | 1.764 | 2.456 | 7.932 | 7.477 | 10.843
|  |  |  |  |  |  |  |  |  |

Note: To save space, we report only the first 3 elements of
$\widehat{\mathbf{A}}$ and $\widehat{\mathbf{B}}$ are reported. The other
results have a similar pattern.

Table 1 shows the mean absolute error (MAE) of the QMLE estimate
$\widehat{\boldsymbol{\theta}}$ for various $p$ and $T$. For the sake of
spatial efficiency, we only documented estimation results with nine
parameters.444The estimations with the rest of parameters provide similar
results, which are available upon request. The result illustrates that MAEs
decreases as the parameter as $T$ or $p$ increases, which supports the
theoretical findings in Theorem S3.1.

Figure 1: Average predicted volatility errors of
$\widehat{\boldsymbol{\Sigma}}_{t+1}-\boldsymbol{\Sigma}_{t+1}$ against $T$
under the Frobenius, spectral, max, and relative Frobenius norms with
$p=20,100,500$.

With the estimated parameter $\widehat{\boldsymbol{\theta}}$, we validated the
one-step ahead volatility prediction $\widehat{\boldsymbol{\Sigma}}_{t+1}$
with several matrix norms. We implemented the thresholding procedure for the
idiosyncratic volatility in (4.1) with the tuning parameters $C_{\tau}=1$ and
$s_{p}=1$. Figure 1 displays the prediction errors of
$\widehat{\boldsymbol{\Sigma}}_{t+1}-\boldsymbol{\Sigma}_{t+1}$ with the
Frobenius, spectral, max, and relative Frobenius norms. We find that the
errors are large at $p=500$ for the Frobenius and spectral norms that depend
on the dimensionality $p$, whereas the error of the relative Frobenius norm
decreases as the dimension p increases. This result thus highlights that the
relative Frobenius norm can lead to the blessing of dimensionality, thereby
helping risk managers enjoy the use of financial big data for volatility
prediction. The max norm shows no significant difference for different $p$,
which can be explained by the fact that a large volatility matrix is more
likely to have a large max norm error, despite a decrease in the element-wise
error. In all cases, the errors decrease as the sample size $T$ increases,
which supports the theoretical findings in Theorem 4.1.

Figure 2: Average predicted volatility errors of the portfolio volatility
matrix $\widehat{\boldsymbol{\Sigma}}_{s,t+1}-\boldsymbol{\Sigma}_{s,t+1}$
against $T$ under the Frobenius, spectral, max, and relative Frobenius norms
with $p=500$ and portfolio size 5.

We compared the performance of predicting one-step ahead volatility matrices
of the assets in the portfolio to clarify the benefits of using financial big
data for small portfolio risk analysis We set the portfolio size as $s=5$ and
the total number of assets as $p=500$. To capture the market dynamics within
the portfolios of interest and compare with our proposed model (P-GARCH), we
considered the CCC model (Bollerslev, , 1990), and BEKK with diagonal-
constraint and variance targeting (Engle and Kroner, , 1995; Pedersen and
Rahbek, , 2014). Additionally, in the empirical study, we often impose the
static or slow time varying covariance assumption, and under this condition,
we adopted POET (Fan et al., , 2013) and historical sample covariance (Hist-
Vol). Figure 2 plots the average predicted volatility errors of the portfolio
volatility matrix
$\widehat{\boldsymbol{\Sigma}}_{s,t+1}-\boldsymbol{\Sigma}_{s,t+1}$ against
the sample size $T$ under the Frobenius, spectral, max, and relative Frobenius
norms. We observe that the P-GARCH model is superior to the other multivariate
volatility matrix models by showing much lower errors across different sample
sizes. This finding can be explained by the fact that the P-GARCH can address
the market dynamics with small portfolios, whereas the others cannot fully
address dynamics.

Finally, we compared the VaR forecasting performance with the small portfolio
risk analysis methods. We set $\alpha$ as 1% and fixed $p=500$ and varied the
portfolio size from 1 to 20. We also calculated the small portfolio matrix
estimators $\widehat{\boldsymbol{\Sigma}}_{s,t+1}$ based on the P-GARCH, CCC,
BEKK, POET, and Hist-Vol for the volatility forecast. We compared the
considered models with the portfolio univariate GARCH(1,1) model, which is a
widely-used approach to modeling a portfolio return (Bollerslev, , 1986). For
the $\alpha$-quantile, we considered the normal and student-t distribution
with 6 degrees of freedom ($\nu=6$) for the parametric method, whereas we used
the non-parametric method with $\lceil\alpha T\rceil$-th smallest values of
the standardized portfolio log returns, as in Section 4.2. The portfolio was
set to be equally weighted. The true VaR is
$z_{\alpha}(\mathbf{w}^{\top}\boldsymbol{\Sigma}_{t+1}\mathbf{w})^{1/2}$,
where $z_{\alpha}$ is the $\alpha$-quantile for the standard normal
distribution.

Figure 3: MAE of the 1%-level VaR forecasts for $p=500$ with 500 replications.
From the left, the plot shows the MAE of VaR with the parametric methods under
normal and t-distribution ($\nu=6$) and non-parametric method with the sample
quantile with the portfolio sizes 1, 5, and 20.

Figure 3 shows the MAEs of the 1%-level VaR forecasts for $p=500$ with 500
replications. We identify that the P-GARCH outperforms the other models for
any $\sigma$-based methods and portfolio size. When comparing the quantile
estimation methods, given that the true model is based on the normal
distribution, the normal distribution-based method shows the best performance.
However, the sample quantile-based method shows similar performance as the
sample size $T$ increases. This is because the sample quantile is a non-
parametric consistent estimation procedure; thus, with a large $T$, we can
estimate the quantile consistently. Finally, the parametric-based methods,
such as P-GARCH, CCC, BEKK, and port-GARCH, usually perform better than the
non-parametric methods. However, as the portfolio size increases, the BEKK
method performs worse. One possible explanation is that as the portfolio size
increases, the BEKK model becomes too complex to estimate model parameters.

###  5.2 Empirical study

In this section, we examined the performance of one-step ahead VaR forecasting
with the empirical data. The data comprise 18 years CRSP daily percentage log
returns of S&P 500 constituents from January 1, 2000, to December 31, 2017
(4,523 days). The log returns are calculated from the percentage return of the
last sale or closing bid/ask price including dividends. We selected firms that
have been ever a constituent of the S&P 500 during this period and filtered
out the stocks that have any missing data. This selection process led us to
include 492 firms ($p=492$) for the empirical study.

Figure 4: Scree plot of the sample covariance matrix of 492 stocks.

To employ the proposed P-GARCH model, we should determine the number of
factors ($r$). Figure 4 shows that the optimal $r$ is less than $3$, and the
rank $\widehat{r}$ estimated by (S1) in the online supplement is also 3. Thus,
we determined the possible values of $r$ as 1, 2, and 3, which is in line with
previous empirical studies (Chan et al., , 1998; Fama and French, , 1992,
1993). For the thresholding step, we used the global industry classification
standard (GICS) sector (Fan et al., , 2016). For example, we maintained
within-sector volatilities but set others to zero. The equally weighted
portfolios of sizes 5 and 20 are randomly sampled from the 492 stocks with 500
repetitions, and the single-asset portfolio contains only one stock from 492
assets. With these portfolios, we predicted VaR with a significance level of
$\alpha=10\%$, 5%, 2%, and 1% using the parametric and non-parametric methods,
as discussed in Section 4.2. For the model fitting, we employed a rolling
window scheme with the window size $T=252$ days (1 year). The parameter is
updated every 10 days, whereas the VaR forecasting is done every day with the
updated log return data. The total forecasting number $N$ is 4,270. The
estimated parameters of P-GARCH(r=3) for the entire period are
$\widehat{\boldsymbol{\theta}}=(\widehat{\boldsymbol{\omega}},\mathrm{vec}(\widehat{\mathbf{A}}),\mathrm{vec}(\widehat{\mathbf{B}}))=10^{-3}\times($16.32,
0.32, 0.42, 83.99, 0.00, 0.08, 6.65, 48.23, 0.90, 87.04, 5.88, 62.87, 892.7,
0.00, 0.00, 0.01, 944.7, 0.23, 0.01, 0.00, 934.99$)$. The estimated parameter
matrix $\widehat{\mathbf{B}}$ is almost diagonal, which indicates that the
recursive structure mainly comes from its own volatility. Similarly, the
estimated parameter matrix $\widehat{\mathbf{A}}$ is also almost diagonal,
except for $A_{13}$. However, when considering the magnitudes of the first and
third eigenvalues, $A_{13}$ does not affect the first eigenvalue
significantly. Thus, the dynamics of the three factors are mainly explained by
their own past volatilities. Figure 5 depicts one sample path of the estimated
one-step ahead 1%-level VaR under t-distribution with a portfolio size of 5.
It shows that predicted VaR of the P-GARCH model can capture the dynamics of
extreme loss.

We compared the out-of-sample VaR forecast performances with the small
portfolio risk analysis methods considered in Section 5.1. For example, we
used the P-GARCH, CCC, BEKK, portfolio GARCH, Hist-Vol, and POET for the
volatility prediction and the standard normal quantile, t-quantile with
degrees of freedom 6, and sample quantile for the $\alpha$-quantile. To test
whether the predicted VaR is correct, we used the unconditional coverage test
($LR_{uc}$) (Kupiec, , 1995), the conditional coverage test ($LR_{cc}$)
(Christoffersen, , 1998), and the dynamic quantile test (DQ${}_{\text{VaR}}$)
(Kuester et al., , 2006). Table 2 reports the averaged hit rates
($N_{\alpha}/N$) and the average $p$-values of the $LR_{uc},LR_{cc}$,
DQ${}_{\text{Hit}}$, and DQ${}_{\text{VaR}}$ for portfolio sizes of 1, 5, 20
and $\alpha=1\%$ with the normal distribution, student-t distribution, and
sample quantile. Figure 6 presents box plots of individual portfolio’s LRcc
$p$-values for portfolio sizes of 1, 5, and 20 and $\alpha=1\%$ with the
normal distribution, t-distribution, and sample quantile. From Table 2 and
Figure 6, we find that when comparing the dynamic models and static models,
the dynamic models, such as the P-GARCH, BEKK, CCC, and portfolio GARCH, show
better performance than the non-parametric models. This indicates that the
financial market is not static and that the GARCH-type models can account for
the market dynamics. When the dynamic models are compared, the P-GARCH has the
highest $p$-values for most of the VaR tests. From this empirical result, we
can conjecture that the P-GARCH accounts for the market dynamics by
incorporating financial big data, and this leads to improved performance in
VaR forecasting for the relatively small portfolio. Finally, when the quantile
methods are compared, the t-distribution shows the best performance. Moreover,
the normal distribution has the smallest $p$-values for all the portfolio
sizes. This pattern coincides with the stylized fact that the log return
follows the conditional heavy-tailed distribution (Cont, , 2001).

Figure 5: Plots of daily log returns and predicted one-step ahead VaRs of the
P-GARCH model.  Table 2: Average hit rate ($N_{\alpha}/N$) and average
$p$-values of the $1\%$-level VaR test statistics under the normal,
t-distribution $\nu=6$ and sample quantile for portfolio sizes of 1, 5, and
20.

Models | $N_{\alpha}/N$ | $\text{LR}_{uc}$ | $\text{LR}_{cc}$ |  DQ${}_{\text{Hit}}$ |  DQ${}_{\text{VaR}}$ |  | $N_{\alpha}/N$ | $\text{LR}_{uc}$ | $\text{LR}_{cc}$ |  DQ${}_{\text{Hit}}$ |  DQ${}_{\text{VaR}}$ |  | $N_{\alpha}/N$ | $\text{LR}_{uc}$ | $\text{LR}_{cc}$ |  DQ${}_{\text{Hit}}$ |  DQ${}_{\text{VaR}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  normal-distribution ($\text{size}=1$) |  |  normal-distribution ($\text{size}=5$) |  |  normal-distribution ($\text{size}=20$)
P-GARCH(r=1) | 0.016 | 0.043 | 0.040 | 0.024 | 0.021 |  | 0.018 | 0.002 | 0.003 | 0.001 | 0.001 |  | 0.019 | 0.000 | 0.000 | 0.000 | 0.000
P-GARCH(r=2) | 0.015 | 0.049 | 0.047 | 0.032 | 0.029 |  | 0.017 | 0.004 | 0.004 | 0.001 | 0.001 |  | 0.019 | 0.000 | 0.000 | 0.000 | 0.000
P-GARCH(r=3) | 0.015 | 0.045 | 0.042 | 0.028 | 0.025 |  | 0.017 | 0.002 | 0.002 | 0.000 | 0.001 |  | 0.019 | 0.000 | 0.000 | 0.000 | 0.000
CCC | 0.018 | 0.002 | 0.004 | 0.003 | 0.001 |  | 0.018 | 0.000 | 0.000 | 0.000 | 0.000 |  | 0.019 | 0.000 | 0.000 | 0.000 | 0.000
BEKK | 0.017 | 0.006 | 0.009 | 0.007 | 0.004 |  | 0.020 | 0.000 | 0.000 | 0.000 | 0.000 |  | 0.022 | 0.000 | 0.000 | 0.000 | 0.000
port-GARCH | 0.018 | 0.003 | 0.004 | 0.003 | 0.001 |  | 0.019 | 0.000 | 0.000 | 0.000 | 0.000 |  | 0.021 | 0.000 | 0.000 | 0.000 | 0.000
POET(r=1) | 0.018 | 0.008 | 0.007 | 0.001 | 0.001 |  | 0.021 | 0.000 | 0.000 | 0.000 | 0.000 |  | 0.023 | 0.000 | 0.000 | 0.000 | 0.000
POET(r=2) | 0.018 | 0.008 | 0.007 | 0.001 | 0.001 |  | 0.020 | 0.000 | 0.000 | 0.000 | 0.000 |  | 0.022 | 0.000 | 0.000 | 0.000 | 0.000
POET(r=3) | 0.018 | 0.008 | 0.007 | 0.001 | 0.001 |  | 0.020 | 0.000 | 0.000 | 0.000 | 0.000 |  | 0.022 | 0.000 | 0.000 | 0.000 | 0.000
Hist-Vol | 0.018 | 0.008 | 0.007 | 0.001 | 0.001 |  | 0.020 | 0.000 | 0.000 | 0.000 | 0.000 |  | 0.022 | 0.000 | 0.000 | 0.000 | 0.000
|  t-distribution ($\text{size}=1$) |  |  t-distribution ($\text{size}=5$) |  |  t-distribution ($\text{size}=20$)
P-GARCH(r=1) | 0.011 | 0.367 | 0.243 | 0.110 | 0.105 |  | 0.012 | 0.339 | 0.279 | 0.039 | 0.044 |  | 0.012 | 0.208 | 0.206 | 0.005 | 0.005
P-GARCH(r=2) | 0.011 | 0.405 | 0.262 | 0.121 | 0.117 |  | 0.012 | 0.377 | 0.309 | 0.060 | 0.061 |  | 0.012 | 0.190 | 0.202 | 0.014 | 0.012
P-GARCH(r=3) | 0.011 | 0.390 | 0.246 | 0.120 | 0.118 |  | 0.012 | 0.353 | 0.272 | 0.053 | 0.053 |  | 0.013 | 0.150 | 0.154 | 0.012 | 0.012
CCC | 0.013 | 0.202 | 0.195 | 0.093 | 0.060 |  | 0.013 | 0.212 | 0.208 | 0.012 | 0.016 |  | 0.013 | 0.125 | 0.153 | 0.001 | 0.001
BEKK | 0.012 | 0.264 | 0.250 | 0.090 | 0.078 |  | 0.014 | 0.065 | 0.066 | 0.000 | 0.001 |  | 0.016 | 0.004 | 0.004 | 0.000 | 0.000
port-GARCH | 0.013 | 0.202 | 0.195 | 0.093 | 0.061 |  | 0.013 | 0.114 | 0.128 | 0.012 | 0.008 |  | 0.014 | 0.048 | 0.072 | 0.001 | 0.001
POET(r=1) | 0.013 | 0.201 | 0.107 | 0.027 | 0.021 |  | 0.015 | 0.031 | 0.023 | 0.000 | 0.000 |  | 0.016 | 0.003 | 0.002 | 0.000 | 0.000
POET(r=2) | 0.013 | 0.201 | 0.107 | 0.027 | 0.021 |  | 0.015 | 0.036 | 0.027 | 0.000 | 0.000 |  | 0.016 | 0.004 | 0.003 | 0.000 | 0.000
POET(r=3) | 0.013 | 0.201 | 0.107 | 0.027 | 0.021 |  | 0.015 | 0.036 | 0.028 | 0.000 | 0.000 |  | 0.016 | 0.004 | 0.003 | 0.000 | 0.000
Hist-Vol | 0.013 | 0.201 | 0.107 | 0.027 | 0.021 |  | 0.015 | 0.032 | 0.024 | 0.000 | 0.000 |  | 0.016 | 0.003 | 0.002 | 0.000 | 0.000
|  sample quantile ($\text{size}=1$) |  |  sample quantile ($\text{size}=5$) |  |  sample quantile ($\text{size}=20$)
P-GARCH(r=1) | 0.013 | 0.172 | 0.130 | 0.050 | 0.032 |  | 0.013 | 0.104 | 0.120 | 0.029 | 0.028 |  | 0.014 | 0.024 | 0.031 | 0.004 | 0.005
P-GARCH(r=2) | 0.013 | 0.183 | 0.144 | 0.066 | 0.041 |  | 0.013 | 0.097 | 0.112 | 0.033 | 0.031 |  | 0.014 | 0.021 | 0.034 | 0.009 | 0.007
P-GARCH(r=3) | 0.013 | 0.152 | 0.125 | 0.050 | 0.031 |  | 0.013 | 0.075 | 0.088 | 0.025 | 0.027 |  | 0.015 | 0.012 | 0.021 | 0.005 | 0.005
CCC | 0.014 | 0.042 | 0.059 | 0.031 | 0.013 |  | 0.014 | 0.043 | 0.058 | 0.003 | 0.004 |  | 0.015 | 0.014 | 0.023 | 0.000 | 0.000
BEKK | 0.013 | 0.107 | 0.125 | 0.047 | 0.027 |  | 0.014 | 0.054 | 0.063 | 0.001 | 0.001 |  | 0.014 | 0.017 | 0.010 | 0.000 | 0.000
port-GARCH | 0.014 | 0.042 | 0.059 | 0.031 | 0.013 |  | 0.014 | 0.024 | 0.039 | 0.005 | 0.003 |  | 0.015 | 0.007 | 0.012 | 0.000 | 0.000
POET(r=1) | 0.013 | 0.107 | 0.071 | 0.015 | 0.006 |  | 0.014 | 0.043 | 0.029 | 0.000 | 0.000 |  | 0.015 | 0.016 | 0.008 | 0.000 | 0.000
POET(r=2) | 0.013 | 0.107 | 0.071 | 0.015 | 0.006 |  | 0.014 | 0.043 | 0.029 | 0.000 | 0.000 |  | 0.015 | 0.016 | 0.008 | 0.000 | 0.000
POET(r=3) | 0.013 | 0.107 | 0.071 | 0.015 | 0.006 |  | 0.014 | 0.043 | 0.029 | 0.000 | 0.000 |  | 0.015 | 0.016 | 0.008 | 0.000 | 0.000
Hist-Vol | 0.013 | 0.107 | 0.071 | 0.015 | 0.006 |  | 0.014 | 0.043 | 0.029 | 0.000 | 0.000 |  | 0.015 | 0.016 | 0.008 | 0.000 | 0.000

Figure 6: Box plots of individual portfolio’s $p$-values for $LR_{cc}$ with
$\alpha=1\%$. Each column shows three different VaR estimation results with
the normal, t-distribution ($\nu=6$), and sample quantile. Each row shows
three different VaR estimates for the portfolio sizes 1, 5, 20.

##  6 Conclusion

In this paper, we aim to develop a dynamic process of portfolio risk
measurement with the use of financial big data. The proposed model benefits
from the use of financial big data in that it leads risk managers to take into
account latent risk factors that are not included in their portfolios. One may
be concerned about the curse of dimensionality with regard to the use of
financial big data in parameter estimation; however, the proposed approach can
overcome this problem with the POET method (Fan et al., , 2013) under the
dynamic factor model in a GARCH form. We showed that a large number of assets
aids in increasing the accuracy of factor volatility estimation. This finding
results from the fact that every asset contains common market information so
that an increase in the number of assets enables one to obtain more
information on latent risk factors. Thus, the proposed model achieves the
blessing of dimensionality in estimating latent common factors.

We tested a finite sample performance of the proposed dynamic volatility model
and the corresponding VaR model with Monte Carlo simulations and showed that
the results support our theoretical findings. We further investigated the
performance of our one-step ahead volatility matrices in the asset portfolio
and concluded that our dynamic process with P-GARCH outperforms other
competitive models for any $\sigma$-based methods and portfolio size. To
guarantee the performance of our predictive model with financial big data, we
used S&P 500 constituents with parametric and non-parametric methods. The
results illustrated that the predicted VaR of our dynamic volatility model can
capture the dynamics of extreme loss, which can be extreme cases from
systematic factors in the market. Finally, the empirical study shows that our
dynamic model with P-GARCH outperforms the other competitive dynamic models in
terms of forecasting the VaR with financial big data for a relatively small
portfolio.

The lesson from the 2008 financial crisis and similar systemic events in the
financial market of history clearly argues that risk managers in the financial
industry may have to take into account latent risks possibly accounting for
systematic factors. The portfolios that they aim to analyze may consist of
idiosyncratic factors with a limited number of risk factors (or assets);
hence, the predicted level of risk in their target portfolios may fail to
address the entire risk landscape. This problem can be also of importance for
financial regulators in that underestimation of risk capitals based on a small
portfolio risk measurement would potentially prevent the financial markets
from being stabilized, thereby leading to negative externality to the economy,
as observed in 2008 (Kaserer and Klein, , 2019)555The negative externality was
realized in 2008 as a huge bailout to, e.g., AIG as described in the
introduction. Most of the $374 billion total bailout was provided to financial
firms, e.g., AIG, Citigroup and Bank of America, and the government rationale
of this financial support was to avoid the collapse of world financial markets
potentially creating the risk of another Great Depression (Harrington, ,
2009). This negative externality to the economy clearly burdened tax-payers
for not only the size of bailout from taxes, but also the negative spillover
effect on real economy.. This study contributes to this regard by providing an
innovative approach to using financial big data in improve the accuracy of
portfolio risks.

In this paper, we only consider the factor dynamics, but several studies have
also shown that the idiosyncratic volatilities are also affected by common
market factors (Barigozzi and Hallin, , 2016; Connor et al., , 2006; Herskovic
et al., , 2016; Rangel and Engle, , 2012; Shin et al., 2021b, ). Thus, it is
interesting but difficult to develop idiosyncratic dynamic models. In
addition, we consider the $\sigma$-based approach to estimate the VaR values
with the parametric and non-parametric quantile structures. However, the
parametric structure cannot account for some stylized features, such as the
skewness of the asset returns, whereas the empirical quantile approach often
has large errors in the extreme tail. Thus, it is important and interesting to
develop a quantile estimation procedure that can accomodate the stylized
features of the asset returns. Finally, we assume that the mean process is
constant over time, but this assumption is often violated in finance practice.
Unlike the volatility process, the mean process does not have a strong linear
time series structure, which makes the modeling of a dynamic mean process
difficult. Thus, a parametric structure of the mean process, which can account
for the heterogeneity of the mean process, should be developed. We leave these
interesting problems for future study.

## Online supplement

Theorems and technical proofs are provided online at Wiley in the
Supplementary Material of this article. Readers may refer to the supplementary
material associated with this article, available at Wiley
(https://onlinelibrary.wiley.com/journal/15396975).

Supplementary Material on “Next Generation Models for Portfolio Risk
Management: An Approach Using Financial Big Data”

This appendix provides micelleneous materials, theorems, and technical proofs
of the paper entitled “Next Generation of Portfolio Risk Management: An
Approach to Using Financial Big Data”.

We first introduce some notations. For a symmetric matrix $\mathbf{A}$,
$\lambda_{k}(\mathbf{A})$ is the $k^{th}$ largest eigenvalue of $\mathbf{A}$,
$\lambda_{\max}(\mathbf{A})$ is the largest eigenvalue, and
$\lambda_{\min}(\mathbf{A})$ is the smallest eigenvalue. $\|\cdot\|_{\infty}$
is the matrix $\ell_{\infty}$ norm. Denote $\mathbf{A}\circ\mathbf{B}$ the
Hadamard (or element-wise) product and
$\mathbf{x}^{2}=\mathbf{x}\circ\mathbf{x}$ the element-wise square of a vector
$\mathbf{x}$. $\operatorname{tr}(\mathbf{A})$ and $\text{vec}(\mathbf{A})$ are
the trace and vectorization of $\mathbf{A}$, respectively. Denote
$\mathbbm{1}_{E}$ the indicator function of the event $E$ and $\mathbf{e}_{i}$
the $i^{th}$ standard basis vector. To simplify the notations, we define
$\partial$ as a partial derivatives: $\partial_{x}y(\mathbf{x})=\partial
y(\mathbf{x})/\partial{x}$;
$\partial_{x}\mathbf{y}(\mathbf{x})=\partial\mathbf{y}(\mathbf{x})/\partial
x=[\partial y_{1}/\partial x,\ldots,\partial y_{\dim(\mathbf{y})}/\partial
x]^{\top}$; $\partial_{\mathbf{x}}y(\mathbf{x})=[\partial y/\partial
x_{1},\ldots,\partial y/\partial x_{\dim(\mathbf{x})}]$;
$\partial_{\mathbf{x}}\mathbf{y}(\mathbf{x})=[\partial\mathbf{y}/\partial
x_{1},\ldots,\partial\mathbf{y}/\partial x_{\dim(\mathbf{x})}]$, where the
functions $y(\mathbf{x})$ and $\mathbf{y}(\mathbf{x})$ are differentiable. $C$
is a generic constant that may differ from each equations.

##  Appendix S1 Choice of the rank

In this paper, we assume that the rank $r$ is known. However, in practice, the
rank is unknown, so it is important to estimate the rank. Ait-Sahalia and Xiu,
(2017) introduced the rank estimation procedure as follows:

$\widehat{r}=\arg\min_{1\leq j\leq
r_{\max}}p^{-1}\widehat{\lambda}_{j}+j\times c_{1}\left\\{\sqrt{\log
p/T}+p^{-1}\log p\right\\}^{c_{2}}-1,$ (S1)

where $r_{\max}$, $c_{1}$, and $c_{2}$ are tuning parameters. Under some
technical conditions, we can show its consistency. We note that the
overestimated $\widehat{r}$ is not harmful because the factor model includes
the important factors and the impact of the misspecified factors is relatively
small (Fan et al., , 2013). In contrast, when the rank is underestimated, the
factor model losses important information.

##  Appendix S2 Theorems for the latent factor estimation

######  Assumption 1.

  1. (a)

$\mathbb{E}\left[f_{it}^{8}\right],\mathbb{E}\left[h_{it}^{4}(\boldsymbol{\theta}_{0})\right],\mathbb{E}\left[\left(u_{it}u_{jt}-\Sigma_{u,ij}\right)^{4}\right]$,
and
$\mathbb{E}\big{[}\big{(}\mathbf{q}^{\top}\mathbf{u}_{t}/\sqrt{\mathbf{q}^{\top}\boldsymbol{\Sigma}_{u}\mathbf{q}}\big{)}^{4}\big{]}$
are bounded by some constant $C$ for all $i,j$ and $\mathbf{q}$ s.t.
$\|\mathbf{q}\|=1$.

  2. (b)

The minimum eigen-gap $\delta_{r}=\min_{k\leq r}\big{|}\lambda_{k}(\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu_{f})-\lambda_{k+1}(\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu_{f})\big{|}$ satisfies $\delta_{r}\geq C$ a.s.

######  Remark 3.

Assumption 1(b) is the pervasive condition that is widely used in the low-rank
matrix inferences (Ait-Sahalia and Xiu, , 2017; Fan et al., , 2013, 2018,
2019; Kim and Fan, , 2019; Stock and Watson, , 2002). The pervasive condition
with the sparsity condition (2.7) helps to distinguish the latent factor from
the idiosyncratic volatility. Additionally, since the common market factor
affects the whole asset, its proportion to the total variation is significant.
Mathematically, this implies that the corresponding eigenvalues have the $p$
order, so the pervasive condition is not restrictive.

The following theorem provides the convergence rates of the latent factor
component estimators $\widehat{\mathbf{V}}$ and $\widehat{\mathbf{f}}_{t}$.

######  Theorem S2.1.

Under the model in Section 2.3, suppose Assumption 1. Then we have

$\displaystyle\min_{\mathbf{O}}{\big{\|}\widehat{\mathbf{V}}-\mathbf{V}\mathbf{O}\big{\|}^{2}}=O_{p}\left(\frac{p}{T}+\frac{s_{p}^{2}}{p}\right),$
(S1)
$\displaystyle\big{\|}\widehat{\mathbf{f}}_{t}^{2}-\mathbf{f}_{t}^{2}\big{\|}=O_{p}\left(\frac{1}{\sqrt{T}}+\sqrt{\frac{s_{p}}{p}}\right),$
(S2)

where $\mathbf{O}$ is a diagonal sign matrix which has value $\pm 1$.

######  Remark 4.

Eigenvectors have the sign problem; for example, $-\mathbf{v}_{i}$ and
$\mathbf{v}_{i}$ are not distinguishable, where $\mathbf{v}_{i}$ is the
$i^{th}$ column vector of $\mathbf{V}$. So we put the sign matrix $\mathbf{O}$
in (S1) to identify the eigenmatrix.

######  Remark 5.

Theorem S2.1 shows that the convergence rate of $\widehat{\mathbf{f}}_{t}^{2}$
is $1/\sqrt{T}+\sqrt{s_{p}/p}$. The first term $1/\sqrt{T}$ is the usual
optimal convergence rate for estimating the mean conditional volatility matrix
$\mkern 1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu$. The second term $\sqrt{s_{p}/p}$ is the cost to identify the latent
factor. This results match with the convergence rate derived in Fan et al.,
(2011, 2013, 2018). We note that unlike the usual i.i.d. assumption, to manage
the heterogeneous volatility structure, we impose the dynamic structure on the
covariance matrix. Theorem S2.1 shows the asymptotic properties under the
heterogeneous structure.

###  S2.1 Proof of Theorem S2.1

######  Lemma S2.1.

Under the assumptions of Theorem S2.1, we have

$\displaystyle\mathbb{E}\left\|\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu}-\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu\right\|_{F}^{4}=O\left(\frac{p^{4}}{T^{2}}\right).$

Proof of Lemma S2.1. We have

$\displaystyle\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu}-\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern 1.5mu$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}{\left((\mathbf{y}_{t}-\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)(\mathbf{y}_{t}-\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)^{\top}-\boldsymbol{\Sigma}_{t}\right)}$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}{\Big{(}\mathbf{V}\left(\mathbf{f}_{t}\mathbf{f}_{t}^{\top}-\operatorname{diag}(\mathbf{h}_{0,t})\right)\mathbf{V}^{\top}+\mathbf{V}\mathbf{f}_{t}\mathbf{u}_{t}^{\top}+\mathbf{u}_{t}\mathbf{f}_{t}^{\top}\mathbf{V}^{\top}+\mathbf{u}_{t}\mathbf{u}_{t}^{\top}-\boldsymbol{\Sigma}_{u}\Big{)}}$
$\displaystyle\quad\;-(\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu-\boldsymbol{\mu})(\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu-\boldsymbol{\mu})^{\top},$

where $\mkern 1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu=(1/T)\sum_{t=1}^{T}{\mathbf{y}_{t}}$ and
$\mathbf{h}_{0,t}=\mathbf{h}_{t}(\boldsymbol{\theta}_{0})$. Then each element
is a martingale difference. Consider
$\frac{1}{T}\sum_{t=1}^{T}\mathbf{V}\left(\mathbf{f}_{t}\mathbf{f}_{t}^{\top}-\operatorname{diag}(\mathbf{h}_{0,t})\right)\mathbf{V}^{\top}$.
We have

$\displaystyle\mathbb{E}\left[\left(\frac{1}{T}\sum_{t=1}^{T}(f_{it}^{2}-h_{0,it})\right)^{4}\right]$
$\displaystyle\leq$ $\displaystyle
C\mathbb{E}\left[\left(\frac{1}{T}\sum_{t=1}^{T}(f_{it}^{2}-h_{0,it})^{2}\right)^{2}\right]$
$\displaystyle\leq$ $\displaystyle CT^{-2},$

where $h_{0,it}$ is the $i$th element of $\mathbf{h}_{0,t}$, and the first
inequality is due to the Burkholder-Davis-Gundy inequality (Chow and Teicher,
, 2012). Thus, we have

$\displaystyle\mathbb{E}\left[\left\|\frac{1}{T}\sum_{t=1}^{T}\mathbf{V}\left(\mathbf{f}_{t}\mathbf{f}_{t}^{\top}-\operatorname{diag}(\mathbf{h}_{0,t})\right)\mathbf{V}^{\top}\right\|_{F}^{4}\right]$
$\displaystyle\leq$ $\displaystyle
Cp^{4}\mathbb{E}\left[\left\|\frac{1}{T}\sum_{t=1}^{T}\left(\mathbf{f}_{t}\mathbf{f}_{t}^{\top}-\operatorname{diag}(\mathbf{h}_{0,t})\right)\
\right\|_{F}^{4}\right]$ $\displaystyle\leq$ $\displaystyle Cp^{4}T^{-2}.$

We can show other elements similarly. $\blacksquare$

Proof of Theorem S2.1. First consider (S1). By the Davis-Kahan’s sine theorem
(Yu et al., , 2015), we have

$\displaystyle\mathbb{E}\left\|\widehat{\mathbf{V}}-\mathbf{V}\mathbf{O}\right\|^{4}$
$\displaystyle=p^{2}\mathbb{E}\left\|\frac{1}{\sqrt{p}}\widehat{\mathbf{V}}-\frac{1}{\sqrt{p}}\mathbf{V}\mathbf{O}\right\|^{4}$
$\displaystyle\leq p^{2}\mathbb{E}\left[\frac{\left\|\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu}-\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu+\boldsymbol{\Sigma}_{u}\right\|^{4}}{p^{4}\delta_{r}^{4}}\right]$
$\displaystyle\leq C\left(\frac{p^{2}}{T^{2}}+\frac{s_{p}^{4}}{p^{2}}\right),$
(S3)

where the last equality is due to Lemma S2.1 and the fact that

$\displaystyle\|\boldsymbol{\Sigma}_{u}\|\leq\|\boldsymbol{\Sigma}_{u}\|_{\infty}$
$\displaystyle\leq$ $\displaystyle\max_{j}\sum_{i=1}^{p}|\Sigma_{u,ij}|$
$\displaystyle=$
$\displaystyle\max_{j}\sum_{i=1}^{p}|\Sigma_{u,ij}|^{q}|\Sigma_{u,ij}|^{1-q}$
$\displaystyle\leq$ $\displaystyle
C\max_{j}\sum_{i=1}^{p}|\Sigma_{u,ij}|^{q}\leq Cs_{p},$

where the last inequality is due to the sparsity condition 2.7.

Consider (S2). Without of loss of generality, we assume that $\mathbf{O}$ is
the identity matrix. Algebraic manipulations show

$\displaystyle\mathbb{E}\left\|\widehat{\mathbf{f}}_{t}^{2}-\mathbf{f}_{t}^{2}\right\|^{2}$
$\displaystyle=\mathbb{E}\left\|\left(\frac{1}{p}\widehat{\mathbf{V}}^{\top}(\mathbf{y}_{t}-\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)\right)^{2}-\left(\frac{1}{p}\mathbf{V}^{\top}(\mathbf{y}_{t}-\boldsymbol{\mu}-\mathbf{u}_{t})\right)^{2}\right\|^{2}$
$\displaystyle=\frac{1}{p^{4}}\mathbb{E}\left\|(\widehat{\mathbf{V}}^{\top}\mathbf{y}_{t})^{2}-2(\widehat{\mathbf{V}}^{\top}\mathbf{y}_{t}\circ\widehat{\mathbf{V}}^{\top}\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)+(\widehat{\mathbf{V}}^{\top}\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)^{2}-(\mathbf{V}^{\top}\mathbf{y}_{t})^{2}+2p\mathbf{f}_{t}\circ\mathbf{V}^{\top}(\boldsymbol{\mu}+\mathbf{u}_{t})+(\mathbf{V}^{\top}(\boldsymbol{\mu}+\mathbf{u}_{t}))^{2}\right\|^{2}$
$\displaystyle=\frac{1}{p^{4}}\mathbb{E}\Big{\|}(\widehat{\mathbf{V}}-\mathbf{V})^{\top}\mathbf{y}_{t}\circ(\widehat{\mathbf{V}}+\mathbf{V})^{\top}\mathbf{y}_{t}+2p\mathbf{f}_{t}\circ\mathbf{V}^{\top}\mathbf{u}_{t}+(\mathbf{V}^{\top}\mathbf{u}_{t})^{2}$
$\displaystyle\qquad\qquad-2(\widehat{\mathbf{V}}^{\top}\mathbf{y}_{t}\circ\widehat{\mathbf{V}}^{\top}\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)+(\widehat{\mathbf{V}}^{\top}\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)^{2}+2p\mathbf{f}_{t}\circ\mathbf{V}^{\top}\boldsymbol{\mu}+(\mathbf{V}^{\top}\boldsymbol{\mu})^{2}+2(\mathbf{V}^{\top}\mathbf{u}_{t}\circ\mathbf{V}^{\top}\boldsymbol{\mu})\Big{\|}^{2}$
$\displaystyle\leq\frac{C}{p^{4}}\bigg{(}\mathbb{E}\left\|(\widehat{\mathbf{V}}-\mathbf{V})^{\top}\mathbf{y}_{t}\circ(\widehat{\mathbf{V}}+\mathbf{V})^{\top}\mathbf{y}_{t}\right\|^{2}+p^{2}\mathbb{E}\left\|\mathbf{f}_{t}\circ\mathbf{V}^{\top}\mathbf{u}_{t}\right\|^{2}+\mathbb{E}\left\|(\mathbf{V}^{\top}\mathbf{u}_{t})^{2}\right\|^{2}$
$\displaystyle\qquad\qquad+\mathbb{E}\Big{\|}-2(\widehat{\mathbf{V}}^{\top}\mathbf{y}_{t}\circ\widehat{\mathbf{V}}^{\top}\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)+(\widehat{\mathbf{V}}^{\top}\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)^{2}+2(\mathbf{V}^{\top}\mathbf{y}_{t}\circ\mathbf{V}^{\top}\boldsymbol{\mu})-(\mathbf{V}^{\top}\boldsymbol{\mu})^{2}\Big{\|}^{2}\bigg{)}$
$\displaystyle=\frac{C}{p^{4}}\Big{(}\text{(I)}+\text{(II)}+\text{(III)}+\text{(IV)}\Big{)}.$

For (I), we have

(I)
$\displaystyle=\mathbb{E}\left\|(\widehat{\mathbf{V}}-\mathbf{V})^{\top}\mathbf{y}_{t}\circ(\widehat{\mathbf{V}}+\mathbf{V})^{\top}\mathbf{y}_{t}\right\|^{2}$
$\displaystyle\leq\mathbb{E}\left[\left\|(\widehat{\mathbf{V}}-\mathbf{V})^{\top}\mathbf{y}_{t}\right\|^{2}\left\|(\widehat{\mathbf{V}}+\mathbf{V})^{\top}\mathbf{y}_{t}\right\|^{2}\right]$
$\displaystyle\leq\mathbb{E}\left[\left\|\widehat{\mathbf{V}}-\mathbf{V}\right\|^{2}\left\|\widehat{\mathbf{V}}+\mathbf{V}\right\|^{2}\|\mathbf{y}_{t}\|^{4}\right]$
$\displaystyle\leq
4p\,\mathbb{E}\left[\left\|\widehat{\mathbf{V}}-\mathbf{V}\right\|^{2}\|\mathbf{y}_{t}\|^{4}\right]$
$\displaystyle\leq
4p\,\sqrt{\mathbb{E}\left\|\widehat{\mathbf{V}}-\mathbf{V}\right\|^{4}}\sqrt{\mathbb{E}\|\mathbf{y}_{t}\|^{8}}$
$\displaystyle\leq C\left(\frac{p^{4}}{T}+p^{2}s_{p}^{2}\right),$

where the fourth inequality is due to the Hölder’s inequality and the last
inequality is due to (S2.1). For (II), we have

(II)
$\displaystyle\leq\mathbb{E}\|\mathbf{f}_{t}\|^{2}\mathbb{E}\left\|\mathbf{V}^{\top}\mathbf{u}_{t}\right\|^{2}$
$\displaystyle\leq C\left(ps_{p}\right),$

where the last inequality is by the fact that
$\mathbb{E}\|\mathbf{V}^{\top}\mathbf{u}_{t}\|^{2}=\operatorname{tr}(\mathbf{V}^{\top}\boldsymbol{\Sigma}_{u}\mathbf{V})=O(p\,s_{p})$.
For (III), we have

(III)
$\displaystyle=\mathbb{E}\left[\sum_{i=1}^{r}(\mathbf{v}_{i}^{\top}\mathbf{u}_{t})^{4}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{i=1}^{r}\frac{(\mathbf{v}_{i}^{\top}\mathbf{u}_{t})^{4}}{(\mathbf{v}_{i}^{\top}\boldsymbol{\Sigma}_{u}\mathbf{v}_{i})^{2}}\left(\mathbf{v}_{i}^{\top}\boldsymbol{\Sigma}_{u}\mathbf{v}_{i}\right)^{2}\right]$
$\displaystyle\leq
Cp^{2}\|\boldsymbol{\Sigma}_{u}\|^{2}\mathbb{E}\left[\sum_{i=1}^{r}\frac{(\mathbf{v}_{i}^{\top}\mathbf{u}_{t})^{4}}{(\mathbf{v}_{i}^{\top}\boldsymbol{\Sigma}_{u}\mathbf{v}_{i})^{2}}\right]$
$\displaystyle\leq Cs_{p}^{2}p^{2},$

where $\mathbf{v}_{i}$ is the $i^{th}$ column of $\mathbf{V}$. Similarly, we
can show

(IV) $\displaystyle\leq
2\mathbb{E}\Big{\|}-\widehat{\mathbf{V}}^{\top}\mathbf{y}_{t}\circ\widehat{\mathbf{V}}\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu+\mathbf{V}^{\top}\boldsymbol{\mu}\circ\mathbf{V}^{\top}\mathbf{y}_{t}\Big{\|}^{2}+\mathbb{E}\Big{\|}(\widehat{\mathbf{V}}^{\top}\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)^{2}-(\mathbf{V}^{\top}\boldsymbol{\mu})^{2}\Big{\|}^{2}$
$\displaystyle\leq C\left(\frac{p^{4}}{T}+p^{2}s_{p}^{2}\right).$

Therefore, we have

$\mathbb{E}\big{\|}\widehat{\mathbf{f}}_{t}^{2}-\mathbf{f}_{t}^{2}\big{\|}^{2}=O\left(1/T+s_{p}^{2}/p^{2}+s_{p}/p\right).$

$\blacksquare$

##  Appendix S3 Theorems for the QMLE estimation procedure

Define

$\displaystyle L_{0,T}(\boldsymbol{\theta})$
$\displaystyle=-\frac{1}{T}\sum_{t=1}^{T}{\sum_{i=1}^{r}{\left(\log{h_{it}(\boldsymbol{\theta})}+h_{it}^{-1}(\boldsymbol{\theta})h_{0,it}\right)}}=\frac{1}{T}\sum_{t=1}^{T}{l_{0,t}(\boldsymbol{\theta})},$
$\displaystyle L_{T}(\boldsymbol{\theta})$
$\displaystyle=-\frac{1}{T}\sum_{t=1}^{T}{\sum_{i=1}^{r}{\left(\log{h_{it}(\boldsymbol{\theta})}+h_{it}^{-1}(\boldsymbol{\theta})f_{it}^{2}\right)}}=\frac{1}{T}\sum_{t=1}^{T}{l_{t}(\boldsymbol{\theta})},$
$\displaystyle\widehat{L}_{T}(\boldsymbol{\theta})$
$\displaystyle=-\frac{1}{T}\sum_{t=1}^{T}{\sum_{i=1}^{r}{\left(\log{\widehat{h}_{it}(\boldsymbol{\theta})}+\widehat{h}_{it}^{-1}(\boldsymbol{\theta})\widehat{f}_{it}^{2}\right)}}=\frac{1}{T}\sum_{t=1}^{T}{\widehat{l}_{t}(\boldsymbol{\theta})},$

where $\mathbf{h}_{0,t}=\mathbf{h}_{t}(\boldsymbol{\theta}_{0})$ and
$\widehat{\mathbf{h}}_{t}(\boldsymbol{\theta})=\boldsymbol{\omega}+\mathbf{A}\widehat{\mathbf{f}}_{t-1}^{2}+\mathbf{B}\widehat{\mathbf{h}}_{t-1}(\boldsymbol{\theta})$.

######  Assumption 2.

  1. (a)

The parameter space $\boldsymbol{\Theta}$ is a compact set such that every
element $\boldsymbol{\theta}\in\boldsymbol{\Theta}$ is positive;
$\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\mathbb{E}[h_{it}^{4}(\boldsymbol{\theta})]$
is bounded for all $i=1,\ldots,r$ and $t=1,\ldots,T$; the eigenvalues of
$\mathbf{B}$ and $\mathbf{A}$ are positive and $\|\mathbf{B}\|<1$; and
$\boldsymbol{\theta}_{0}$ is the interior point.

  2. (b)

$\mathbf{f}_{t}^{2}$’s are non-degenerate random variables.

The below theorem illustrates the convergence rate of the QMLE estimator
$\widehat{\boldsymbol{\theta}}$.

######  Theorem S3.1.

Under the model in Section 2.3, suppose Assumptions 1 and 2 in the online
supplement. Then we have

$\displaystyle\big{\|}\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\big{\|}_{\max}=O_{p}\left(\frac{1}{\sqrt{T}}+\sqrt{\frac{s_{p}}{p}}\right).$
(S1)

######  Remark 6.

Theorem S3.1 shows that the QMLE estimator $\widehat{\boldsymbol{\theta}}$ has
the convergence rate $1/\sqrt{T}+\sqrt{s_{p}/p}$. The term $\sqrt{s_{p}/p}$
originates from the latent factors estimation in (3.2), which is the cost to
identify the latent factor. Thus, when the factors are observable, the
convergence rate will be $1/\sqrt{T}$.

###  S3.1 Proof of Theorem S3.1

######  Lemma S3.1.

Under the assumptions of Theorem S3.1, let $\mathbf{h}_{t}$ to be a function
of $\boldsymbol{\theta}$ and $\mathbf{h}_{1}$,

$\mathbf{h}_{t}(\boldsymbol{\theta},\mathbf{h}_{1})=\mathbf{B}^{t-1}\mathbf{h}_{1}+(\mathbf{I}_{r}-\mathbf{B})^{-1}\left(\mathbf{I}_{r}-\mathbf{B}^{t-1}\right)\boldsymbol{\omega}+\sum_{k=0}^{t-2}{\mathbf{B}^{k}\mathbf{A}\mathbf{f}_{t-1-k}^{2}}.$

Then, for some true initial value $\mathbf{h}_{0,1}$, we have

$\displaystyle\left\|\mathbf{h}_{t}(\boldsymbol{\theta},\mathbf{h}_{1})-\mathbf{h}_{t}(\boldsymbol{\theta},\mathbf{h}_{0,1})\right\|$
$\displaystyle=O_{p}\left(\|\mathbf{B}\|^{t-1}\right),$
$\displaystyle\left|L_{T}(\boldsymbol{\theta},\mathbf{h}_{1})-L_{T}(\boldsymbol{\theta},\mathbf{h}_{0,1})\right|$
$\displaystyle=O_{p}\left(\frac{1}{T}\right).$

Proof of Lemma S3.1. By the sub-multiplicativity, we have

$\displaystyle\left\|\mathbf{h}_{t}(\boldsymbol{\theta},\mathbf{h}_{1})-\mathbf{h}_{t}(\boldsymbol{\theta},\mathbf{h}_{0,1})\right\|$
$\displaystyle=\left\|\mathbf{B}^{t-1}\mathbf{h}_{1}-\mathbf{B}^{t-1}\mathbf{h}_{0,1}\right\|$
$\displaystyle\leq\|\mathbf{B}\|^{t-1}\|\mathbf{h}_{1}-\mathbf{h}_{0,1}\|$
$\displaystyle=C\|\mathbf{B}\|^{t-1}.$

This implies that

$\left|L_{T}(\boldsymbol{\theta},\mathbf{h}_{1})-L_{T}(\boldsymbol{\theta},\mathbf{h}_{0,1})\right|=O_{p}\left(\frac{1}{T}\right).$

$\blacksquare$

Lemma S3.1 shows that the effect of the initial value is negligible. Thus, for
the mathematical convenience, we assume that
$\mathbf{h}_{1}(\boldsymbol{\theta})=(\mathbf{I}_{r}-\mathbf{B})^{-1}\boldsymbol{\omega}+\sum_{k=0}^{\infty}{\mathbf{B}^{k}\mathbf{A}\mathbf{f}_{-k}^{2}}$.
Then we have
$\mathbf{h}_{t}(\boldsymbol{\theta})=(\mathbf{I}_{r}-\mathbf{B})^{-1}\boldsymbol{\omega}+\sum_{k=0}^{\infty}{\mathbf{B}^{k}\mathbf{A}\mathbf{f}_{t-1-k}^{2}}$.
To obtain the consistent $\widehat{\boldsymbol{\theta}}$, it is sufficient to
show that $\widehat{L}_{T}(\boldsymbol{\theta})\xrightarrow{p}L_{0,T}$
uniformly in $\boldsymbol{\theta}$ (Lemma S3.2) and
$L_{0,T}(\boldsymbol{\theta})$ attains a unique global maximum at
$\boldsymbol{\theta}_{0}$ (Amemiya, , 1985).

######  Lemma S3.2 (Uniform Convergence).

Under the assumptions of Theorem S3.1, we have

$\displaystyle\widehat{L}_{T}(\boldsymbol{\theta})\xrightarrow{p}L_{0,T}(\boldsymbol{\theta})\text{
uniformly in }\boldsymbol{\theta}.$

Proof of Lemma S3.2. We have

$\displaystyle|\widehat{L}_{T}(\boldsymbol{\theta})-L_{0,T}(\boldsymbol{\theta})|\leq|\widehat{L}_{T}(\boldsymbol{\theta})-L_{T}(\boldsymbol{\theta})|+|L_{T}(\boldsymbol{\theta})-L_{0,T}(\boldsymbol{\theta})|.$

For the simplicity, we omit the parameter $\boldsymbol{\theta}$. Consider
$|\widehat{L}_{T}(\boldsymbol{\theta})-L_{T}(\boldsymbol{\theta})|$. We have

$\displaystyle\left|\widehat{L}_{T}(\boldsymbol{\theta})-L_{T}(\boldsymbol{\theta})\right|$
$\displaystyle\leq\frac{1}{T}\sum_{t=1}^{T}{\left(\Big{|}\sum_{i=1}^{r}{\log{\widehat{h}_{it}}-\log{h_{it}}}\Big{|}+\big{|}\widehat{\mathbf{f}}_{t}^{2}-\mathbf{f}_{t}^{2}\big{|}^{\top}\widehat{\mathbf{h}}_{t}^{-1}+{\mathbf{f}_{t}^{2}}^{\top}\big{|}\widehat{\mathbf{h}}_{t}^{-1}-\mathbf{h}_{t}^{-1}\big{|}\right)}$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}{\Big{(}\text{(I)}+\text{(II)}+\text{(III)}\Big{)}},$

where $\widehat{\mathbf{h}}_{t}^{-1}$ is the element-wise inverse of
$\widehat{\mathbf{h}}_{t}$. Recall
$\mathbf{h}_{t}(\boldsymbol{\theta})=(\mathbf{I}_{r}-\mathbf{B})^{-1}\boldsymbol{\omega}+\sum_{k=0}^{\infty}{\mathbf{B}^{k}\mathbf{A}\mathbf{f}_{t-1-k}^{2}}$.
Then by the fact that $\log{(1+x)}\leq x$ for all $x>-1$, we have

$\displaystyle\sup_{\boldsymbol{\theta}}\text{(I)}$
$\displaystyle=\sup_{\boldsymbol{\theta}}\left|\sum_{i=1}^{r}{\log{\left(1+\frac{\Delta\widehat{h}_{it}}{h_{it}}\right)}}\right|$
$\displaystyle\leq\sup_{\boldsymbol{\theta}}\sum_{i=1}^{r}{\left|\frac{\Delta\widehat{h}_{it}}{h_{it}}\right|}$
$\displaystyle\leq\frac{1}{\min_{i\leq
r}{\omega_{\min,i}}}\boldsymbol{1}^{\top}\sup_{\boldsymbol{\theta}}\left|\sum_{k=0}^{\infty}{\mathbf{B}^{k}\mathbf{A}\big{(}\Delta\widehat{\mathbf{f}}_{t-1-k}^{2}\big{)}}\right|$
$\displaystyle=o_{p}(1),$

where
$\Delta\widehat{\mathbf{h}}_{t}=\widehat{\mathbf{h}}_{t}-\mathbf{h}_{t}$,
$\Delta\widehat{\mathbf{f}}_{t}^{2}=\widehat{\mathbf{f}}_{t}^{2}-\mathbf{f}_{t}^{2}$,
and the last equality is due to Theorem S2.1. Similarly, we can show that (II)
and (III) uniformly converges to zero. Therefore, we have

$\sup_{\boldsymbol{\theta}}|\widehat{L}_{T}(\boldsymbol{\theta})-L_{T}(\boldsymbol{\theta})|=o_{p}(1).$

Now we consider $L_{T}(\boldsymbol{\theta})-L_{0,T}(\boldsymbol{\theta})$. By
Theorem 2.1 in Newey,  (1991), uniform convergence is equivalent to the
pointwise convergence and stochastic equicontinuity, where the stochastic
equicontinuity is satisfied by the Lipschitz condition
$\big{|}\partial_{\theta}\big{(}L_{T}(\boldsymbol{\theta})-L_{0,T}(\boldsymbol{\theta})\big{)}\big{|}\leq
O_{p}(1)$. Thus, it is enough to show

  1. 1.

$L_{T}(\boldsymbol{\theta})\xrightarrow{p}L_{0,T}(\boldsymbol{\theta})$
pointwise in $\boldsymbol{\theta}$;

  2. 2.

$\mathbb{E}\big{[}\sup_{\boldsymbol{\theta}}\big{|}\partial_{\theta}\big{(}L_{T}(\boldsymbol{\theta})-L_{0,T}(\boldsymbol{\theta})\big{)}\big{|}\big{]}<\infty$
for all $\theta\in\boldsymbol{\theta}$.

Since $\mathbf{f}_{t}^{2}-\mathbf{h}_{0,t}$ is a martingale difference, we can
show

$\displaystyle\mathbb{E}\left[\left(L_{T}(\boldsymbol{\theta})-L_{0,T}(\boldsymbol{\theta})\right)^{2}\right]=o(1).$

Now, we show the Lipschitz condition. We have

$\displaystyle\mathbb{E}\big{[}\sup_{\boldsymbol{\theta}}\big{|}\partial_{\theta}\big{(}L_{T}(\boldsymbol{\theta})-L_{0,T}(\boldsymbol{\theta})\big{)}\big{|}\big{]}$
$\displaystyle=\mathbb{E}\bigg{[}\sup_{\boldsymbol{\theta}}\bigg{|}\frac{1}{T}\sum_{t=1}^{T}{\sum_{i=1}^{r}{-h_{it}^{-2}(f_{it}^{2}-h_{0,it})\partial_{\theta}h_{it}}}\bigg{|}\bigg{]}$
$\displaystyle\leq
C\frac{1}{T}\sum_{t=1}^{T}{\sum_{i=1}^{r}\mathbb{E}\left[\sup_{\boldsymbol{\theta}}(\partial_{\theta}h_{it})\left(f_{it}^{2}+h_{0,it}\right)\right]}$
$\displaystyle=\frac{C}{T}\sum_{t=1}^{T}{\sum_{i=1}^{r}\mathbb{E}\left[2\sup_{\boldsymbol{\theta}}(\partial_{\theta}h_{it})\left(h_{0,it}\right)\right]}$
$\displaystyle\leq\frac{C}{T}\sum_{t=1}^{T}{\sum_{i=1}^{r}{\mathbb{E}\left[h_{0,it}^{2}\right]^{\frac{1}{2}}\mathbb{E}\left[\sup_{\boldsymbol{\theta}}(\partial_{\theta}h_{it})^{2}\right]^{\frac{1}{2}}}}$
$\displaystyle<\infty,$

where the first inequality is due to $\partial_{\theta}h_{it}>0$ for all
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$, and the second inequality is due
to Hölder’s inequality. $\blacksquare$

######  Lemma S3.3 (Uniqueness of $\boldsymbol{\theta}$).

Under the assumptions of Theorem S3.1,

$\boldsymbol{\theta}^{*}=\arg\max_{\boldsymbol{\theta}}{L_{0,T}(\boldsymbol{\theta})}$

is unique almost surely, and
$\boldsymbol{\theta}^{*}=\boldsymbol{\theta}_{0}$. Moreover, we have

$\displaystyle\widehat{\boldsymbol{\theta}}\xrightarrow{p}\boldsymbol{\theta}_{0}.$
(S2)

Proof of Lemma S3.3. Since $\log{x}+t/x$ has a unique minimizer at $x=t$,
$\mathbf{h}_{it}(\boldsymbol{\theta}_{0})$ is a unique maximizer of
$l_{0,t}(\mathbf{h}_{it})$ for all $t$. If $\boldsymbol{\theta}_{0}$ is not a
unique parameter to have $\mathbf{h}_{it}(\boldsymbol{\theta}_{0})$, then
there exists $\boldsymbol{\theta}^{*}\neq\boldsymbol{\theta}_{0}$ such that
$\mathbf{h}_{t}(\boldsymbol{\theta}^{*})=\mathbf{h}_{t}(\boldsymbol{\theta}_{0})$.
Then
$\\{\mathbf{h}_{t}(\boldsymbol{\theta}^{*})-\mathbf{h}_{t}(\boldsymbol{\theta}_{0})\\}_{t\leq
T}=\\{\mathbf{0}\\}_{t\leq T}$, and

$\displaystyle\begin{pmatrix}\mathbf{h}_{2}(\boldsymbol{\theta}^{*})-\mathbf{h}_{2}(\boldsymbol{\theta}_{0}),\ldots,\mathbf{h}_{T}(\boldsymbol{\theta}^{*})-\mathbf{h}_{T}(\boldsymbol{\theta}_{0})\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\boldsymbol{\omega}^{*}-\boldsymbol{\omega}_{0}&\mathbf{A}^{*}-\mathbf{A}_{0}&\mathbf{B}^{*}-\mathbf{B}_{0}\end{pmatrix}\begin{pmatrix}1&1&\cdots&1\\\
\mathbf{f}_{1}^{2}&\mathbf{f}_{2}^{2}&\cdots&\mathbf{f}_{T-1}^{2}\\\
\mathbf{h}_{1}(\boldsymbol{\theta}_{0})&\mathbf{h}_{2}(\boldsymbol{\theta}_{0})&\cdots&\mathbf{h}_{T-1}(\boldsymbol{\theta}_{0})\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\boldsymbol{\omega}^{*}-\boldsymbol{\omega}_{0}&\mathbf{A}^{*}-\mathbf{A}_{0}&\mathbf{B}^{*}-\mathbf{B}_{0}\end{pmatrix}\mathbf{M}=\mbox{\large
0},$

where 0 is a zero matrix. By Assumption 2(b), $\mathbf{M}\mathbf{M}^{\top}$ is
invertible a.s., which implies

$\begin{pmatrix}\boldsymbol{\omega}^{*}-\boldsymbol{\omega}_{0}&\mathbf{A}^{*}-\mathbf{A}_{0}&\mathbf{B}^{*}-\mathbf{B}_{0}\end{pmatrix}=\mbox{\large
0}.$

Therefore, $\boldsymbol{\theta}^{*}=\boldsymbol{\theta}_{0}$ a.s. which is a
contradiction.

With the uniqueness solution result, by Lemma S3.2 and Theorem 4.1.2 in
Amemiya,  (1985), we can show (S2). $\blacksquare$

######  Lemma S3.4.

Under the assumptions of Theorem S3.1, for all fixed $c\geq 1$, $l\leq r$ and
$\theta\in\boldsymbol{\theta}$,

$\mathbb{E}\left[\left|\sup_{\boldsymbol{\theta}}{\frac{\partial_{\theta}\widehat{h}_{lt}}{\widehat{h}_{lt}}}\right|^{c}\right]<\infty,\;\mathbb{E}\left[\left|\sup_{\boldsymbol{\theta}}{\frac{\partial_{\theta}^{2}\widehat{h}_{lt}}{\widehat{h}_{lt}}}\right|^{c}\right]<\infty,\;\mathbb{E}\left[\left|\sup_{\boldsymbol{\theta}}{\frac{\partial_{\theta}^{3}\widehat{h}_{lt}}{\widehat{h}_{lt}}}\right|^{c}\right]<\infty.$

Proof of Lemma S3.4. Notice that

$\partial_{\theta}\widehat{\mathbf{h}}_{t}=\partial_{\theta}\boldsymbol{\varpi}+\sum_{k=0}^{\infty}{\mathbf{B}^{k}(\partial_{\theta}\mathbf{A})\widehat{\mathbf{f}}_{t-1-k}^{2}}+\sum_{k=1}^{\infty}{\left(\sum_{\xi=0}^{k-1}{\mathbf{B}^{\xi}(\partial_{\theta}\mathbf{B})\mathbf{B}^{k-1-\xi}}\right)\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}}\geq\mathbf{0}.$

where
$\boldsymbol{\varpi}=(\mathbf{I}_{r}-\mathbf{B})^{-1}\boldsymbol{\omega}$. For
the case of $\theta=\varpi_{i}$,
$\mathbb{E}\big{[}\big{|}\sup_{\boldsymbol{\theta}}{\partial_{\theta}\widehat{h}_{lt}/\widehat{h}_{lt}}\big{|}^{c}\big{]}<\infty$
is trivial. With the case $\theta=A_{ij}$,

$\displaystyle\frac{\partial_{\theta}\widehat{h}_{lt}}{\widehat{h}_{lt}}$
$\displaystyle\leq\frac{\mathbf{e}_{l}^{\top}\sum_{k=0}^{\infty}{\mathbf{B}^{k}(\partial_{\theta}\mathbf{A})\widehat{\mathbf{f}}_{t-1-k}^{2}}}{\mathbf{e}_{l}^{\top}\sum_{k=0}^{\infty}{\mathbf{B}^{k}\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}}}$
$\displaystyle\leq\frac{\mathbf{e}_{l}^{\top}\sum_{k=0}^{\infty}{\mathbf{B}^{k}\mathbf{e}_{i}\mathbf{e}_{j}^{\top}\widehat{\mathbf{f}}_{t-1-k}^{2}}}{\mathbf{e}_{l}^{\top}\sum_{k=0}^{\infty}{\mathbf{B}^{k}\mathbf{e}_{i}A_{ij}\mathbf{e}_{j}^{\top}\widehat{\mathbf{f}}_{t-1-k}^{2}}}$
$\displaystyle=\frac{1}{A_{ij}}\leq\frac{1}{A_{\min,ij}},$

where $\mathbf{e}_{i}$ is the $i^{th}$ standard basis vector. Thus,
$\mathbb{E}\big{[}\big{|}\sup_{\boldsymbol{\theta}}{\partial_{\theta}\widehat{h}_{lt}/\widehat{h}_{lt}}\big{|}^{c}\big{]}<\infty$.
For the case of $\theta=B_{ij}$, we have

$\displaystyle\frac{\partial_{\theta}\widehat{h}_{lt}}{\widehat{h}_{lt}}$
$\displaystyle=\frac{\mathbf{e}_{l}^{\top}\sum_{k=1}^{\infty}{\left(\sum_{\xi=0}^{k-1}{\mathbf{B}^{\xi}(\partial_{\theta}\mathbf{B})\mathbf{B}^{k-1-\xi}}\right)\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}}}{\varpi_{l}+\mathbf{e}_{l}^{\top}\sum_{k=0}^{\infty}{\mathbf{B}^{k}\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}}}$
$\displaystyle\leq\sum_{k=1}^{\infty}{\frac{\mathbf{e}_{l}^{\top}\left(\sum_{\xi=0}^{k-1}{\mathbf{B}^{\xi}\mathbf{e}_{i}\mathbf{e}^{\top}_{j}\mathbf{B}^{k-1-\xi}}\right)\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}}{\varpi_{l}+\mathbf{e}_{l}^{\top}\left(\sum_{\xi=0}^{k-1}\mathbf{B}^{\xi}\mathbf{e}_{i}B_{ij}\mathbf{e}^{\top}_{j}\mathbf{B}^{k-1-\xi}\right)\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}}}$
$\displaystyle=\sum_{k=0}^{\infty}{\frac{\mathbf{e}_{l}^{\top}\left(\sum_{\xi=0}^{k-1}{\mathbf{B}^{\xi}\mathbf{e}_{i}\mathbf{e}^{\top}_{j}\mathbf{B}^{k-1-\xi}}\right)\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}}{\varpi_{l}+B_{ij}\mathbf{e}_{l}^{\top}\left(\sum_{\xi=0}^{k-1}{\mathbf{B}^{\xi}\mathbf{e}_{i}\mathbf{e}^{\top}_{j}\mathbf{B}^{k-1-\xi}}\right)\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}}}$
$\displaystyle=\frac{1}{B_{ij}}\sum_{k=0}^{\infty}{\frac{(B_{ij}/\varpi_{l})\mathbf{e}_{l}^{\top}\left(\sum_{\xi=0}^{k-1}{\mathbf{B}^{\xi}\mathbf{e}_{i}\mathbf{e}^{\top}_{j}\mathbf{B}^{k-1-\xi}}\right)\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}}{1+(B_{ij}/\varpi_{l})\mathbf{e}_{l}^{\top}\left(\sum_{\xi=0}^{k-1}{\mathbf{B}^{\xi}\mathbf{e}_{i}\mathbf{e}^{\top}_{j}\mathbf{B}^{k-1-\xi}}\right)\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}}}$
$\displaystyle\leq\frac{1}{B_{ij}}\frac{B_{ij}^{\alpha}}{\varpi_{l}^{\alpha}}\sum_{k=0}^{\infty}{\left[\mathbf{e}_{l}^{\top}\left(\sum_{\xi=0}^{k-1}{\mathbf{B}^{\xi}\mathbf{e}_{i}\mathbf{e}^{\top}_{j}\mathbf{B}^{k-1-\xi}}\right)\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}\right]^{\alpha}}$
$\displaystyle\leq\frac{1}{B_{ij}}\frac{1}{\varpi_{l}^{\alpha}}\sum_{k=0}^{\infty}{\left[\mathbf{e}_{l}^{\top}\mathbf{B}^{k}\mathbf{A}\widehat{\mathbf{f}}_{t-1-k}^{2}\right]^{\alpha}},$

where the second inequality is due to $x/(1+x)\leq x^{\alpha}$ for
$\forall{\alpha}\in[0,1]$. By choosing $\alpha=1/c$, we have

$\displaystyle\left(\mathbb{E}\left[\sup_{\boldsymbol{\theta}}\left|\frac{\partial_{\theta}\widehat{h}_{lt}}{\widehat{h}_{lt}}\right|^{c}\right]\right)^{\frac{1}{c}}$
$\displaystyle\leq
C\sum_{k=0}^{\infty}\|\mathbf{B}_{\max}\|^{k}\mathbb{E}\left[\|\widehat{\mathbf{f}}_{t-1-k}^{2}\|\right]$
$\displaystyle<\infty,$

where first inequality is due to the Minkowski’s inequality, and the last
inequality follows from $\mathbb{E}\big{[}|f_{it}|^{8}\big{]}\leq C$ and
$\mathbb{E}\big{\|}\widehat{\mathbf{f}}_{t}^{2}-\mathbf{f}_{t}^{2}\big{\|}^{2}=o(1)$.
Similarly, we can show the higher order derivatives bounds. $\blacksquare$

######  Lemma S3.5.

Under the assumptions of Theorem S3.1, let
$\mathbb{E}\left\|\Delta\mathbf{f}_{t}^{2}\right\|^{2}=O(\beta_{T}^{2})$. Then
we have

$\left|\partial_{\theta}\widehat{L}_{T}(\boldsymbol{\theta}_{0})-\partial_{\theta}L_{0,T}(\boldsymbol{\theta}_{0})\right|=O_{p}\left(\beta_{T}+\frac{1}{\sqrt{T}}\right).$

Proof of Lemma S3.5. Denote
$\widehat{\mathbf{h}}_{0,t}=\widehat{\mathbf{h}}_{t}(\boldsymbol{\theta}_{0})$
and
$\partial_{\theta}\widehat{\mathbf{h}}_{0,t}=\partial_{\theta}\widehat{\mathbf{h}}_{t}(\boldsymbol{\theta}_{0})$.
Since
$\mathbb{E}\left\|\Delta\mathbf{f}_{t}^{2}\right\|^{2}=O(\beta_{T}^{2})$, we
have

  1. 1.

$\mathbb{E}\left\|\widehat{\mathbf{h}}_{0,t}-\mathbf{h}_{0,t}\right\|^{2}=O(\beta_{T}^{2})$;

  2. 2.

$\mathbb{E}\left\|\widehat{\mathbf{h}}_{0,t}^{-1}-\mathbf{h}_{0,t}^{-1}\right\|^{2}=O(\beta_{T}^{2})$;

  3. 3.

$\mathbb{E}\left\|\partial_{\theta}\widehat{\mathbf{h}}_{0,t}-\partial_{\theta}\mathbf{h}_{0,t}\right\|^{2}=O(\beta_{T}^{2})$;

  4. 4.

$\mathbb{E}\left[{\mathbf{f}_{t}^{2}}^{\top}(\widehat{\mathbf{h}}_{0,t}^{-2}-\mathbf{h}_{0,t}^{-2})\circ\partial_{\theta}\widehat{\mathbf{h}}_{0,t}\right]=O(\beta_{T})$,

where the last equation is by the fact that

$\displaystyle\mathbb{E}\left[{\mathbf{f}_{t}^{2}}^{\top}(\widehat{\mathbf{h}}_{0,t}^{-2}-\mathbf{h}_{0,t}^{-2})\circ\partial_{\theta}\widehat{\mathbf{h}}_{0,t}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{l=1}^{r}{\left(\frac{h_{0,lt}^{2}-\widehat{h}_{0,lt}^{2}}{h_{0,lt}^{2}\widehat{h}_{0,lt}^{2}}\right)f_{lt}^{2}\partial_{\theta}\widehat{h}_{0,lt}}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{l=1}^{r}{\left(\frac{h_{0,lt}^{2}-\widehat{h}_{0,lt}^{2}}{h_{0,lt}\widehat{h}_{0,lt}}\right)\epsilon_{lt}^{2}\frac{\partial_{\theta}\widehat{h}_{0,lt}}{\widehat{h}_{0,lt}}}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{l=1}^{r}{\left(\frac{h_{0,lt}-\widehat{h}_{0,lt}}{\widehat{h}_{0,lt}}\right)\epsilon_{lt}^{2}\frac{\partial_{\theta}\widehat{h}_{0,lt}}{\widehat{h}_{0,lt}}+\left(\frac{h_{0,lt}-\widehat{h}_{0,lt}}{h_{0,lt}}\right)\epsilon_{lt}^{2}\frac{\partial_{\theta}\widehat{h}_{0,lt}}{\widehat{h}_{0,lt}}}\right]$
$\displaystyle\leq
C\mathbb{E}\left[\sum_{l=1}^{r}{\left|h_{0,lt}-\widehat{h}_{0,lt}\right|\epsilon_{lt}^{2}\frac{\partial_{\theta}\widehat{h}_{0,lt}}{\widehat{h}_{0,lt}}}\right]$
$\displaystyle\leq
C\sum_{l=1}^{r}{\left(\mathbb{E}\left[\left(h_{0,lt}-\widehat{h}_{0,lt}\right)^{2}\right]\right)^{\frac{1}{2}}\left(\mathbb{E}\left[\left(\frac{\partial_{\theta}\widehat{h}_{0,lt}}{\widehat{h}_{0,lt}}\right)^{2}\right]\right)^{\frac{1}{2}}}$
$\displaystyle=O(\beta_{T}),$

where the last equality is due to Lemma S3.4. Thus, we have

$\displaystyle\partial_{\theta}\widehat{L}_{T}(\boldsymbol{\theta}_{0})-\partial_{\theta}L_{T}(\boldsymbol{\theta}_{0})$
$\displaystyle=-\frac{1}{T}\sum_{t=1}^{T}{(\widehat{\mathbf{h}}_{0,t}-\widehat{\mathbf{f}}_{t}^{2})^{\top}(\widehat{\mathbf{h}}_{0,t}^{-2}\circ\partial_{\theta}\widehat{\mathbf{h}}_{0,t})-(\mathbf{h}_{0,t}-\mathbf{f}_{t}^{2})^{\top}(\mathbf{h}_{0,t}^{-2}\circ\partial_{\theta}\mathbf{h}_{0,t})}$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\bigg{\\{}(\widehat{\mathbf{f}}_{t}^{2}-\mathbf{f}_{t}^{2})^{\top}(\widehat{\mathbf{h}}_{0,t}^{-2}\circ\partial_{\theta}\widehat{\mathbf{h}}_{0,t})+{\mathbf{f}_{t}^{2}}^{\top}(\widehat{\mathbf{h}}_{0,t}^{-2}-\mathbf{h}_{0,t}^{-2})\circ\partial_{\theta}\widehat{\mathbf{h}}_{0,t}$
$\displaystyle\qquad\qquad+{\mathbf{f}_{t}^{2}}^{\top}\mathbf{h}_{0,t}^{-2}\circ(\partial_{\theta}\widehat{\mathbf{h}}_{0,t}-\partial_{\theta}\mathbf{h}_{0,t})-(\widehat{\mathbf{h}}_{0,t}^{-1}-\mathbf{h}_{0,t}^{-1})^{\top}\partial_{\theta}\widehat{\mathbf{h}}_{0,t}$
$\displaystyle\qquad\qquad-{\mathbf{h}_{0,t}^{-1}}^{\top}(\partial_{\theta}\widehat{\mathbf{h}}_{0,t}-\partial_{\theta}\mathbf{h}_{0,t})\bigg{\\}}$

and

$\displaystyle\mathbb{E}\left|\partial_{\theta}\widehat{L}_{T}(\boldsymbol{\theta}_{0})-\partial_{\theta}L_{T}(\boldsymbol{\theta}_{0})\right|$
$\displaystyle\leq\frac{C}{T}\sum_{t=1}^{T}\left\\{\beta_{T}+4\beta_{T}+\beta_{T}+\beta_{T}+\beta_{T}\right\\}$
$\displaystyle=O(\beta_{T}).$ (S3)

Moreover, we have

$\displaystyle\partial_{\theta}L_{T}(\boldsymbol{\theta}_{0})-\partial_{\theta}L_{0,T}(\boldsymbol{\theta}_{0})$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}{(\mathbf{f}_{t}^{2}-\mathbf{h}_{0,t})^{\top}(\mathbf{h}_{0,t}^{-2}\circ\partial_{\theta}\mathbf{h}_{0,t})}$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}{(\mbox{\boldmath$\varepsilon$}_{t}^{2}-\boldsymbol{1})^{\top}(\mathbf{h}_{0,t}^{-1}\circ\partial_{\theta}\mathbf{h}_{0,t})}$

and

$\displaystyle\mathbb{E}[(\partial_{\theta}L_{T}(\boldsymbol{\theta}_{0})-\partial_{\theta}L_{0,T}(\boldsymbol{\theta}_{0}))^{2}]$
$\displaystyle=\frac{1}{T^{2}}\sum_{t=1}^{T}{\operatorname{tr}{\left((\mathbb{E}[\mbox{\boldmath$\varepsilon$}_{t}^{2}{\mbox{\boldmath$\varepsilon$}_{t}^{2}}^{\top}]-\boldsymbol{1}\boldsymbol{1}^{\top})\mathbb{E}\left[(\mathbf{h}_{0,t}^{-1}\circ\partial_{\theta}\mathbf{h}_{0,t})(\mathbf{h}_{0,t}^{-1}\circ\partial_{\theta}\mathbf{h}_{0,t})^{\top}\right]\right)}}$
$\displaystyle=O\left(\frac{1}{T}\right),$ (S4)

where the last equality is due to Lemma S3.4. By combining (S3.1) and (S3.1),
we complete the proof. $\blacksquare$

######  Lemma S3.6.

Under the assumptions of Theorem S3.1, we have

$\partial_{\boldsymbol{\theta}}^{2}\widehat{L}_{T}(\boldsymbol{\theta}^{*})\xrightarrow{p}\partial_{\boldsymbol{\theta}}^{2}L_{0,T}(\boldsymbol{\theta}_{0}),$

and $-\partial_{\boldsymbol{\theta}}^{2}L_{0,T}(\boldsymbol{\theta}_{0})$ is a
positive definite matrix.

Proof of Lemma S3.6. We have

$\displaystyle\partial_{\boldsymbol{\theta}}^{2}\widehat{L}_{T}(\boldsymbol{\theta}^{*})-\partial_{\boldsymbol{\theta}}^{2}L_{0,T}(\boldsymbol{\theta}_{0})$
$\displaystyle=\left(\partial_{\boldsymbol{\theta}}^{2}\widehat{L}_{T}(\boldsymbol{\theta}^{*})-\partial_{\boldsymbol{\theta}}^{2}\widehat{L}_{T}(\boldsymbol{\theta}_{0})\right)+\left(\partial_{\boldsymbol{\theta}}^{2}\widehat{L}_{T}(\boldsymbol{\theta}_{0})-\partial_{\boldsymbol{\theta}}^{2}L_{0,T}(\boldsymbol{\theta}_{0})\right)$
$\displaystyle=\text{(I)}+\text{(II)}.$

By the mean value theorem and (S2), we have
$\text{(I)}=\partial_{\boldsymbol{\theta}}^{3}\widehat{L}_{T}(\boldsymbol{\theta}^{**})(\boldsymbol{\theta}^{*}-\boldsymbol{\theta}_{0})$
and
$\boldsymbol{\theta}^{*}-\boldsymbol{\theta}_{0}\leq\big{|}\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\big{|}\xrightarrow{p}0$,
respectively. Thus, it is enough to show
$\partial_{\boldsymbol{\theta}}^{3}\widehat{L}_{T}(\boldsymbol{\theta}^{**})=O_{p}(1)$.
Denote
$\partial_{ijk}^{3}=\partial_{\theta_{i}}\partial_{\theta_{j}}\partial_{\theta_{j}}$.
Then, for all $i,j,k\leq\dim(\boldsymbol{\theta})$, we have

$\displaystyle\partial_{ijk}^{3}\widehat{L}_{T}(\boldsymbol{\theta})$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\bigg{\\{}\left(\widehat{\mathbf{h}}_{t}-\widehat{\mathbf{f}}_{t}^{2}\right)^{\top}\widehat{\mathbf{h}}_{t}^{\circ-2}\circ\partial_{ijk}^{3}\widehat{\mathbf{h}}_{t}$
$\displaystyle\qquad\qquad+\left(2\widehat{\mathbf{f}}_{t}^{2}-\widehat{\mathbf{h}}_{t}\right)^{\top}\widehat{\mathbf{h}}_{t}^{-3}\circ\left(\partial_{i}\widehat{\mathbf{h}}_{t}\circ\partial_{jk}\widehat{\mathbf{h}}_{t}+\partial_{j}\widehat{\mathbf{h}}_{t}\circ\partial_{ki}\widehat{\mathbf{h}}_{t}+\partial_{k}\widehat{\mathbf{h}}_{t}\circ\partial_{ij}\widehat{\mathbf{h}}_{t}\right)$
$\displaystyle\qquad\qquad+\left(2\widehat{\mathbf{h}}_{t}-6\widehat{\mathbf{f}}_{t}^{2}\right)^{\top}\widehat{\mathbf{h}}_{t}^{-4}\circ\partial_{i}\widehat{\mathbf{h}}_{t}\circ\partial_{j}\widehat{\mathbf{h}}_{t}\circ\partial_{k}\widehat{\mathbf{h}}_{t}\bigg{\\}}.$

By
$\mathbb{E}\big{[}\widehat{\mathbf{f}}_{t}\widehat{\mathbf{f}}_{t}^{\top}\big{]}<\infty$
and Lemma S3.4, we can show
$\mathbb{E}\big{|}\partial_{ijk}^{3}\widehat{L}_{T}(\theta^{**})\big{|}<\infty$.
Thus, $\text{(I)}\xrightarrow{p}0$. Similar to the proof of Lemma S3.5, we can
show

(II)
$\displaystyle=\partial_{ij}^{2}\widehat{L}_{T}(\boldsymbol{\theta}_{0})-\partial_{ij}^{2}L_{0,T}(\boldsymbol{\theta}_{0})$
$\displaystyle=\partial_{ij}^{2}\widehat{L}_{T}(\boldsymbol{\theta}_{0})-\partial_{ij}^{2}L_{T}(\boldsymbol{\theta}_{0})+\partial_{ij}^{2}L_{T}(\boldsymbol{\theta}_{0})-\partial_{ij}^{2}L_{0,T}(\boldsymbol{\theta}_{0})$
$\displaystyle=O_{p}(\beta_{T})+O_{p}\left(\frac{1}{\sqrt{T}}\right).$

Moreover, since $\mathbf{h}_{t}^{-2}>\mathbf{0}$ and $\mathbf{f}_{t}^{2}$’s
are non-degenerate, we can show

$-\partial_{\boldsymbol{\theta}}^{2}L_{0,T}(\boldsymbol{\theta}_{0})=\frac{1}{T}\sum_{t=1}^{T}\sum_{i=1}^{r}{{h_{0,it}^{-2}}\left(\partial_{\boldsymbol{\theta}}h_{0,it}\
{\partial_{\boldsymbol{\theta}}h_{0,it}}^{\top}\right)}\succ 0.$

$\blacksquare$

Proof of Theorem S3.1. By the mean value theorem, there is
$\boldsymbol{\theta}^{*}$ between $\widehat{\boldsymbol{\theta}}$ and
$\boldsymbol{\theta}_{0}$ such that

$\displaystyle\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}$
$\displaystyle=\partial^{2}_{\boldsymbol{\theta}}\widehat{L}_{T}(\boldsymbol{\theta}^{*})^{-1}\left(\partial_{\boldsymbol{\theta}}\widehat{L}_{T}(\widehat{\boldsymbol{\theta}})-\partial_{\boldsymbol{\theta}}\widehat{L}_{T}(\boldsymbol{\theta}_{0})\right)$
$\displaystyle=\partial^{2}_{\boldsymbol{\theta}}\widehat{L}_{T}(\boldsymbol{\theta}^{*})^{-1}\left(\partial_{\boldsymbol{\theta}}L_{0,T}(\boldsymbol{\theta}_{0})-\partial_{\boldsymbol{\theta}}\widehat{L}_{T}(\boldsymbol{\theta}_{0})\right),$

where the last equality is by the fact that
$\partial_{\boldsymbol{\theta}}\widehat{L}_{T}(\widehat{\boldsymbol{\theta}})=\partial_{\boldsymbol{\theta}}L_{0,T}(\boldsymbol{\theta}_{0})=0$.
Then, by Lemmas S3.5 and S3.6, we have

$\displaystyle\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}=O_{p}\left(\beta_{T}+\frac{1}{\sqrt{T}}\right),$

and thus, the result of Theorem S2.1 completes the proof.$\blacksquare$

##  Appendix S4 Proof of Theorem 4.1

######  Assumption 3.

  1. (a)

The sample covariance estimator satisfies the concentration inequality, for
any given $a>0$,

$\mathbb{P}\left\\{\max_{ij}\Big{|}\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\Sigma\mkern-1.5mu}\mkern 1.5mu}_{ij}-\mkern
1.5mu\overline{\mkern-1.5mu\Sigma\mkern-1.5mu}\mkern 1.5mu_{ij}\Big{|}\geq
C_{a}\sqrt{\frac{\log p}{T}}\right\\}\leq p^{-a},$

where $C_{a}$ is a constant only depending on $a$.

  2. (b)

There is a constant $C$ such that

$\frac{1}{r}\max_{i\leq p}\sum_{j=1}^{r}V_{ij}^{2}\leq C.$

  3. (c)

The smallest eigenvalue of $\boldsymbol{\Sigma}_{u}$ stays away from $0$, and
$\big{|}\Sigma_{u,ij}\big{|}\leq C$ for all $i$, $j$.

  4. (d)

$p=o(T^{2})$.

######  Remark 7.

Assumption 3(a) is the sub-Gaussian condition. When investigating the high-
dimensional inferences, the sub-Gaussian condition is essential. Recently,
several robust variance estimation methods have been developed, which can
satisfy the sub-Gaussian condition Assumption 3(a) under the finite forth
moment condition (Catoni, , 2012; Fan and Kim, , 2018; Minsker, , 2018; Shin
et al., 2021a, ). Thus, the sub-Gaussian condition Assumption 3(a) is not
restrictive. On the other hand, we can investigate the asymptotic behaviors
under sub-Weilbull conditions (Kuchibhotla and Chakrabortty, , 2018;
Vladimirova et al., , 2020).

######  Remark 8.

Assumption 3(b) is called the incoherence condition, which is widely assumed
in analyzing the low-rank matrix (Candès and Recht, , 2009; Fan et al., ,
2017). The basic intuition is that the factor loading matrix $\mathbf{V}$ is
not to be sparse. That is, the factor affects almost all the stock returns.
Thus, under the factor model, the incoherence condition is acceptable.

######  Lemma S4.1.

Under the assumptions of Theorem 4.1, we have

$\displaystyle\left\|\widehat{\mathbf{V}}\widehat{\boldsymbol{\Sigma}}_{f,t+1}\widehat{\mathbf{V}}^{\top}-\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}\mathbf{V}^{\top}\right\|_{\boldsymbol{\Sigma}_{t+1}}=O_{p}\left(\frac{1}{\sqrt{T}}+\frac{\sqrt{p}}{T}+\frac{s_{p}^{2}}{p\sqrt{p}}+\frac{s_{p}}{p}\right).$
(S1)

Proof of Lemma S4.1. Without loss of generality, we assume that $\mathbf{O}$
is the identity matrix. Similar to the proofs of Theorem 4.2 in Kim et al.,
(2018), we have

$\displaystyle\left\|\widehat{\mathbf{V}}\widehat{\boldsymbol{\Sigma}}_{f,t+1}\widehat{\mathbf{V}}^{\top}-\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}\mathbf{V}^{\top}\right\|_{\boldsymbol{\Sigma}_{t+1}}$
$\displaystyle\leq\left\|\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}\left(\boldsymbol{\Sigma}_{f,t+1}^{-\frac{1}{2}}\widehat{\boldsymbol{\Sigma}}_{f,t+1}\boldsymbol{\Sigma}_{f,t+1}^{-\frac{1}{2}}-\mathbf{I}_{r}\right)\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}\mathbf{V}^{\top}\right\|_{\boldsymbol{\Sigma}_{t+1}}$
$\displaystyle\quad+\left\|(\widehat{\mathbf{V}}-\mathbf{V})\widehat{\boldsymbol{\Sigma}}_{f,t+1}(\widehat{\mathbf{V}}-\mathbf{V})^{\top}\right\|_{\boldsymbol{\Sigma}_{t+1}}+2\left\|\mathbf{V}\widehat{\boldsymbol{\Sigma}}_{f,t+1}(\widehat{\mathbf{V}}-\mathbf{V})^{\top}\right\|_{\boldsymbol{\Sigma}_{t+1}}$
$\displaystyle=\text{(I)}+\text{(II)}+\text{(III)}.$

First consider (I). By the Sherman-Morrison-Woodbury formula (Higham, , 2002),
we have

$\displaystyle(\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}\mathbf{V}^{\top}+\boldsymbol{\Sigma}_{u})^{-1}=\boldsymbol{\Sigma}_{u}^{-1}-\boldsymbol{\Sigma}_{u}^{-1}\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}(\mathbf{I}_{r}+\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}\mathbf{V}^{\top}\boldsymbol{\Sigma}_{u}^{-1}\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}})^{-1}\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}\mathbf{V}^{\top}\boldsymbol{\Sigma}_{u}^{-1}$

and by denoting
$\mathbf{X}=\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}\mathbf{V}^{\top}\boldsymbol{\Sigma}_{u}^{-1}\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}$,

$\displaystyle\left\|\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}\mathbf{V}^{\top}\boldsymbol{\Sigma}_{t+1}^{-1}\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}\right\|$
$\displaystyle=\left\|\mathbf{X}-\mathbf{X}(\mathbf{I}_{r}+\mathbf{X})^{-1}\mathbf{X}\right\|$
$\displaystyle=\left\|\mathbf{X}(\mathbf{I}_{r}+\mathbf{X})^{-1}\right\|$
$\displaystyle=\left\|\mathbf{I}_{r}-(\mathbf{I}_{r}+\mathbf{X})^{-1}\right\|$
$\displaystyle\leq 2.$

Then, since $\|\mathbf{A}\mathbf{B}\|_{F}\leq\|\mathbf{A}\|\|\mathbf{B}\|_{F}$
and
$\big{\|}\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}\mathbf{V}^{\top}\boldsymbol{\Sigma}_{t+1}^{-\frac{1}{2}}\big{\|}\leq\sqrt{2}$,
we have

(I)
$\displaystyle=\frac{1}{\sqrt{p}}\left\|\boldsymbol{\Sigma}_{t+1}^{-\frac{1}{2}}\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}\left(\boldsymbol{\Sigma}_{f,t+1}^{-\frac{1}{2}}\widehat{\boldsymbol{\Sigma}}_{f,t+1}\boldsymbol{\Sigma}_{f,t+1}^{-\frac{1}{2}}-\mathbf{I}_{r}\right)\boldsymbol{\Sigma}_{f,t+1}^{\frac{1}{2}}\mathbf{V}^{\top}\boldsymbol{\Sigma}_{t+1}^{-\frac{1}{2}}\right\|_{F}$
$\displaystyle\leq\frac{2}{\sqrt{p}}\left\|\boldsymbol{\Sigma}_{f,t+1}^{-\frac{1}{2}}\widehat{\boldsymbol{\Sigma}}_{f,t+1}\boldsymbol{\Sigma}_{f,t+1}^{-\frac{1}{2}}-\mathbf{I}_{r}\right\|_{F}$
$\displaystyle=\frac{2}{\sqrt{p}}\left(\sum_{i=1}{\frac{\left|\widehat{h}_{it+1}(\widehat{\boldsymbol{\theta}})-h_{0,it+1}\right|^{2}}{h_{0,it+1}^{2}}}\right)^{\frac{1}{2}}$
$\displaystyle=O_{p}\left(\frac{1}{\sqrt{pT}}+\frac{\sqrt{s_{p}}}{p}\right),$

where the last equality is by the fact that

$\displaystyle\widehat{h}_{it+1}(\widehat{\boldsymbol{\theta}})-h_{0,it+1}$
$\displaystyle=\widehat{h}_{it+1}(\widehat{\boldsymbol{\theta}})-\widehat{h}_{it+1}(\boldsymbol{\theta}_{0})+\widehat{h}_{it+1}(\boldsymbol{\theta}_{0})-h_{0,it+1}$
$\displaystyle=O_{p}\left(\frac{1}{\sqrt{T}}+\sqrt{\frac{s_{p}}{p}}\right).$

For (II) and (III), similarly, by Theorem S2.1, we can show

(II)
$\displaystyle=\frac{1}{\sqrt{p}}\left\|\boldsymbol{\Sigma}_{t+1}^{-\frac{1}{2}}(\widehat{\mathbf{V}}-\mathbf{V})\widehat{\boldsymbol{\Sigma}}_{f,t+1}(\widehat{\mathbf{V}}-\mathbf{V})^{\top}\boldsymbol{\Sigma}_{t+1}^{-\frac{1}{2}}\right\|_{F}$
$\displaystyle\leq\frac{1}{\sqrt{p}}\left\|\boldsymbol{\Sigma}_{t+1}^{-1}\right\|\left\|\widehat{\boldsymbol{\Sigma}}_{f,t+1}\right\|_{F}\left\|\widehat{\mathbf{V}}-\mathbf{V}\right\|^{2}$
$\displaystyle=O_{p}\left(\frac{\sqrt{p}}{T}+\frac{s_{p}^{2}}{p\sqrt{p}}\right)$

and

(III)
$\displaystyle=\frac{1}{\sqrt{p}}\left\|\boldsymbol{\Sigma}_{t+1}^{-\frac{1}{2}}\mathbf{V}\widehat{\boldsymbol{\Sigma}}_{f,t+1}(\widehat{\mathbf{V}}-\mathbf{V})^{\top}\boldsymbol{\Sigma}_{t+1}^{-\frac{1}{2}}\right\|_{F}$
$\displaystyle=\frac{\sqrt{2}}{\sqrt{p}}\left\|\widehat{\boldsymbol{\Sigma}}_{f,t+1}^{\frac{1}{2}}(\widehat{\mathbf{V}}-\mathbf{V})^{\top}\boldsymbol{\Sigma}_{t+1}^{-\frac{1}{2}}\right\|_{F}$
$\displaystyle\leq\frac{C}{\sqrt{p}}\left\|\widehat{\boldsymbol{\Sigma}}_{f,t+1}^{\frac{1}{2}}\right\|_{F}\left\|\widehat{\mathbf{V}}-\mathbf{V}\right\|$
$\displaystyle=O_{p}\left(\frac{1}{\sqrt{T}}+\frac{s_{p}}{p}\right).$

$\blacksquare$

######  Lemma S4.2.

Under the assumptions of Theorem 4.1, suppose that
$\|\widehat{\boldsymbol{\Sigma}}_{u}-\boldsymbol{\Sigma}_{u}\|_{\max}=O_{p}(\gamma_{T})$,
and the thresholding level satisfies the condition $\tau_{T}=C\gamma_{T}$ such
that
$|\widehat{\Sigma}_{u,ij}-\Sigma_{u,ij}|<\tau_{T}(\widehat{\Sigma}_{u,ii}\widehat{\Sigma}_{u,jj})^{1/2}/2$.
Then we have

$\left\|{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})-\boldsymbol{\Sigma}_{u}\right\|=O_{p}\left(s_{p}\gamma_{T}^{1-q}\right).$

Proof of Lemma S4.2. Let
$\tau_{ij}=\tau_{T}(\widehat{\Sigma}_{u,ii}\widehat{\Sigma}_{u,jj})^{1/2}$.
Similar to the proofs of Section 3 in Fan and Kim,  (2018), under the event
$|\widehat{\Sigma}_{u,ij}-\Sigma_{u,ij}|<\tau_{ij}/2$, we have

$\displaystyle\left\|{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})-\boldsymbol{\Sigma}_{u}\right\|$
$\displaystyle\leq\left\|{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})-\boldsymbol{\Sigma}_{u}\right\|_{\infty}$
$\displaystyle=\max_{i}\sum_{j=1}^{p}{\left|s_{ij}(\widehat{\Sigma}_{u,ij})\mathbbm{1}_{|\widehat{\Sigma}_{u,ij}|\geq\tau_{ij}}-\Sigma_{u,ij}\mathbbm{1}_{|\Sigma_{u,ij}|\geq\tau_{ij}}-\Sigma_{u,ij}\mathbbm{1}_{|\Sigma_{u,ij}|<\tau_{ij}}\right|}$
$\displaystyle\leq\max_{i}\sum_{j=1}^{p}\bigg{\\{}\left|s_{ij}(\widehat{\Sigma}_{u,ij})-\Sigma_{u,ij}\right|\mathbbm{1}_{|\widehat{\Sigma}_{u,ij}|\geq\tau_{ij}}+\left|\Sigma_{u,ij}\right|\left|\mathbbm{1}_{|\widehat{\Sigma}_{u,ij}|\geq\tau_{ij}}-\mathbbm{1}_{|\Sigma_{u,ij}|\geq\tau_{ij}}\right|$
$\displaystyle\qquad\qquad\quad+\left|\Sigma_{u,ij}\right|\mathbbm{1}_{|\Sigma_{u,ij}|<\tau_{ij}}\bigg{\\}}$
$\displaystyle\leq\max_{i}\sum_{j=1}^{p}\bigg{\\{}\frac{3}{2}\tau_{ij}\mathbbm{1}_{|\Sigma_{u,ij}|\geq\frac{1}{2}\tau_{ij}}+\left|\Sigma_{u,ij}\right|\mathbbm{1}_{|\Sigma_{u,ij}|\leq\frac{3}{2}\tau_{ij}}+\left|\Sigma_{u,ij}\right|^{q}\tau_{ij}^{1-q}\bigg{\\}}$
$\displaystyle\leq
C\max_{i}\sum_{j=1}^{p}{|\Sigma_{u,ij}|^{q}\tau_{ij}^{1-q}}$
$\displaystyle=O\left(s_{p}\gamma_{T}^{1-q}\right),$

where the last equality is due to the sparse condition. $\blacksquare$

Proof of Theorem 4.1. By Assumption (a)(a), we have

$\displaystyle\left\|\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu}-\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu\right\|_{\max}=O_{p}\left(\sqrt{\frac{\log{p}}{T}}\right).$ (S2)

We denote $\widehat{\mathbf{v}}_{i}$ and $\mathbf{v}_{i}$ the $i^{th}$ column
vectors of $\widehat{\mathbf{V}}$ and $\mathbf{V}$, respectively. Then by
Theorem 2.1 in Fan et al.,  (2017),

$\displaystyle\sum_{i=1}^{r}{\left\|\widehat{\mathbf{v}}_{i}\widehat{\mathbf{v}}_{i}^{\top}-\mathbf{v}_{i}\mathbf{v}_{i}^{\top}\right\|_{\max}}$
$\displaystyle\leq\sum_{i=1}^{r}{\left\\{\left\|(\widehat{\mathbf{v}}_{i}-\mathbf{v}_{i})\widehat{\mathbf{v}}_{i}^{\top}\right\|_{\max}+\left\|\mathbf{v}_{i}(\widehat{\mathbf{v}}_{i}-\mathbf{v}_{i})^{\top}\right\|_{\max}\right\\}}$
$\displaystyle\leq\sum_{i=1}^{r}{\left\\{\left\|\widehat{\mathbf{v}}_{i}-\mathbf{v}_{i}\right\|_{\infty}\left\|\widehat{\mathbf{v}}_{i}-\mathbf{v}_{i}+\mathbf{v}_{i}\right\|_{\infty}+\left\|\widehat{\mathbf{v}}_{i}-\mathbf{v}_{i}\right\|_{\infty}\left\|\mathbf{v}_{i}\right\|_{\infty}\right\\}}$
$\displaystyle\leq\sum_{i=1}^{r}{\left\\{\left\|\widehat{\mathbf{v}}_{i}-\mathbf{v}_{i}\right\|_{\infty}^{2}+2\left\|\widehat{\mathbf{v}}_{i}-\mathbf{v}_{i}\right\|_{\infty}\left\|\mathbf{v}_{i}\right\|_{\infty}\right\\}}$
$\displaystyle\leq C\left(\frac{\left\|\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu}-\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu+\boldsymbol{\Sigma}_{u}\right\|_{\infty}^{2}}{p^{2}\delta_{r}^{2}}+\frac{\left\|\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu}-\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu+\boldsymbol{\Sigma}_{u}\right\|_{\infty}}{p\delta_{r}}\right)$
$\displaystyle=O_{p}\left(\sqrt{\frac{\log{p}}{T}}+\frac{s_{p}}{p}\right),$

where the last equality is due to $\|\mathbf{A}\|_{\infty}\leq
p\|\mathbf{A}\|_{\max}$. Thus, we have

$\displaystyle\left\|\widehat{\mathbf{V}}\widehat{\boldsymbol{\Sigma}}_{f,t+1}\widehat{\mathbf{V}}^{\top}-\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}\mathbf{V}^{\top}\right\|_{\max}$
$\displaystyle\leq\sum_{i=1}^{r}{\left\\{\left\|\widehat{h}_{it+1}\left(\widehat{\mathbf{v}}_{i}\widehat{\mathbf{v}}_{i}^{\top}-\mathbf{v}_{i}\mathbf{v}_{i}^{\top}\right)\right\|_{\max}+\left\|(\widehat{h}_{it+1}-h_{0,it+1})\mathbf{v}_{i}\mathbf{v}_{i}^{\top}\right\|_{\max}\right\\}}$
$\displaystyle\leq\left\|\widehat{\mathbf{h}}_{t+1}-\mathbf{h}_{0,t+1}+\mathbf{h}_{0,t+1}\right\|_{\infty}\sum_{i=1}^{r}{\left\|\widehat{\mathbf{v}}_{i}\widehat{\mathbf{v}}_{i}^{\top}-\mathbf{v}_{i}\mathbf{v}_{i}^{\top}\right\|_{\max}}+\left\|\widehat{\mathbf{h}}_{t+1}-\mathbf{h}_{0,t+1}\right\|_{\infty}\sum_{i=1}^{r}{\left\|\mathbf{v}_{i}\mathbf{v}_{i}^{\top}\right\|_{\max}}$
$\displaystyle\leq
O_{p}\left(\sqrt{\frac{\log{p}}{T}}+\frac{s_{p}}{p}\right)+O_{p}\left(\frac{1}{\sqrt{T}}+\sqrt{\frac{s_{p}}{p}}\right).$
(S3)

By combining (S2) and (S4), we have

$\displaystyle\left\|\widehat{\boldsymbol{\Sigma}}_{u}-\boldsymbol{\Sigma}_{u}\right\|_{\max}$
$\displaystyle\leq\left\|\widehat{\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu}-\mkern
1.5mu\overline{\mkern-1.5mu\boldsymbol{\Sigma}\mkern-1.5mu}\mkern
1.5mu\right\|_{\max}+\left\|\widehat{\mathbf{V}}\widehat{\boldsymbol{\Sigma}}_{f,t+1}\widehat{\mathbf{V}}^{\top}-\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}\mathbf{V}^{\top}\right\|_{\max}$
$\displaystyle=O_{p}\left(\sqrt{\frac{\log{p}}{T}}+\sqrt{\frac{s_{p}}{p}}\right),$
(S4)

which leads to $\gamma_{T}=C\big{(}\sqrt{\log{p}/T}+\sqrt{s_{p}/p}\big{)}$ and
$\|\widehat{\boldsymbol{\Sigma}}_{t+1}-\boldsymbol{\Sigma}_{t+1}\|_{\max}=O_{p}(\gamma_{T})$.
Combining Lemmas S4.1 and S4.2, and (S4), we have

$\displaystyle\left\|\widehat{\boldsymbol{\Sigma}}_{t+1}-\boldsymbol{\Sigma}_{t+1}\right\|_{\boldsymbol{\Sigma}_{t+1}}$
$\displaystyle\leq\left\|\widehat{\mathbf{V}}\widehat{\boldsymbol{\Sigma}}_{f,t+1}\widehat{\mathbf{V}}^{\top}-\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}\mathbf{V}^{\top}\right\|_{\boldsymbol{\Sigma}_{t+1}}+\left\|{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})-\boldsymbol{\Sigma}_{u}\right\|_{\boldsymbol{\Sigma}_{t+1}}$
$\displaystyle=O_{p}\left(\frac{\sqrt{p}}{T}+s_{p}\gamma_{T}^{1-q}\right),$

where the last equality is by the fact

$\displaystyle\left\|{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})-\boldsymbol{\Sigma}_{u}\right\|_{\boldsymbol{\Sigma}_{t+1}}$
$\displaystyle\leq\frac{1}{\sqrt{p}}\left\|{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})-\boldsymbol{\Sigma}_{u}\right\|\left\|\boldsymbol{\Sigma}_{t+1}^{-1}\right\|_{F}$
$\displaystyle\leq\left\|{\cal
T}(\widehat{\boldsymbol{\Sigma}}_{u})-\boldsymbol{\Sigma}_{u}\right\|\left\|\boldsymbol{\Sigma}_{t+1}^{-1}\right\|$
$\displaystyle=O_{p}\left(s_{p}\gamma_{T}^{1-q}\right).$

$\blacksquare$

##  Appendix S5 Proof of Theorem 4.2

######  Assumption 4.

  1. (a)

The estimated standardized return
$\widehat{x}_{t}=\mathbf{w}^{\top}(\mathbf{y}_{t}-\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu)/(\mathbf{w}^{\top}\widehat{\boldsymbol{\Sigma}}_{t}\mathbf{w})^{1/2}$
satisfies the concentration inequality, for any given $a>0$,

$\max_{t\leq T}\mathbb{P}\left\\{\left|\widehat{x}_{t}-x_{t}\right|\geq
C_{a}\sqrt{\log{T}}\left(\sqrt{\frac{\log{p}}{T}}+\sqrt{\frac{s_{p}}{p}}\right)\right\\}\leq
T^{-a},$

where
$x_{t}=\mathbf{w}^{\top}(\mathbf{y}_{t}-\boldsymbol{\mu})/(\mathbf{w}^{\top}\boldsymbol{\Sigma}_{t}\mathbf{w})^{1/2}$
and $C_{a}$ is a constant only depending on $a$.

  2. (b)

The cumulative density function of $x_{t}$ and its inverse function are
continuous.

  3. (c)

The variance of the portfolio return $r_{t}$ is bounded and strictly positive.

######  Remark 9.

The sub-Gaussian condition Assumption (a)(a) is true under some sub-Gaussian
condition for $\mathbf{y}_{t}$. Assumption (b)(b) is usually imposed to study
quantile estimations (Chen and Tang, , 2005).

Proof of Theorem 4.2. Consider the parametric VaR estimator case. By the
definition of VaR in (2.3), we have

$\displaystyle\left|\widehat{\text{VaR}}_{\alpha,t+1}-\text{VaR}_{\alpha,t+1}\right|$
$\displaystyle\leq
C\left|c_{\alpha}\sqrt{\mathbf{w}^{\top}\widehat{\boldsymbol{\Sigma}}_{t+1}\mathbf{w}}-c_{\alpha}\sqrt{\mathbf{w}^{\top}\boldsymbol{\Sigma}_{t+1}\mathbf{w}}\right|+C\left|\mathbf{w}^{\top}(\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu-\boldsymbol{\mu})\right|.$

Then, since $\|\mkern 1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu-\boldsymbol{\mu}\|_{\max}^{2}=O_{p}(\log{p}/T)$ and $\|\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu-\boldsymbol{\mu}\|^{2}=O_{p}(p/T)$, we have

$\displaystyle\left(\mathbf{w}^{\top}(\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu-\boldsymbol{\mu})\right)^{2}=O_{p}\left(\min{\left\\{\frac{\log{p}}{T},\|\mathbf{w}\|^{2}\frac{p}{T}\right\\}}\right).$
(S1)

Consider
$\Big{|}c_{\alpha}\sqrt{\mathbf{w}^{\top}\widehat{\boldsymbol{\Sigma}}_{t+1}\mathbf{w}}-c_{\alpha}\sqrt{\mathbf{w}^{\top}\boldsymbol{\Sigma}_{t+1}\mathbf{w}}\Big{|}$.
We have

$\displaystyle\left|c_{\alpha}\sqrt{\mathbf{w}^{\top}\widehat{\boldsymbol{\Sigma}}_{t+1}\mathbf{w}}-c_{\alpha}\sqrt{\mathbf{w}^{\top}\boldsymbol{\Sigma}_{t+1}\mathbf{w}}\right|$
$\displaystyle\leq
C\left|\mathbf{w}^{\top}\left(\widehat{\boldsymbol{\Sigma}}_{t+1}-\boldsymbol{\Sigma}_{t+1}\right)\mathbf{w}\right|$
$\displaystyle\leq
C\sum_{i,j=1}^{p}{|w_{i}w_{j}|}\left\|\widehat{\boldsymbol{\Sigma}}_{t+1}-\boldsymbol{\Sigma}_{t+1}\right\|_{\max}$
$\displaystyle=O_{p}\left(\sqrt{\frac{\log{p}}{T}}+\sqrt{\frac{s_{p}}{p}}\right),$
(S2)

where the last is due to Theorem 4.1 and $\|\mathbf{w}\|_{1}\leq C$. Combining
(S1) and (S5), we have

$\displaystyle\left|\widehat{\text{VaR}}_{\alpha,t+1}-\text{VaR}_{\alpha,t+1}\right|=O_{p}\left(\sqrt{\frac{\log{p}}{T}}+\sqrt{\frac{s_{p}}{p}}\right).$

Consider the non-parametric $\sigma$-based VaR estimator case. Define

$\displaystyle\widehat{F}_{T}(x)=\frac{1}{T}\sum_{t=1}^{T}{\mathbbm{1}_{\\{\widehat{x}_{t}\leq
x\\}}},\qquad F_{T}(x)=\frac{1}{T}\sum_{t=1}^{T}{\mathbbm{1}_{\\{x_{t}\leq
x\\}}},$

where $x$ is an $\alpha$-quantile value, and $F(x)$ as a cumulative
distribution function of $x_{t}$. Then the expected value of the absolute
difference between $\widehat{F}_{T}(x)$ and $F_{T}(x)$ is

$\displaystyle\mathbb{E}\left|\widehat{F}_{T}(x)-F_{T}(x)\right|\leq\frac{1}{T}\sum_{t=1}^{T}{\big{(}\mathbb{P}\\{\widehat{x}_{t}\leq
x,x_{t}>x\\}+\mathbb{P}\\{\widehat{x}_{t}>x,x_{t}\leq x\\}\big{)}}.$

Let $\vartheta_{T}=\sqrt{\log{p}/T}+\sqrt{s_{p}/p}$. Then, since Assumption
(a)(a) implies $|\widehat{x}_{t}-x_{t}|=O_{p}(\vartheta_{T}\sqrt{\log{T}})$,
we have, for large enough $C$,

$\displaystyle\mathbb{P}\\{\widehat{x}_{t}\leq x,x_{t}>x\\}$
$\displaystyle\leq\mathbb{P}\\{x_{t}\leq x+|\widehat{x}_{t}-x_{t}|,x_{t}>x\\}$
$\displaystyle=\mathbb{P}\\{x<x_{t}\leq x+|\widehat{x}_{t}-x_{t}|\\}$
$\displaystyle\leq\mathbb{P}\\{x<x_{t}\leq
x+C\vartheta_{T}\sqrt{\log{T}}\\}+\mathbb{P}\\{|\widehat{x}_{t}-x_{t}|>C\vartheta_{T}\sqrt{\log{T}}\\}$
$\displaystyle\leq F(x+C\vartheta_{T}\sqrt{\log{T}})-F(x)+\frac{1}{\sqrt{T}}$
$\displaystyle=O\left(\vartheta_{T}\sqrt{\log{T}}+\frac{1}{\sqrt{T}}\right),$

where the last equality is due to Assumption (b)(b). Similarly, the bound for
$\mathbb{P}\\{\widehat{x}_{t}>x,x_{t}\leq x\\}$ can be found. Then
$\big{|}\widehat{F}_{T}(x)-F_{T}(x)\big{|}=O_{p}\big{(}\vartheta_{T}\sqrt{\log{T}}+1/\sqrt{T}\big{)}$.
By the Dvoretzky–Kiefer-Wolfowitz inequality (Dvoretzky et al., , 1956), we
have $\big{|}F_{T}(x)-F(x)\big{|}=O_{p}\big{(}1/\sqrt{T}\big{)}$. Thus,

$\displaystyle\big{|}\widehat{F}_{T}(x)-F(x)\big{|}$
$\displaystyle\leq\big{|}\widehat{F}_{T}(x)-F_{T}(x)\big{|}+\big{|}F_{T}(x)-F(x)\big{|}$
$\displaystyle=O_{p}\big{(}\vartheta_{T}\sqrt{\log{T}}+1/\sqrt{T}\big{)}.$

Define

$\displaystyle F^{-1}(y)=\inf{\big{\\{}x:F(x)>y\big{\\}}}.$

Then the $\lceil\alpha T\rceil$-th smallest value of $\\{x_{t}\\}_{t=1}^{T}$
is bounded as follows:

$\displaystyle\widehat{F}^{-1}_{T}(\alpha)$
$\displaystyle=\inf{\big{\\{}x:\widehat{F}_{T}(x)>\alpha\big{\\}}}$
$\displaystyle\leq\inf{\big{\\{}x:F(x)>\alpha+\big{|}\widehat{F}_{T}(x)-F(x)\big{|}\big{\\}}}$
$\displaystyle\leq
F^{-1}\left(\alpha+\big{|}\widehat{F}_{T}(x)-F(x)\big{|}\right)$

and

$\displaystyle\widehat{F}^{-1}_{T}(\alpha)$
$\displaystyle\geq\inf{\big{\\{}x:F(x)>\alpha-\big{|}\widehat{F}_{T}(x)-F(x)\big{|}\big{\\}}}$
$\displaystyle\geq
F^{-1}\left(\alpha-\big{|}\widehat{F}_{T}(x)-F(x)\big{|}\right).$

Thus, we have

$\displaystyle\big{|}\widehat{F}_{T}^{-1}(\alpha)-F^{-1}(\alpha)\big{|}\leq
C\big{|}\widehat{F}_{T}(x)-F(x)\big{|}.$ (S3)

We also have

$\displaystyle\mathbf{w}^{\top}\widehat{\boldsymbol{\Sigma}}_{t+1}\mathbf{w}$
$\displaystyle=\mathbf{w}^{\top}\boldsymbol{\Sigma}_{t+1}\mathbf{w}+\mathbf{w}^{\top}(\widehat{\boldsymbol{\Sigma}}_{t+1}-\boldsymbol{\Sigma}_{t+1})\mathbf{w}$
$\displaystyle=\mathbf{w}^{\top}\mathbf{V}\boldsymbol{\Sigma}_{f,t+1}\mathbf{V}^{\top}\mathbf{w}+\mathbf{w}^{\top}\boldsymbol{\Sigma}_{u}\mathbf{w}+o_{p}(1)$
$\displaystyle\leq\max_{i\leq
r}{h_{it+1}}\left\|\mathbf{V}^{\top}\mathbf{w}\right\|^{2}+C+o_{p}(1)$
$\displaystyle=O_{p}(1),$ (S4)

where the inequality is due to the gross exposure condition and
$|\Sigma_{u,ij}|\leq C$, and the last equality is due to
$\|\mathbf{V}^{\top}\mathbf{w}\|^{2}=\|\sum_{i=1}^{p}{w_{i}\widetilde{\mathbf{v}}_{i}}\|^{2}\leq\sum_{i=1}^{p}{|w_{i}|\|\widetilde{\mathbf{v}}_{i}\|^{2}}\leq
C$, where $\widetilde{\mathbf{v}}_{i}$ is the $i^{th}$ row vector of
$\mathbf{V}$. Combining results from (S1) to (S5), we have

$\displaystyle\left|\widehat{\text{VaR}}_{\alpha,t+1}-\text{VaR}_{\alpha,t+1}\right|$
$\displaystyle\leq
C\sqrt{\mathbf{w}^{\top}\widehat{\boldsymbol{\Sigma}}_{t+1}\mathbf{w}}\left|\widehat{F}_{T}^{-1}(\alpha)-F^{-1}(\alpha)\right|+C\left|\mathbf{w}^{\top}(\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{y}\mkern-1.5mu}\mkern
1.5mu-\boldsymbol{\mu})\right|$
$\displaystyle\quad+CF^{-1}(\alpha)\left|\sqrt{\mathbf{w}^{\top}\widehat{\boldsymbol{\Sigma}}_{t+1}\mathbf{w}}-\sqrt{\mathbf{w}^{\top}\boldsymbol{\Sigma}_{t+1}\mathbf{w}}\right|$
$\displaystyle=O_{p}\left(\vartheta_{T}\sqrt{\log{T}}\right).$

$\blacksquare$

## References

  * Acerbi and Tasche, (2002) Acerbi, C. and Tasche, D. (2002).  Expected shortfall: a natural coherent alternative to value at risk.  Economic notes, 31(2):379–388.
  * Ahmadi-Javid, (2012) Ahmadi-Javid, A. (2012).  Entropic value-at-risk: A new coherent risk measure.  Journal of Optimization Theory and Applications, 155(3):1105–1123.
  * Ait-Sahalia and Xiu, (2017) Ait-Sahalia, Y. and Xiu, D. (2017).  Using principal component analysis to estimate a high dimensional factor model with high-frequency data.  Journal of Econometrics, 201(2):384–399.
  * Amemiya, (1985) Amemiya, T. (1985).  Advanced econometrics.  Harvard university press.
  * Bai, (2003) Bai, J. (2003).  Inferential theory for factor models of large dimensions.  Econometrica, 71(1):135–171.
  * Bai and Li, (2012) Bai, J. and Li, K. (2012).  Statistical analysis of factor models of high dimension.  The Annals of Statistics, 40(1):436–465.
  * Bai and Ng, (2013) Bai, J. and Ng, S. (2013).  Principal components estimation and identification of static factors.  Journal of Econometrics, 176(1):18–29.
  * Barigozzi and Hallin, (2016) Barigozzi, M. and Hallin, M. (2016).  Generalized dynamic factor models and volatilities: recovering the market volatility shocks.  The Econometrics Journal, 19(1).
  * Bauwens et al., (2006) Bauwens, L., Laurent, S., and Rombouts, J. V. (2006).  Multivariate garch models: a survey.  Journal of Applied Econometrics, 21(1):79–109.
  * Bensaid et al., (1992) Bensaid, B., Lesne, J.-P., Pages, H., and Scheinkman, J. (1992).  Derivative asset pricing with transaction costs 1.  Mathematical Finance, 2(2):63–86.
  * Bickel and Levina, (2008) Bickel, P. J. and Levina, E. (2008).  Covariance regularization by thresholding.  The Annals of Statistics, 36(6):2577–2604.
  * Billio and Pelizzon, (2000) Billio, M. and Pelizzon, L. (2000).  Value-at-risk: a multivariate switching regime approach.  Journal of Empirical Finance, 7(5):531–554.
  * Boivin and Ng, (2006) Boivin, J. and Ng, S. (2006).  Are more data always better for factor analysis?  Journal of Econometrics, 132(1):169–194.
  * Bollerslev, (1986) Bollerslev, T. (1986).  Generalized autoregressive conditional heteroskedasticity.  Journal of Econometrics, 31(3):307–327.
  * Bollerslev, (1990) Bollerslev, T. (1990).  Modelling the coherence in short-run nominal exchange rates: a multivariate generalized arch model.  Review of Economics and Statistics, 72(3):498–505.
  * Broda and Paolella, (2009) Broda, S. A. and Paolella, M. S. (2009).  Chicago: A fast and accurate method for portfolio risk calculation.  Journal of Financial Econometrics, 7(4):412–436.
  * Brooks and Persand, (2003) Brooks, C. and Persand, G. (2003).  Volatility forecasting for risk management.  Journal of Forecasting, 22(1):1–22.
  * Cai and Liu, (2011) Cai, T. and Liu, W. (2011).  Adaptive thresholding for sparse covariance matrix estimation.  Journal of the American Statistical Association, 106(494):672–684.
  * Cajueiro and Tabak, (2004) Cajueiro, D. O. and Tabak, B. M. (2004).  The hurst exponent over time: testing the assertion that emerging markets are becoming more efficient.  Physica A: Statistical Mechanics and its Applications, 336(3-4):521–537.
  * Candès and Recht, (2009) Candès, E. J. and Recht, B. (2009).  Exact matrix completion via convex optimization.  Foundations of Computational Mathematics, 9(6):717.
  * Carhart, (1997) Carhart, M. M. (1997).  On persistence in mutual fund performance.  The Journal of Finance, 52(1):57–82.