EconBase
← Back to paper

Unified Discrete-Time Factor Stochastic Volatility and Continuous-Time Ito Models for Combining Inference Based on Low-Frequency and High-Frequency

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

81,595 characters · 16 sections · 46 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Unified Discrete-Time Factor Stochastic Volatility and Continuous-Time It\^o Models for Combining Inference Based on Low-Frequency and High-Frequency

abstractThis paper introduces unified models for high-dimensional factor-based It\^o process, which can accommodate both continuous-time It\^o diffusion and discrete-time stochastic volatility (SV) models by embedding the discrete SV model in the continuous instantaneous factor volatility process. We call it the SV-It\^o model. Based on the series of daily integrated factor volatility matrix estimators, we propose quasi-maximum likelihood and least squares estimation methods. Their asymptotic properties are established. We apply the proposed method to predict future vast volatility matrix whose asymptotic behaviors are studied. A simulation study is conducted to check the finite sample performance of the proposed estimation and prediction method. An empirical analysis is carried out to demonstrate the advantage of the SV-It\^o model in volatility prediction and portfolio allocation problems.

Keywords: Factor model, high dimensionality, POET, quasi-maximum likelihood estimation, stochastic volatility model

\eject \baselineskip=20pt

Introduction

Volatility analysis with high-frequency data is a vibrant research area in finance. Researchers have devoted to developing volatility estimation methods under the continuous-time frameworks such as It\^{o} diffusion process. Example estimators for daily volatility include: two-time scale realized volatility zhang2005tale, multi-scale realized volatility zhang2006efficient, zhang2011estimating, kernel realized volatility barndorff2008designing, barndorff2011multivariate, pre-averaging realized volatility christensen2010pre, jacod2009microstructure, quasi-maximum likelihood estimator ait2010high, xiu2010quasi, local method of moments bibinger2014estimating, and robust pre-averaging realized volatility fan2018robust. These non-parametric estimation methods can estimate historical volatilities very well. However, in financial practices, we often need to predict future volatilities for the purpose of risk management and portfolio allocation, and these non-parametric methods cannot capture the market dynamics effectively for volatility prediction. Therefore, we shift our focus to parametric models that are known to be able to account for the financial market dynamics.

Given low-frequency data such as daily log returns, researchers introduced well-performing discrete-time models such as GARCH, stochastic volatility (SV), and vector autoregression (VAR) models to explain market dynamics by using historical volatilities and returns as innovations. These models are widely used for empirical financial analyses as they are easy to implement and can explain the market dynamics effectively. Thus, it is natural for researchers to harness the discrete-time model structures in the continuous-time frameworks. See engle2006multiple, hansen2012realized, kim2016unified, shephard2010realising for related research works. In their proposed methods, researchers used non-parametric realized volatility estimators from high-frequency data to make inferences for parametric discrete-time models at the low-frequency. Their empirical studies demonstrated that for a finite number of assets, by combining the low- and high-frequency methods, the proposed models perform better in volatility prediction tasks than traditional models with low- or high-frequency data alone. In financial applications, we often encounter a large number of assets, and models designed for the finite dimension become inconsistent and suffer from the curse of dimensionality. Thus, in this paper, we explore the approach to unify the discrete-time and continuous-time models appropriately under the high-dimensional set-up for the purpose of vast volatility matrix estimation and prediction.

To overcome the curse of dimensionality, it is common to impose sparsity on the entire vast volatility matrix bickel2008covariance, cai2011adaptive, kim2016asymptotic, tao2013optimal, wang2010vast. However, in finance, there exist common market factors such as industry sectors, inflation reports, Fed rate hikes, and oil prices, which affect the entire market. Therefore, assets are widely correlated, and the sparse condition imposed on the entire volatility matrix is not appropriate. To model this stylized feature, approximate factor models which indicate that the volatility matrix consists of the low-rank factor volatility and the sparse idiosyncratic volatility matrices were often used ait2017using, fan2018robust, fan2013large, fan2015incorporating, kim2019factor. Under this set-up, kim2019factor further described the eigenvalue process of the latent factor volatility matrix by a unified GARCH-It\^o model structure kim2016unified. In their proposal, the daily integrated eigenvalues of the latent factor volatility matrix process have historical squared factor returns as the innovations.

In light of the pioneer works, we model the high-frequency data from a large number of assets by the high-dimensional factor-based It\^o process that consists of the latent factor and the idiosyncratic diffusion processes, and further introduce the SV model structure to daily factor integrated volatility matrices that have the finite rank. In specific, we develop a continuous instantaneous factor volatility process which has an autoregressive (AR) structure at integer times so that it is a form of some interpolation of the AR structure. The instantaneous factor volatility process reflects the current market dynamics while its daily integrated volatility matrices retain the exact AR structure as in the SV model. We name our proposal the SV-It\^o model. When estimating the factor volatility matrix that is latent, we face an identifiability issue. To overcome this, researchers often impose some structures on the latent factor loading and factor volatility matrices such as the orthonormal and diagonal conditions used in kim2019factor. These conditions are restrictive in the sense that dynamics among market factors cannot be studied. In this paper, instead of imposing the strong diagonal condition on the factor volatility matrices, we assume stationary condition and follow the factor volatility matrix estimation procedure described in tao2011large. Based on the series of daily factor volatility matrix estimators and the imposed AR structure, we develop quasi-maximum likelihood estimation (QMLE) and least squares estimation (LSE) methods for model parameters whose asymptotic properties are established. With the proposed QMLE or LSE, as well as the imposed AR structure, we can estimate future factor volatility matrices effectively. On the other hand, the idiosyncratic volatilities are related to firm-specific risk so that the correlations among assets are weak. Thus, we impose the common sparse condition on the idiosyncratic volatility matrix and assume it to be time-invariant. When estimating the idiosyncratic volatility matrix, we harness its sparse structure and follow the principal orthogonal complement thresholding (POET) estimation procedure proposed by fan2013large, fan2015incorporating. Combining the future factor and idiosyncratic volatility matrix estimators, we propose an estimator for future vast volatility matrix and examine its asymptotic properties.

This paper is organized as follows. Section (ref) introduces a unified SV-It\^o model under the high-dimensional factor-based It\^o diffusion process and investigates its properties. Section (ref) develops non-parametric estimation methods for the latent factor loading and factor volatility matrices. We propose quasi-maximum likelihood estimation and least squares estimation procedures based on the low-rank factor volatilities. The asymptotic properties of proposed estimation methods are established. Section (ref) proposes an estimator for the conditional expected vast volatility matrix and investigates its asymptotic properties. Section (ref) presents numerical illustrations on the proposed methodologies and Section (ref) concludes the paper. Proofs are collected in Section (ref).

Unified discrete-time and continuous-time models

Discrete-time and continuous-time models

Both the discrete-time models such as GARCH and SV as well as the continuous-time models such as OU and CIR provide stochastic methods for financial data analyses. Discrete-time models are relatively simple parametric models and are often adopted to model the dynamic evolution of the volatility process based on the low-frequency data. Continuous-time models are described by more complicated stochastic differential equations instead and can provide non-parametric estimators for daily volatility based on the high-frequency data. These two types of models have unique characteristics and are not compatible. Since the low- and high-frequency data hold trading information for the same asset, it is natural to develop unified models and draw combined inferences. Some of the recent attempts include engle2006multiple, hansen2012realized, kim2016unified, kim2019factor, shephard2010realising, tao2011large. In this paper, we introduce unified discrete-time SV and continuous-time factor-based It\^o diffusion models by allowing both the number of low-frequency and high-frequency observations, as well as the number of assets, go to infinity.

Notation

For any given $p_1$-by-$p_2$ matrix $\bfm A = \left(A_{ij}\right)$, we denote its spectral norm by $\|\bfm U\|_2$, its Frobenius norm by $\|\bfm A\|_F = \sqrt{ \mathrm{Tr} (\bfm A^{\top} \bfm A) }$, and its max norm by $\| \bfm A \| _{\max} = \max_{i,j} | A_{ij}|$. Moreover, $\mathrm{vec} (\cdot)$ denotes the operator that stacks the columns of a matrix, $\mathrm{vec}^{-1} (\cdot)$ is the inverse operator of $\mathrm{vec}(\cdot)$, and $\mathrm{vech}(\cdot )$ is a column vector obtained by vectorizing only the lower triangular part of a matrix. $\mathrm{Diag} (\cdot)$ returns a square diagonal matrix with the elements of a vector on the main diagonal, and $\mathrm{diag} (\cdot)$ returns a column vector of the main diagonal elements of a matrix. Also let $\mathrm{det}(\bfm A)$ be the determinant of a matrix $\bfm A$ and $ \mathrm{Tr}(\bfm A)$ be the trace of $\bfm A$. Let $C$ be a generic constant whose values are free of $n, m,$ and $p$, and may change from appearance to appearance.

Unified models

Let $\bfm X _t = \(X_{1,t} , \ldots, X_{p,t} \)^{\top}$ be the vector of true underlying log prices of $p$ assets at time $t$. In finance, we usually assume that high-frequency data $\bfm X_t$ obey some continuous diffusion process. To account for common market factors in financial industry, we consider the following factor-based diffusion process:

equation[equation omitted — 115 chars of source]

where $\ensuremath{\boldsymbol{\mu}}_t \in \mathbb{R}^{p}$ is a drift vector, $\bfm L$ is a $p$-by-$r$ factor loading matrix and $r$ is the total number of market factors. Moreover, $\bfm f_t$ and $\bfm u_t$ are the factor and idiosyncratic diffusion processes, respectively, and obey the following

equation*[equation* omitted — 164 chars of source]

where $\bfsym \sigma_t$ is an $r$-by-$r$ matrix, $\ensuremath{\boldsymbol{\vartheta}}_t$ is a $p$-by-$p$ matrix, $\bfm B_t$ and $\bfm W_t$ are independent $r$-dimensional and $p$-dimensional Brownian motions, respectively. The stochastic processes $\ensuremath{\boldsymbol{\mu}}_t$, $\bfsym \sigma_t$, and $\ensuremath{\boldsymbol{\vartheta}}_t$ are defined on a filtered probability space $\( \Omega, {\cal F}, \{{\cal F}_t, t\in [0, \infty)\}, P \)$ with filtration ${\cal F}_t$ satisfying the usual conditions. The daily integrated volatility matrices follow to be

equation*[equation* omitted — 152 chars of source]

where $\ensuremath{\boldsymbol{\Psi}}_k = \int_{k-1}^k \bfsym \sigma_t ^{\top} \bfsym \sigma_t dt $ and $\bfsym \Gamma_k ^s = \int_{k-1} ^k \ensuremath{\boldsymbol{\vartheta}}_t ^{\top} \ensuremath{\boldsymbol{\vartheta}}_t dt$.

We note that the idiosyncratic process $\bfm u_t$ corresponds to the firm-specific risk so that its volatility matrix may be sparse. Moreover, since firm-specific risk is generally unpredictable, investors seek to minimize its negative impact on an investment portfolio by diversification or hedging. On the other hand, the latent factor process $\bfm f_t$ corresponds to the systematic risk or undiversifiable risk that affects the whole market. Thus, it is natural to capture market dynamics by modeling the latent factor process. In light of this, we propose the following latent factor diffusion process $\bfm f_t$ that embeds the SV model. Define the instantaneous volatility process of $\bfm f_t$ by

equation*[equation* omitted — 77 chars of source]

and

eqnarray[eqnarray omitted — 394 chars of source]

where $\bfm Z_t = (Z_{i,t}) _{i=1,\ldots,r} = \int_{[t]} ^{t} \ensuremath{\boldsymbol{\nu}} ^{\top} d \bfm B_s ^1$, $\bfm B_t^1$ is an $r$-dimensional standard Brownian motion, $[t]$ denotes the integer part of $t$ except that $[t] = t -1$ when $t$ is an integer, and $\bfsym \alpha_j$'s are $r$-by-$r$ matrices. We name our proposal the SV-It\^{o} model.

The instantaneous volatility process $\bfsym \Sigma_t$ is almost surely continuous with respect to time $t$. When restricted the integer time points, $\bfsym \Sigma_t$ of the SV-It\^o model retains the following AR structure,

equation*[equation* omitted — 210 chars of source]

Thus, the instantaneous volatility process is formed by some interpolation of the AR structure where the random fluctuation is explained by the $\bfm Z_t$ process. Also, current market dynamics are reflected through the terms $\int_{[t]} ^{t} \bfsym \Sigma_s ds$ and $\bfm Z_t$.

For statistical inferences, we have an interest in developing parametric models based on the series of daily integrated factor volatility matrices, i.e., $\int_{k-1}^k \bfsym \Sigma_t dt$, $k=1,\ldots,n$. In the following proposition, we show that the daily integrated factor volatility matrices have the discrete-time AR structure.

propositionUnder the SV-It\^o model, we have the following iterative relations. \begin{enumerate} • For any $n, k \in \mathbb{N}$, we have \begin{eqnarray*} &&\bfm F(k) \equiv \int_{n-1}^n \frac{(n-t)^k}{k!} \mathrm{vec}( \bfsym \Sigma_t ) dt \cr &&= \frac{ \mathrm{vec} (\bfsym \alpha_0 \bfsym \alpha_0 ^{\top}) + (k+1) \mathrm{vec}( \bfsym \Sigma _{n-1} ) + \sum_{j=1}^{q-1} \bfm A_{j+1} \mathrm{vec}(\ensuremath{\boldsymbol{\Psi}} _{n-j } ) }{(k+2)!} + \frac{k+1}{(k+3)!} \mathrm{vec}(\ensuremath{\boldsymbol{\nu}} ^{\top} \ensuremath{\boldsymbol{\nu}} ) \cr && + \mathrm{vec} \( \( \int_{n-1}^n \frac{(k+1) (n-t) ^{k+2}}{ (k+2)!} Z_{j,t} d Z_{i,t} + \int_{n-1}^n \frac{(k+1) (n-t) ^{k+2}}{ (k+2)!} Z_{i,t} d Z_{j,t} \) _{i,j=1,\ldots, r} \) \cr && + \bfm A_{1} \bfm F(k+1), \end{eqnarray*} where $\bfm A_j =\bfsym \alpha_{j} \otimes \bfsym \alpha_{j} $ for $j=1,\ldots,q$, and the operator $\otimes$ denotes the Kronecker product. • For $\mathrm{det}(\bfsym \alpha_1) \neq 0$ and $\|\bfsym \alpha_1\|_2 <1$, we have \begin{eqnarray*} \int_{n-1}^n \mathrm{vec}( \bfsym \Sigma_t ) dt &=& \ensuremath{\boldsymbol{\varrho}}_1 \mathrm{vec} (\bfsym \alpha_0 \bfsym \alpha_0 ^{\top} ) + \( \ensuremath{\boldsymbol{\varrho}}_2 - 2 \ensuremath{\boldsymbol{\varrho}}_3 \) \mathrm{vec}(\ensuremath{\boldsymbol{\nu}} ^{\top} \ensuremath{\boldsymbol{\nu}} ) \cr &&+ \sum_{j=1}^{q-1} \{ (\ensuremath{\boldsymbol{\varrho}}_1- \ensuremath{\boldsymbol{\varrho}}_2) \bfm A_{j } + \ensuremath{\boldsymbol{\varrho}}_2 \bfm A_{j+1} \} \mathrm{vec}(\ensuremath{\boldsymbol{\Psi}} _{n-j} ) \cr &&+ (\ensuremath{\boldsymbol{\varrho}}_1- \ensuremath{\boldsymbol{\varrho}}_2) \bfm A_{q } \mathrm{vec}(\ensuremath{\boldsymbol{\Psi}} _{n-q} ) + \bfm D_n a.s., \end{eqnarray*} where $\ensuremath{\boldsymbol{\varrho}}_1= \bfm A_1 ^{-1} \( e ^{\bfm A_1} - \bfm I_{r^2} \)$, $\ensuremath{\boldsymbol{\varrho}}_2= \bfm A_1 ^{-2} \( e ^{\bfm A_1} - \bfm I_{r^2} - \bfm A_1 \)$, $\ensuremath{\boldsymbol{\varrho}}_3= \bfm A_1 ^{-3} ( e ^{\bfm A_1} - \bfm I_{r^2} - \bfm A_1 - \frac{\bfm A_1 ^2}{2})$, $\bfm I_{r^2}$ is the $r^2$-dimensional identity matrix, $e^{\bfm A} = \sum_{k=0}^\infty \frac{\bfm A^k} {k!}$, and \begin{eqnarray*} \bfm D_n &=& \sum_{k=0}^{\infty} \bfm A_1 ^k \mathrm{vec} \Big ( \Big ( \int_{n-1}^n \frac{(k+1) (n-t) ^{k+2}}{ (k+2)!} Z_{j,t} d Z_{i,t} \cr && \qquad \qquad \qquad \qquad \qquad + \int_{n-1}^n \frac{(k+1) (n-t) ^{k+2}}{ (k+2)!} Z_{i,t} d Z_{j,t} \Big ) _{i,j=1,\ldots, r} \Big ). \end{eqnarray*} \end{enumerate}

Proposition (ref) shows that the daily integrated volatility matrices retain some AR structure. Moreover, by the construction, there exist $r(r+1)/2$ vector $\bfsym \beta_0$ and $r(r+1)/2$-by-$r(r+1)/2$ matrices $\bfsym \beta_j, j=1,\ldots, q$, such that

equation[equation omitted — 206 chars of source]

where $ \bfsym \beta_0 =\mathrm{vech} \( \mathrm{vec} ^{-1} \( \ensuremath{\boldsymbol{\varrho}}_1 \mathrm{vec} (\bfsym \alpha_0 \bfsym \alpha_0 ^{\top} ) + \( \ensuremath{\boldsymbol{\varrho}}_2 - 2 \ensuremath{\boldsymbol{\varrho}}_3 \) \mathrm{vec}(\ensuremath{\boldsymbol{\nu}} ^{\top} \ensuremath{\boldsymbol{\nu}} ) \) \)$, and the $i$th row of $\bfsym \beta_j$ is obtained by $ \mathrm{vech} ( \mathrm{vec} ^{-1} \( \bfm A \) + \mathrm{vec} ^{-1} \( \bfm A ^{\top} \) - \mathrm{Diag} \( \mathrm{diag} \( \mathrm{vec} ^{-1} \( \bfm A \) \) \) )$ for the $i$th row of the coefficient matrix $\bfm A$ corresponding to $\mathrm{vec} (\ensuremath{\boldsymbol{\Psi}} _{n-j})$, also $\bfm e_n = \mathrm{vech} \( \mathrm{vec}^{-1} ( \bfm D_n) \)$ a.s. Given the SV-It\^o model, volatility dynamics can be explained by the previous volatilities $\ensuremath{\boldsymbol{\Psi}}_k$'s and random fluctuation is modeled by the martingale difference term $\bfm e_n$ which comes from the $\bfm Z_t$ process in (ref). Furthermore, under some stationary conditions, the expectation of the daily integrated factor volatility is

equation*[equation* omitted — 189 chars of source]

where the largest eigenvalue of $\sum_{j=1}^q \bfsym \beta_j $ should be strictly less than 1.

remarkRecently, kim2019factor introduced a factor GARCH-It\^o model to capture the market dynamics with high-frequency data and for a large number of assets. In specific, they imposed some GARCH-type dynamic structure on the eigenvalue process of the factor volatility matrix process, so the daily integrated eigenvalues are functions of historical squared factor returns. The proposed SV-It\^o model in this paper also imposes some dynamic AR structure on the factor volatility process so that the two models share similarities in their approaches. However, to explain the market dynamics at the low-frequency, the factor GARCH-It\^o model employs low-frequency factor return information while the SV-It\^o model uses the integrated factor volatilities over low-frequency periods. Our empirical study shows that incorporating integrated factor volatilities helps to capture the market dynamics promptly (see Section (ref)). Moreover, the factor GARCH-It\^o model is restricted to the eigenvalue structure so that it cannot capture the cross-sectional dynamics. On the other hand, the SV-It\^o model adopts a more general structure so that the dynamics for correlations among market factors can be modeled as well (see Section (ref)).

Parameter estimation

In this section, we introduce estimation procedures for the factor loading and factor volatility matrices, as well as for the model parameteres, and establish their asymptotic properties.

The model set-up and realized volatility matrix estimators

Let $p$ be the total number of assets and $n$ be the total number of low-frequency observations. The high-frequency prices for the $i$th asset during the $k$th low-frequency period are observed at times $t_{i,k, \ell } \in (k-1, k]$, $i=1, \ldots, p$, $k=1, \ldots,n$ and $\ell=1, \ldots, m_{i,k}$. Denote $m_i$ the averaged number of high-frequency observations during each low-frequency period for the $i$th asset, that is, $m_i=\sum_{k=1}^n m_{i,k}/n$. Further let $m=\sum_{i=1}^p m_i/p$. Let $Y_{i, t_{i,k,\ell}}$ be the observed log price of the $i$th asset at time $t_{i,k,\ell}$. High-frequency data are normally non-synchronized so that $t_{i_1,k, \ell} \neq t_{i_2, k, \ell}$ for $i_1 \neq i_2$. Moreover, due to imperfections of the trading mechanisms ait2008high, high-frequency data are often contaminated by market microstructure noises so that the observed log price $Y_{i,t_{i,k,\ell}}$ is a noisy version of the corresponding true log price $X_{i, t_{i,k,\ell}}$, that is,

equation[equation omitted — 184 chars of source]

where $\varepsilon_{i, t_{i,k,\ell}}$'s are stationary noises with mean zero and variance $\eta_i$. We further assume that $\varepsilon_i$ and $X_i$ are independent with each other.

Given non-synchronized and noisy high-frequency data, researchers constructed nonparametric realized volatility matrix estimators that take advantage of sub-sampling and local-averaging techniques to remove the effect of market microstructure noises so that the integrated volatility matrix $\bfsym \Gamma_k$ can be estimated consistently and efficiently. Examples include multi-scale realized volatility matrix (MSRVM) zhang2011estimating, pre-averaging realized volatility matrix (PRVM) christensen2010pre, and kernel realized volatility matrix (KRVM) barndorff2011multivariate estimators. See also ait2010high, bibinger2014estimating, fan2018robust, kim2018large, wang2010vast. When the number of assets is $p$ finite, these estimators can achieve the optimal convergence rate of $m^{-1/4}$ in the presence of market microstructure noises for estimating $\bfsym \Gamma_{k}$.

Non-parametric factor volatility matrix estimation

When the number of assets $p$ is finite, we may view the factor process $\bfm f_t$ as the log price process $\bfm X_t$ and estimate the daily integrated factor volatility matrices $\ensuremath{\boldsymbol{\Psi}}_k$'s with the well-performing estimators such as the MSRVM, the PRVM, and the KRVM. The estimators of $\ensuremath{\boldsymbol{\Psi}}_k$'s and the AR structure described in (ref) can be employed to estimate model parameter $\bfsym \beta_j$'s directly in this case (see Sections (ref) and (ref)). However, in practice, we often encounter a large number of assets and the factor process $\bfm f_t$ is latent. In this section, we first discuss how to estimate the latent factor volatility matrices $\ensuremath{\boldsymbol{\Psi}}_k$'s in the high-dimensional set-up.

Given the factor-based It\^o diffusion process in (ref), the daily integrated volatility is

equation*[equation* omitted — 156 chars of source]

The idiosyncratic risk is related to the firm-specific risk so that the corresponding co-volatility matrices are sparse. In light of this, we impose the following sparse condition on the idiosyncratic volatility matrix $\bfsym \Gamma_{k}^s = \(\Gamma_{k, ij} ^s \) _{i,j=1,\ldots, p}$:

equation[equation omitted — 211 chars of source]

where $\delta \in [0,1)$, $M$ is a positive bounded random variable and the sparsity level $\pi(p)$ diverges very slowly such as $\log p$. On the other hand, the factor process $\bfm f_t$ often depends on a few common market factors such as industry sectors, inflation reports, Fed rate hikes, and oil prices. Thus, the number of market factors $r$ takes a much smaller value than the number of assets $p$, so we assume that the rank $r$ is finite. In this paper, we further assume that $r$ is known. For the latent factor model, we face the identifiability issue. To manage this issue, researchers often impose some structures on the factor loading matrix $\bfm L$ such as $\bfm L ^{\top} \bfm L=p \bfm I _r$ and also assume that the factor volatility matrices $\ensuremath{\boldsymbol{\Psi}}_k$'s are diagonal. It follows that $\bfm L$ and $\ensuremath{\boldsymbol{\Psi}}_k$'s are corresponding to the eigenvectors and eigenvalues of the factor volatility matrices. Under these assumptions, kim2019factor proposed the factor GARCH-It\^o model for the eigenvalues of the factor volatility matrices. However, in this case, the correlation structure among assets is constant, which makes it difficult to investigate the cross-sectional market dynamics.

To account for the dynamics among assets, we consider the following structure that is more general. First, we assume that the idiosyncratic volatility matrix $\bfsym \Gamma_k^s$ satisfies

equation[equation omitted — 136 chars of source]

Moreover, we assume that the factor volatility matrices $\ensuremath{\boldsymbol{\Psi}}_k$'s form a stationary process such that

equation*[equation* omitted — 218 chars of source]

as $n$ goes to infinity, where the liming variable $ \ensuremath{\boldsymbol{\Psi}}_{\infty}^2 $ may be $\operatorname{E} \[ \left \{ \ensuremath{\boldsymbol{\Psi}}_1 -\operatorname{E} \( \ensuremath{\boldsymbol{\Psi}}_1\)\right \} ^2 \]$. For the factor loading matrix $\bfm L$, we assume that it satisfies $\bfm L ^{\top} \bfm L = p \bfm I_r$ and $\bfm L \bfm V = \bfm L$, where $\bfm V$ denotes the eigenvectors of $\ensuremath{\boldsymbol{\Psi}}_{\infty}^2$.

Under these conditions, we propose the following procedure to estimate the factor loading and volatility matrices. First, we consider

eqnarray*[eqnarray* omitted — 422 chars of source]

where $\bar{\bfsym \Gamma}=\frac{1}{n} \sum_{k=1}^n \bfsym \Gamma_{k}$ and $\bar{\ensuremath{\boldsymbol{\Psi}}} =\frac{1}{n} \sum_{k=1}^n \ensuremath{\boldsymbol{\Psi}}_{k}$. Note that $\mathbb{S}_n$ is the rank $r$ matrix and is free of the idiosyncratic volatility matrix $\bfsym \Gamma^s$. As $n \to \infty$, note that we have

equation*[equation* omitted — 124 chars of source]

Thus, the scaled factor loading matrix $p^{-1/2}\bfm L$ relates to the eigenvectors of $\mathbb{S}_{\infty}= \bfm L \ensuremath{\boldsymbol{\Psi}}_{\infty}^2 \bfm L^{\top}$. We estimate $p^{-1/2}\bfm L$ by the first $r$ eigenvectors of $\mathbb{S}_{n}$. However, in practice, $\mathbb{S}_n$ is not observable so that we use its estimator,

equation[equation omitted — 168 chars of source]

where $\bar{\widehat{\bfsym \Gamma}} = \frac{1}{n} \sum_{k=1}^n \widehat{\bfsym \Gamma}_k$ and $\widehat{\bfsym \Gamma}_{k}$ is the MSRVM, the PRVM, or the KRVM estimator for daily integrated volatility matrix $\bfsym \Gamma_{k}$. We then estimate the factor loading matrix $p^{-1/2} \bfm L$ by the first $r$ eigenvectors of $\widehat{\mathbb{S}}_{n,m}$ and denote the estimator as $p^{-1/2}\widehat{\bfm L}$. Finally, the factor volatility matrix estimators can be obtained by

equation*[equation* omitted — 164 chars of source]
remarkIn the literature on approximate factor models, researchers often assume that the factor volatility matrices $\ensuremath{\boldsymbol{\Psi}}_k$'s are diagonal so that their entries are related to the eigenvalues. Our model structure includes the diagonal structure and when it is assumed, we do not need to require the stationary condition.

To investigate the asymptotic behaviors of the proposed non-parametric estimators, we need the following technical conditions.

assumption\begin{enumerate} • For some given $b\geq 2$, $ \sup_{k \in \mathbb{N}} \max_{1\leq i,j \leq p} \operatorname{E} \( |\widehat{\Gamma}_{k, ij} -\Gamma_{k,ij} | ^{2b} \) \leq C m^{-b/2}$; • $ \operatorname{E} \( \|\frac{1}{n} \sum_{k=1}^n \ensuremath{\boldsymbol{\Psi}}_k ^2 - \( \frac{1}{n} \sum_{k=1}^n \ensuremath{\boldsymbol{\Psi}}_k \)^2 - \ensuremath{\boldsymbol{\Psi}}_{\infty}^2 \| _{F} ^b\) \leq C n^{-b/2}$; • For all $j=1,\ldots, r$, $\lambda_{j}( \mathbb{S}_{\infty})- \lambda_{j+1} (\mathbb{S}_{\infty}) \geq C p$ for some fixed constant $C$, where $\lambda_{j} (\bfm A)$ is the $j$th largest eigenvalue of the square matrix $\bfm A$. \end{enumerate}
remarkAssumption (ref) (a) is satisfied when the instantaneous volatility processes and noise have the finite $4b$th moment kim2016asymptotic, tao2013fast. Assumption (ref) (b) is required to obtain the consistent estimator for $\bfm L$ that is uniquely defined. Finally, Assumption (ref) (c) is the so-called pervasive condition which is often imposed when investigating the approximate factor models fan2013large.

We present the theorem that investigates the asymptotic behaviors of the proposed non-parametric estimators for the factor loading and factor integrated volatility matrix.

thmUnder the models (ref) and (ref), when Assumption (ref), (ref), and the sparsity condition (ref) are met, we have \begin{eqnarray} && \operatorname{E} \( \| \widehat{\mathbb{S}}_{n,m} -\mathbb{S}_{\infty} \|_F ^b \) \leq C p^b \( m^{-b/4} + n^{-b/2} \), \\ && \max_{1\leq i\leq r} \operatorname{E} \( \| p^{-1/2} \widehat{\bfm L}_i - p^{-1/2} sign (\widehat{\bfm L}_i ^{\top} \bfm L_i) \bfm L_i \|_F ^b \) \leq C \( m^{-b/4} + n^{-b/2} \), \\ && \sup_{k \in \mathbb{N}} \operatorname{E} \( \| \widehat{\ensuremath{\boldsymbol{\Psi}}}_{k} -\ensuremath{\boldsymbol{\Psi}}_{k} \|_F ^{2b/3} \) \leq C \left \{ n^{-b/3} + m^{-b/6} + \( \pi(p) /p \) ^{2b/3} \right \}, \end{eqnarray} where $\bfm L_i$ is the $i$th column of $\bfm L$.
remarkTheorem (ref) shows that the latent factor process can be estimated consistently with the convergence rate $n^{-1/2} + m^{-1/4}+\pi(p)/p$. The term $n^{-1/2}$ is coming from identifying the latent factor loading matrix $\bfm L$ under the stationary condition. However, if we do impose the diagonal structure on the factor volatility matrix $\ensuremath{\boldsymbol{\Psi}}_k$, the term $n^{-1/2}$ is removed. The term $m^{-1/4}$ is coming from estimating the daily integrated volatility matrix $\bfsym \Gamma_k$, which is known as the optimal convergence rate in the presence of the market microstructure noises. Finally, the term $\pi(p)/p$ is the cost to identify the latent factor volatility matrix $\ensuremath{\boldsymbol{\Psi}}_k$ from the integrated volatility matrix $\bfsym \Gamma_k$. These results imply Assumption (ref) (d) and helps to establish the convergence rate in Theorem (ref).

Quasi-maximum likelihood estimation

In this section, we propose a quasi-maximum likelihood estimation procedure for the true parameter $\ensuremath{\boldsymbol{\theta}}_0 =(\bfsym \beta_{0,0}^{\top}, \mathrm{vec}(\bfsym \beta_{0,1})^{\top},\ldots, \mathrm{vec} (\bfsym \beta_{0,q})^{\top})^{\top} \in \mathbb{R} ^d$, where $\bfsym \beta_j$'s are defined in (ref) and $d=r (r +1)\{ 2+q r (r+1) \}/4$.

We first develop the estimation procedure by pretending that $\ensuremath{\boldsymbol{\Psi}}_k$'s are known and consider the following quasi-likelihood function:

equation*[equation* omitted — 264 chars of source]

where $\mathrm{vech} \( \bfm H_k (\ensuremath{\boldsymbol{\theta}}) \) = \bfsym \beta_0 + \sum_{j=1}^q \bfsym \beta_j \mathrm{vech}(\ensuremath{\boldsymbol{\Psi}}_{k-j})$ and $\bfm H_k (\ensuremath{\boldsymbol{\theta}}) $ is symmetric. The difference between $\ensuremath{\boldsymbol{\Psi}}_k$ and $\bfm H_k (\ensuremath{\boldsymbol{\theta}}_0)$ under the SV-It\^o model is the martingale difference whose vectorization is $\bfm e_n$ defined in (ref). Then, under some technical conditions, we can show that the maximizer of $L_{n} (\ensuremath{\boldsymbol{\theta}})$ converges to the true parameter $\ensuremath{\boldsymbol{\theta}}_0$ with the convergence rate of $n^{-1/2}$. However, the daily integrated factor volatility matrix $\ensuremath{\boldsymbol{\Psi}}_k$'s are unobservable so that we need to follow the procedure developed in Section (ref) to obtain their estimators $\widehat{\ensuremath{\boldsymbol{\Psi}}}_k$'s. Given the estimators $\widehat{\ensuremath{\boldsymbol{\Psi}}}_k$'s, we let

equation[equation omitted — 219 chars of source]

where $\widehat{\bfm H}_k (\ensuremath{\boldsymbol{\theta}})$ is symmetric, and define the following quasi-likelihood function for parameter estimation:

equation[equation omitted — 341 chars of source]

The true model parameters $\ensuremath{\boldsymbol{\theta}}_0$ can be obtained by maximizing the quasi-likelihood function in (ref), that is,

equation*[equation* omitted — 200 chars of source]

where $\ensuremath{\boldsymbol{\Theta}}$ is the parametric space of $\ensuremath{\boldsymbol{\theta}}$.

To investigate the asymptotic behaviors of the proposed QMLE, we need the following technical conditions.

assumption\begin{enumerate} • $\ensuremath{\boldsymbol{\Theta}}$ is compact; $\ensuremath{\boldsymbol{\theta}}_0$ is an interior point of $\ensuremath{\boldsymbol{\Theta}}$; • $\min_{k \in \mathbb{N}} \inf_{\ensuremath{\boldsymbol{\theta}} \in \ensuremath{\boldsymbol{\Theta}}} \lambda_{\min} (\bfm H_k (\ensuremath{\boldsymbol{\theta}})) > c_\lambda $ and $\min_{k \in \mathbb{N}} \inf_{\ensuremath{\boldsymbol{\theta}} \in \ensuremath{\boldsymbol{\Theta}}} \lambda_{\min} (\widehat{\bfm H}_k (\ensuremath{\boldsymbol{\theta}})) > c_\lambda$ a.s. for some fixed positive constant $c_{\lambda}$, $\lambda_{\min} (\bfm A)$ is the smallest eigenvalue of a square matrix $\bfm A$; • There exists some fixed constants $C_1$ and $C_2$ such that $C_1 m_i \leq m_{i,k} \leq C_2 m_i$ for all $i$ and $k$, and $C_1 m \leq m_i \leq C_2 m$ for all $i$. • There is some fixed sequence $\tau_m$ such that \begin{equation*} \sup_{k \in \mathbb{N}} \operatorname{E} \( \| \widehat{\ensuremath{\boldsymbol{\Psi}}}_k -\ensuremath{\boldsymbol{\Psi}}_k \|_F^4 \) \leq \tau_m ^{4} =o(1). \end{equation*} \end{enumerate}
remarkAssumption (ref) (a) is required to define the parameter uniquely. Assumption (ref) (b) is often obtained by imposing some positive definitiveness requirement on the intercept part $\bfsym \beta_0$. Assumption (ref) (d) is required to investigate the asymptotic behavior of QMLE $\widehat{\ensuremath{\boldsymbol{\theta}}}$ since the latent factor volatility matrices $\ensuremath{\boldsymbol{\Psi}}_k$'s are not observable, and we need to estimate them. When the number of assets $p$ is fixed, a simple way to establish the rate $\tau_m$ is to treat the factor process $\bfm f_t$ as the log price process $\bfm X_t$, then the integrated volatility matrices $\ensuremath{\boldsymbol{\Psi}}_k$'s can be estimated by the MSRVM, the PRVM, or the KRVM estimator. These well-performing estimators can achieve the optimal convergence rate of $m^{-1/4}$ in estimating $\ensuremath{\boldsymbol{\Psi}}_k$'s when observed stock prices are contaminated by market microstructure noises. $\tau_m$ follows to be $m^{-1/4}$ in this case. On the other hand, in the high-dimensional set-up, we need to identify the latent factor volatility matrix and are required to impose some structures on the factor loading matrix so that $\tau_m$ can be appropriately established. In this case, the estimation procedure presented in Section (ref) has $\tau_m=n^{-1/2} + m^{-1/4}+\pi(p)/p$ (see Theorem (ref)).

The following theorem provides the convergence rate of the QMLE $\widehat{\ensuremath{\boldsymbol{\theta}}}$.

thmUnder Assumption (ref), we have \begin{equation} \left \| \widehat{\ensuremath{\boldsymbol{\theta}}} -\ensuremath{\boldsymbol{\theta}}_0 \right \|_{\max} = O_p (\tau_m + n^{-1/2}). \end{equation}
remarkTheorem (ref) shows that the convergence rate of the QMLE $\widehat{\ensuremath{\boldsymbol{\theta}}}$ is $\tau_m + n^{-1/2}$. The first term $\tau_m$ is the cost to estimate the integrated factor volatility matrices $\ensuremath{\boldsymbol{\Psi}}_k$'s and takes the value $m^{-1/4}+ n^{-1/2} + p/\pi(p)$ (see Section (ref)). The second term $n^{-1/2}$ is the usual convergence rate for parametric estimation based on the low-frequency structure. The asymptotic result obtained in Theorem (ref) holds as long as the daily integrated factor volatility matrices $\ensuremath{\boldsymbol{\Psi}}_k$'s satisfy (ref). That is, the asymptotic result does not depend on the specific instantaneous volatility process described in (ref). For example, we can develop an instantaneous volatility process that can capture intraday volatility dynamics such as the U-shape pattern admati1988theory, andersen1997intraday, andersen2018time, hong2000trading in the following: \begin{eqnarray*} \bfsym \Sigma_t &=& \bfsym \Sigma _{[t]} + (t-[t])^2 \bfsym \alpha_0 ^{\prime}\( \bfm I_r + \sum_{j=1}^{q-1} \bfsym \alpha_{j+1} \ensuremath{\boldsymbol{\Psi}} _{[t]-j+1} \bfsym \alpha_{j+1} ^{ \top} \) \bfsym \alpha_0 ^{\prime \top} \cr && - (t-[t]) \( \bfsym \alpha_0 \bfsym \alpha_0 ^{\top} + \bfsym \Sigma _{[t]} + \sum_{j=1}^{q-1} \bfsym \alpha_{j+1} \ensuremath{\boldsymbol{\Psi}} _{[t]-j+1} \bfsym \alpha_{j+1} ^{\top} \) \cr &&+ \bfsym \alpha_{1} \(\int_{[t]} ^{t} \bfsym \Sigma_s ds \) \bfsym \alpha_1 ^{\top} + ([t] +1 -t) \bfm Z _t \bfm Z_t^{\top}. \end{eqnarray*} We then still obtain the AR structure as described in (ref). If consistent estimators for the instantaneous volatility $\bfsym \Sigma_t$'s were available, we could also study the intraday dynamics. However, this is not the focus of this paper so that we leave it for future research.

Least squares estimation

When considering the vector autoregression form (ref) directly, one of the natural ways to estimate the parameters is the well-known least squares estimation method. More specifically, we define the square loss function as:

equation*[equation* omitted — 289 chars of source]

The true model parameters $\ensuremath{\boldsymbol{\theta}}_0$ can be obtained by minimizing the above square loss function, that is,

equation*[equation* omitted — 210 chars of source]

We call the estimator the LSE. Its asymptotic behavior can be showed similar to the proofs of Theorem (ref), and the convergence rate is

equation*[equation* omitted — 157 chars of source]

It is easy to implement the LSE method as it has the closed form which is the well-known ordinary least squares estimator. However, the vector autoregression equation (ref) results in heterogeneous noise described by the term $\bfm e_n$, which may cause inefficiency in estimating model parameters using the LSE method. On the other hand, the QMLE method adjusts the heterogeneous variance term using the inverse of the conditional variance matrix function $\widehat{\bfm H}^{-1} (\ensuremath{\boldsymbol{\theta}})$. However, the implementation of the QMLE is demanding and its performance may depend on the choice of the initial value of the optimization. To overcome these issues, we suggest adopting the LSE as the initial value for the QMLE procedure.

Large volatility matrix prediction

In finance applications such as portfolio allocation, we are often required to predict future large volatility matrix. In this section, we demonstrate how to harness the proposed SV-It\^o model for constructing a predictor. It is known that the best predictor for future large volatility matrix is the conditional expected value of the integrated volatility matrix given current information. Under the SV-It\^o model and the condition (ref), we have

equation*[equation* omitted — 198 chars of source]

That is, the conditional expectation consists of the conditional expected factor volatility and idiosyncratic volatility matrices.

We first discuss the estimation of the idiosyncratic volatility matrix. Note that

eqnarray*[eqnarray* omitted — 191 chars of source]

where the idiosyncratic volatility matrix $\bfsym \Gamma^s$ satisfies the sparsity condition (ref), and the averaged factor volatility matrix $\bfsym \Phi_n=\bfm L \(\frac{1}{n} \sum_{k=1}^{n} \ensuremath{\boldsymbol{\Psi}} _k \) \bfm L^{\top}$ has the finite rank $r$. Therefore, the averaged integrated volatility matrix $\bar{\bfsym \Gamma}$ also retains a low-rank plus sparse structure. Given this structure, we can apply the POET procedure introduced by fan2013large to estimate the idiosyncratic volatility matrix $\bfsym \Gamma^s$. More specifically, the input of idiosyncratic volatility estimator is

equation*[equation* omitted — 232 chars of source]

where $\widehat{\lambda}_j$ is the $j$th largest eigenvalue of $\bar{\widehat{\bfsym \Gamma}}$ and $\widehat{\bfm q}_j$ is the corresponding eigenvector. We then apply the adaptive threshold scheme to the input of idiosyncratic volatility matrix estimator as follows:

eqnarray*[eqnarray* omitted — 348 chars of source]

where the thresholding function $s_{ij} (\widetilde{\Gamma}_{ij}^s)$ satisfies $|s_{ij} (\widetilde{\Gamma}_{ij}^s) - \widetilde{\Gamma}_{ij} ^s | \leq \varpi_{ij}$, and we use the adaptive threshold level $\varpi_{ij}= \varpi _m \sqrt{(\widetilde{\Gamma}_{ii} ^s\vee 0 )(\widetilde{\Gamma}_{jj}^s \vee 0 )} $ which is the same as applying the threshold $\varpi_m$ to the correlation.

On the other hand, with the QMLE $\widehat{\ensuremath{\boldsymbol{\theta}}}$, we estimate the factor volatility matrix by

equation*[equation* omitted — 129 chars of source]

where $\widehat{\bfm L}$ is estimated by the first $r$ eigenvectors of $\widehat{\mathbb{S}}_{n,m}$ defined in (ref) and the AR structure for daily integrated volatility matrices $\widehat{\bfm H}_k (\ensuremath{\boldsymbol{\theta}})$ is defined in (ref). Combining the idiosyncratic volatility matrix estimator $\widehat{\bfsym \Gamma}^s$, we estimate the future large volatility matrix by

equation*[equation* omitted — 191 chars of source]

We name the proposed estimator the SV-It\^o POET (SV-POET) estimator.

To investigate the asymptotic behaviors of the proposed estimator, we need the following assumptions.

assumption\begin{enumerate} • For some fixed positive constant $C_1$, we have \begin{equation*} \frac{p}{r} \max_{1\leq i \leq p} \sum_{j=1}^r q_{ij} ^2 \leq C_1 a.s., \end{equation*} where $\bfm q_j = (q_{1j}, \ldots, q_{pj})^{\top}$ is the eigenvector of the averaged factor volatility matrix $\bfsym \Phi_n$ corresponding to the $j$th largest eigenvalue; • We have $D_{\lambda} \geq C_2 p$ and $\lambda_1 / D_\lambda \leq C_3 $ a.s., where $D_\lambda= \min \{\bar{ \lambda}_i -\bar{\lambda}_{i+1}: i=1,\ldots,r \}$, $\bar{\lambda}_i$ is the $i$th largest eigenvalue of $\bfsym \Phi_n$, $\lambda_1$ is the largest eigenvalue of $\bar{\bfsym \Gamma}$, and the smallest eigenvalue of $\bfsym \Gamma^s$ stays away from zero; • $\pi(p) /p^{1/2} + \sqrt{ \log p / (n m^{1/2}+m) }= o(1)$. \end{enumerate}
remarkAssumption (ref) (a) and (b) are called the incoherence condition and pervasive condition, respectively, which are often imposed in analyzing low-rank matrix and approximate factor models ait2017using, candes2011robust, fan2013large.

The following theorem investigates the asymptotic behaviors of the SV-POET.

thmUnder the models (ref), (ref), and (ref), the following concentration inequality, \begin{equation} \Pr \left \{\max_{1 \leq i,j \leq p} | \bar{\widehat{\Gamma}}_{ij} -\bar{\Gamma}_{ij} | \geq C \sqrt{\frac{\log p}{m^{1/2} n +m}} \right \} \leq p^{-1}, \end{equation} Assumptions (ref)--(ref), (ref), and the sparsity condition (ref) are met. Take $\varpi_m = C_{\varpi} ( \pi(p) /p + \sqrt{\log p / ( n m^{1/2 } +m ) } )$ for some large fixed constant $C_{\varpi}$, then we have \begin{eqnarray} && \|\widehat{\bfsym \Gamma} ^s - \bfsym \Gamma ^s \| _2 = O_p ( \pi (p) \varpi_m ^{1-\delta} ), \\ && \|\widehat{\bfsym \Gamma} ^s - \bfsym \Gamma ^s \| _{\max} = O_p ( \varpi_m ), \\ && \|\widetilde{\bfsym \Gamma} _{n+1} - \bfsym \Gamma^* \| _{\bfsym \Gamma^*} = O_p \Big ( \tau_m + n^{-1/2} + m^{-1/4} + p^{1/2} (\tau_m ^2+ n^{-1} + m^{-1/2} ) \cr &&\qquad \qquad \qquad \qquad \qquad + \pi(p) \varpi_m ^{1-\delta}\Big ) , \end{eqnarray} where $\bfsym \Gamma^* = \operatorname{E} \(\bfsym \Gamma_{n+1} \middle | \mathcal{F}_n \)$, and the relative Frobenius norm is $\| \bfm A - \bfsym \Gamma^* \| _{\bfsym \Gamma^*} = p^{-1/2} \| \bfsym \Gamma^{*-1/2} (\bfm A - \bfsym \Gamma^* ) \bfsym \Gamma^{*-1/2} \|_F$.
remarkThe concentration inequality condition (ref) has the convergence rate $n^{-1/2} m^{-1/4}+ m^{-1/2}$ which is faster than the usual convergence rate $m^{-1/4}$. The reason is as follows. Usually, we investigate the asymptotic behavior with finite sample period, that is, $n$ is not allowed to go to infinity, and so the convergence rate merely depends on the high-frequency sample size $m$. However, in our setting, we allow the low-frequency sample size $n$ to go to infinity as well. Then the low-frequency summation employs some martingale structure, and we therefore can enjoy the faster convergence rate $n^{-1/2} m^{-1/4}$. The additional term $m^{-1/2}$ is coming from some non-martingale terms such as the drift term. To obtain the sub-Gaussian concentration inequality, we need some sub-Gaussian condition on the observed log stock prices $\bfm Y_t$ such as the bounded instantaneous volatility condition tao2013optimal. Recently, fan2018robust proposed the robust pre-averaged volatility estimation scheme and the sub-Gaussian concentration inequality can be obtained even when the observed log stock prices are heavy-tailed. Thus, this condition is not restrictive.
remarkTheorem (ref) shows that the estimator for future large volatility matrix, $\widetilde{\bfsym \Gamma}_{n+1}$, has the convergence rate of $\tau_m + n^{-1/2} + m^{-1/4} + p^{1/2} (\tau_m ^2+ n^{-1} + m^{-1/2} ) + \pi(p) \varpi_m ^{1-\delta}$. Note that $\tau_m$ depends on the non-parametric estimators in Section (ref) and $\tau_m$ may be $n^{-1/2} + m^{-1/4} + \pi(p) /p$. In this case, the convergence rate will be $n^{-1/2} + m^{-1/4} +\pi(p)/p + p^{1/2} ( n^{-1} + m^{-1/2} ) + \pi(p) \varpi_m ^{1-\delta}$ and the SV-POET estimator $\widetilde{\bfsym \Gamma}_{n+1}$ is consistent as long as $p= o (n^2)$ and $p= o(m)$.

Numerical analysis

A simulation study

In this section, a simulation study was conducted to check finite sample performance of the proposed parameter estimators $\widehat{\ensuremath{\boldsymbol{\theta}}}$ and $\widehat{\ensuremath{\boldsymbol{\theta}}}^{ls}$, as well as to investigate the prediction performance of the proposed SV-POET estimator $\widetilde{\bfsym \Gamma}$, which was also compared with the performance of the estimator for future large volatility matrix proposed in kim2019factor. Let $p$ be the total number of assets studied, $n$ be the total number of low-frequency observations, and $m$ be the total number of high-frequency observations during each low-frequency period. Log prices $\bfm X_t=(X_{1,t}, \ldots, X_{p,t})^{\top}$ at discrete time points $t_{i,j}=i-1+j/m$, $i=1, \ldots, n$ and $j=1, \ldots, m$, were generated according to (ref) with $\ensuremath{\boldsymbol{\mu}}_t=0$ by the Euler scheme. Standard Brownian motions such as $\bfm B_t$ and $\bfm W_t$ were simulated by the normalized partial sums of independent standard normal random variables. We considered a scenario where $r=3$ that suggests three market factors exert an impact on all trading stocks. For the instantaneous factor volatility process in (ref), we considered a diffusion process with $q=1$ so that when the continuous process is restricted to integer times, it retains an AR(1) structure. For the parameters in (ref), we took the following set of values within this simulation study: $\bfsym \alpha_0=\mathrm{Diag}(0.5,0.4,0.3)$, $\mathrm{vec}(\bfsym \alpha_1)=(0.2,0,0,0.5,0.5,-0.2,0.8,-0.5,0.3)^\top$, $\nu=\mathrm{Diag}(0.5,0.5,0.5)$. This set of model parameters results in the following target parameters for estimation:

equation*[equation* omitted — 93 chars of source]
equation*[equation* omitted — 355 chars of source]

Initial values for the instantaneous factor volatility process were $\operatorname{E} \[ \ensuremath{\boldsymbol{\Psi}}_k \]$. The factor loading matrix $\bfm L$ is a $p$-by-$r$ matrix, where the first column takes values $\sqrt{2} \cos \left( 2i \pi/p \right)$, $i=1, \ldots, p$, the second column takes values $\sqrt{2} \sin \left( 2i \pi/p \right)$, $i=1, \ldots, p$, and the third column entries share the same value 1 so that the factor loading matrix retains the structure such that $\bfm L^{\top} \bfm L=p \bfm I_r$. On the other hand, to generate the idiosyncratic diffusion process that has a sparse structure in its daily integrated co-volatility $\bfsym \Gamma^s$, we took $\bfsym \Gamma^s=\left(\Gamma^s_{ij} \right)_{1 \leq i,j \leq p}$, where

equation*[equation* omitted — 88 chars of source]

for the off-diagonal elements and $\Gamma^s_{ii}=0.1, i=1, \ldots, p$, for the diagonal elements. For the high-frequency data $Y_{i,t}$ observed between integer times, we added market microstructure noises to the simulated log price $X_{i,t}$ where the noises were modeled by independent normal random variables with mean 0 and standard deviation $0.005$. Given the simulated log prices $Y_{i,t}$, we employed the PRVM christensen2010pre estimator to obtain daily integrated volatility matrix estimator $\widehat{\bfsym \Gamma}_k$, $k=1,\ldots,n$. The sample variance of PRVM estimators, $\widehat{\mathbb{S}}_{n,m}$, was then computed and its first $r=3$ eigenvectors were adopted to estimate factor loading matrix $\bfm L$. Parameter matrices $\bfsym \beta_{0}$ and $\bfsym \beta_{1}$ were estimated by either maximizing the proposed likelihood function $\widehat{L}_{n,m} (\ensuremath{\boldsymbol{\theta}})$ or minimizing the proposed loss function $\widehat{L}^{ls}_{n,m} (\ensuremath{\boldsymbol{\theta}})$. Parameter estimates from the LSE method were used to initialize the optimization algorithm for the QMLE method. We took $n=125,250,500$ and $m=390,780,2340$ with $p=200$. For each combination of $n$ and $m$, we repeated the simulation for 500 times.

Tables (ref) and (ref) summarize the mean spectral norms, Frobenius norms, and max norms of $\widehat{\bfsym \beta}_0 - \bfsym \beta_{0}$ and $\widehat{\bfsym \beta}_1 - \bfsym \beta_{1}$ given both the QMLE and LSE methods. The results show that as the number of low-frequency or high-frequency observations increases, the estimation performance becomes better, which support the theoretical results derived in Section (ref). Moreover, the QMLE method provides more accurate estimation results than the LSE method. The underlying reason may be that the QMLE method is capable of adjusting the heterogeneous volatility.

table[table omitted — 1,167 chars of source]
table[table omitted — 1,203 chars of source]

The major motivation of our model proposal is to predict future large volatility matrix by taking advantage of the imposed AR model structure at the low-frequency. So we examined the finite sample performance of the proposed estimator $\widetilde{\bfsym \Gamma}_{n+1}$ for the conditional integrated volatility matrix $\operatorname{E} \( \bfsym \Gamma_{n+1} \middle | \mathcal{F}_{n} \)$ based on the procedure described in Section (ref). When estimating the idiosyncratic volatility matrix $\bfsym \Gamma^s$, we applied the threshold $\sqrt{2 \log p/ (n m^{1/2}+ m) }$ on its input $\widetilde{\bfsym \Gamma}^s$. For each simulation, we computed the matrix estimation errors in spectral, max, and relative Frobenius norms respectively:

equation*[equation* omitted — 633 chars of source]

For comparison purpose, we as well examined the prediction performance of the factor and aggregated factor GARCH-It\^o model proposed by kim2019factor. In specific, kim2019factor modeled the eigenvalues of factor volatility matrices by some GARCH-type structure. They also proposed to estimate the factor loading matrix $\bfm L$ in some aggregated form and named the corresponding model as aggregated factor GARCH-It\^o model. For both models, we used $r=3$ and applied threshold $\sqrt{2 \log p/ m^{1/2}}$ for the idiosyncratic volatility matrix estimation. On the other hand, $\operatorname{E} \( \bfsym \Gamma_{n+1} \middle | \mathcal{F}_{n} \)$ has the structure of low-rank plus sparse, thus, we considered the POET procedure introduced by fan2013large to account for such structure. In specific, we chose threshold $\sqrt{2 \log p/ m^{1/2}}$ and used the POET estimator from the previous period $\widehat{\bfsym \Gamma}^{POET}_n$ to estimate $\operatorname{E} \( \bfsym \Gamma_{n+1} \middle | \mathcal{F}_{n} \)$ since when the parametric models are not considered, we often assume martingale structure instead. For the benchmark, we also considered the PRVM estimator $\widehat{\bfsym \Gamma}_n$ from the previous period.

Table (ref) summarizes the mean matrix estimation errors in the spectral, max, and relative Frobenius norms while Figure (ref) plots the mean estimation errors in the relative Frobenius norms against the number, $m$, of high-frequency observations. The proposed SV-POET estimator outperforms the factor and aggregated factor GARCH-It\^o, the POET, and the PRVM methods. The QMLE method provides a bit more accurate prediction results than the LSE method. As the number of low-frequency or high-frequency observations increases, the mean estimation errors decrease for the SV-POET method, which supports the theoretical results in Section (ref). Moreover, the prediction performance of the POET and PRVM only consistently improve given an increasing number of high-frequency observations. This may be because that these estimators are obtained using only the previous period high-frequency observations.

table[table omitted — 3,220 chars of source]
figure[figure omitted — 287 chars of source]

An empirical study

In this section, we demonstrate the proposed prediction methodology with real trading stock prices recorded in minute of $p=200$ companies from January 1st, 2013 to December 31st, 2013. The total number of low-frequency periods follows to be $n=252$ while the daily number of high-frequency returns is $m=390$. We estimated the daily integrated volatility matrix by the PRVM estimator christensen2010pre and projected the obtained PRVM estimators onto the positive semi-definite cone in the spectral norm to ensure their positive semi-definiteness. That is, we set the negative eigenvalues to be $0$. The corresponding PRVM estimates are denoted as $\widehat{\bfsym \Gamma}_k$, $k=1, \ldots, 252$. The sample variance of all PRVM estimators, $\widehat{\mathbb{S}}_{n,m}$, was obtained and its ordered eigenvalues are presented in Figure (ref). Moreover, let $\widehat{\lambda}_{k,j}$ be the $j$th largest eigenvalue of $\widehat{\bfsym \Gamma}_k$, Figure (ref) presents the scree plot based on $\sum_{k=1}^n \widehat{\lambda}_{k,1}/n$, $\sum_{k=1}^n \widehat{\lambda}_{k,2}/n$, $\ldots$, $\sum_{k=1}^n \widehat{\lambda}_{k,p}/n$. Both plots suggest that possible candidates for the number of market factors $r$ is 1,2,3,4 . To determine the rank $r$ specifically, we adopted the procedure as described in ait2017using in the following:

equation*[equation* omitted — 230 chars of source]

where we used $r_{\max}=30$, $c_1=0.02 \widehat{\lambda}_{k,30}$, and $c_2=0.5$. The procedure chose $\widehat{r}=3$.

figure[figure omitted — 219 chars of source]
figure[figure omitted — 198 chars of source]

The AR order $q$ was selected based on standard criteria such as the AIC or BIC, and the fitted model yields $q=1$. We optimized the proposed quasi-likelihood $\widehat{L}_{n,m}(\ensuremath{\boldsymbol{\theta}})$ and the loss function $\widehat{L}^{ls}_{n,m}(\ensuremath{\boldsymbol{\theta}})$ to obtain model parameter estimates as the follows

equation*[equation* omitted — 123 chars of source]
equation*[equation* omitted — 411 chars of source]

and

equation*[equation* omitted — 126 chars of source]
equation*[equation* omitted — 416 chars of source]

The parameter $\bfsym \beta_0$ denotes the intercept term in the factor volatility dynamics and its small estimated values reflect the overall level of daily factor volatilities.

To examine the model prediction performance, we carried out an out-of-sample analysis. In specific, we computed the proposed SV-POET estimator $\widetilde{\bfsym \Gamma}_k$ given observed data from low-frequency period $1$ to $k-1$. To obtain $\widetilde{\bfsym \Gamma}_k$, we first need to estimate the idiosyncratic volatility matrix $\bfsym \Gamma^s$ and in the thresholding step, we used global industry classification standard (GICS) for sectors as guidance fan2015incorporating. Specifically, given the idiosyncratic volatility matrix estimator input $\widetilde{\bfsym \Gamma}^s$, we kept the volatilities within the same sector, but set the others to be zero. The relative prediction errors in various matrix norms: $\| \widetilde{\bfsym \Gamma}_k - \widehat{\bfsym \Gamma}_k \|_2 / \| \widehat{\bfsym \Gamma}_k \|_2$, $\| \widetilde{\bfsym \Gamma}_k - \widehat{\bfsym \Gamma}_k \|_F / \| \widehat{\bfsym \Gamma}_k \|_F$, and $\| \widetilde{\bfsym \Gamma}_k - \widehat{\bfsym \Gamma}_k \|_\text{max} / \| \widehat{\bfsym \Gamma}_k \|_\text{max}$ were examined. Given any forecast origin $h$, we repeated the procedure for the remaining $n-h$ periods and obtained the mean relative prediction errors (MPEs). For comparison purpose, we also studied the factor GARCH-It\^o estimator and the aggregated factor GARCH-It\^o estimator for future volatility matrix $\bfsym \Gamma _k$ as propoesd in kim2019factor. We also used $r=3$ and employed the GICS for the thresholding step. For the benchmark, we as well considered the POET and PRVM methods, and predicted the future volatility matrix $\bfsym \Gamma_k$ by the current volatility matrix estimators $\widehat{\bfsym \Gamma}^{POET}_{k-1}$ and $\widehat{\bfsym \Gamma}_{k-1}$. The GICS was employed for the thresholding step of the POET method fan2013large, fan2015incorporating.

Table (ref) summarizes the MPE values given the SV-POET (LSE or QMLE) estimators, the factor and aggregated factor GARCH-It\^o estimators, the POET and PRVM estimators. To study the dependency of model prediction performance on split points, we report the results for $h=146,168,188$ that correspond to the last trading days of July, August, September in the year 2013. In general, the SV-POET method outperforms the other benchmarks in predicting future volatility matrix. The results are consistent across various split points.

sidewaystable\begin{tabular}{cccccccc} \hline \hline Matrix & Forecast & \multicolumn{2}{c}{SV-POET} & Factor & Aggregated Factor & \multirow{2}{*}{POET} & \multirow{2}{*}{\textbf{PRVM}} \\ \textbf{norms} & \textbf{origins} & QMLE & LSE & \textbf{GARCH-It\^o} & \textbf{GARCH-It\^o} && \\ \hline \multirow{3}{*}{\textbf{Spectral}} & $h=146$ & 0.853 & 0.868 & 1.163 & 0.796 & 1.055 & 1.037 \\ & $h=168$ & 0.831 & 0.846 & 1.248 & 0.838 & 1.085 & 1.066 \\ & $h=188$ & 0.811 & 0.828 & 1.233 & 0.844 & 1.045 & 1.027 \\ &&&&&&& \\ \multirow{3}{*}{\textbf{Frobenius}} & $h=146$ & 0.893 & 0.903 & 1.177 & 0.868 & 1.135 & 1.157 \\ & $h=168$ & 0.882 & 0.892 & 1.229 & 0.898 & 1.157 & 1.179 \\ & $h=188$ & 0.865 & 0.876 & 1.215 & 0.897 & 1.122 & 1.145 \\ &&&&&&& \\ \multirow{3}{*}{\textbf{Max}} & $h=146$ & 0.775 & 0.776 & 1.438 & 0.812 & 1.494 & 1.493 \\ & $h=168$ & 0.785 & 0.787 & 1.290 & 0.811 & 1.325 & 1.324 \\ & $h=188$ & 0.794 & 0.795 & 1.349 & 0.837 & 1.397 & 1.397 \\ \hline \hline \end{tabular} \caption{Mean relative prediction error values for the empirical data set via the (QMLE or LSE), the factor and aggregated factor GARCH-It\^o kim2019factor, the POET fan2013large, fan2015incorporating and PRVM methods for forecast origin $h=146,168,188$. }

We now consider the constrained portfolio allocation problem fan2012vast with the SV-POET estimator $\widetilde{\bfsym \Gamma}_{k}$. Specifically, we minimized the following portfolio risk

equation*[equation* omitted — 184 chars of source]

where $\mathbf{J}=(1, \ldots, 1)^{\top} \in \mathbb{R}^p$ and $c_0$ is the gross exposure constraint which varies from 1 to 2. The portfolio associated with $\widehat{\bfm w}_k$ that minimizes above function is the so-called optimal portfolio. We computed the following out-of-sample portfolio risk for the optimal portfolio in the annualized form

equation*[equation* omitted — 113 chars of source]

where $\widehat{\bfsym \Gamma}_k^{*}$ is the realized variance of day $k$. Given any forecast origin $h$, we repeated the procedure for the remaining $n-h$ periods and obtained the mean out-of-sample portfolio risk as

equation*[equation* omitted — 52 chars of source]

For the benchmarks, we also considered the factor and aggregated factor GARCH-It\^o estimators, the POET estimator, as well as the PRVM estimator.

figure[figure omitted — 410 chars of source]

Figure (ref) plots the annualized mean out-of-sample portfolio risk for various forecast origins $h=146,168,188$ and against different gross exposure constraint $c_0$ values. The results are consistent across the split points. The proposed SV-POET method results in smaller portfolio risk given all split points and exposure constraints comparing to the factor GARCH-It\^o, the POET and PRVM methods. The aggregated factor GARCH-It\^o method performs well when the exposure constraint $c_0$ is small but is unstable when $c_0$ is large. The results suggest that our proposed SV-POET model can capture the market dynamics well by utilizing the AR structure and adopting the historical integrated factor volatilities as innovations for modeling the factor volatility process. On the other hand, the factor and aggregated factor GARCH-It\^o methods assume diagonal factor volatility that is rather restrictive, moreover, they use the GARCH type structure and adopt the squared factor returns as innovations, which may be the reasons for their suboptimal performance. kim2019factor further showed that the performance of their proposed methods could be improved by modeling the idiosyncratic volatility dynamics in their empirical analysis. However, we found that the prediction performance of our proposed SV-POET method could not be enhanced by modeling the idiosyncratic volatility dynamics or by incorporating additional exogenous variables such as overnight factor returns or trading volumes when modeling the factor volatility dynamics. These results indicate that our proposed model alone is sufficient for capturing the market dynamics, while additional information may not be very helpful. Finally, the POET and PRVM methods do not model the volatility matrix process dynamically, which may cause a lack in their empirical performance.

Conclusion

In this paper, we introduce a new method for vast volatility matrix estimation by employing a unified model that can accommodate both the discrete-time stochastic volatility and continuous-time It\^o diffusion models. The proposed SV-It\^o model is capable of studying high-frequency based volatility process in the high-dimensional set-up through a low-dimensional latent factor volatility process that has an autoregressive structure. When estimating the latent factor volatility matrices, the SV-It\^o method assumes a more general structure that is able to account for the cross-sectional market dynamics. We note that this is important as traditional approaches that impose the diagonal assumption are rather restrictive while the corresponding models employed to study the diagonal factor volatility dynamics may not be sufficient for capturing the market dynamics. Model parameters in the SV-It\^{o} model are estimated by either maximizing a quasi-likelihood function or minimizing a squared loss function. The proposed LSE method is easy to implement, however, performs slightly worse than the proposed QMLE method. However, the performance of the QMLE method depends on the initial value selection for the optimization algorithm and may not be very stable. We show that the proposed method presents good performance in predicting future vast volatility matrix and constructing minimum-variance portfolios through our empirical study. When comparing to other existing vast volatility estimation and prediction methods, the proposed model employs the autoregressive structure and historical factor volatilities to explain the factor volatility dynamics, which is more natural and is better supported by the empirical data.

Acknowledgements

The research of Donggyu Kim was supported in part by KAIST Settlement/Research Subsidies for Newly-hired Faculty grant G04170049 and KAIST Basic Research Funds by Faculty (A0601003029). The research of Xinyu Song was supported by the Fundamental Research Funds for the Central Universities (2018110128), China Scholarship Council (201806485017) and National Natural Science Foundation of China (Grant No. 11871323). The research of Yazhen Wang was supported in part by NSF Grants DMS-1528735, DMS-1707605, and DMS-1913149.