Estimating Factor-Based Spot Volatility Matrices with Noisy and Asynchronous High-Frequency Data

Estimating Factor-Based Spot Volatility Matrices with Noisy and Asynchronous High-Frequency Data

\centerline{\bf Abstract}

We propose a new estimator of high-dimensional spot volatility matrices satisfying a low-rank plus sparse structure from noisy and asynchronous high-frequency data collected for an ultra-large number of assets. The noise processes are allowed to be temporally correlated, heteroskedastic, asymptotically vanishing and dependent on the efficient prices. We define a kernel-weighted pre-averaging method to jointly tackle the microstructure noise and asynchronicity issues, and we obtain uniformly consistent estimates for latent prices. We impose a continuous-time factor model with time-varying factor loadings on the price processes, and estimate the common factors and loadings via a local principal component analysis. Assuming a uniform sparsity condition on the idiosyncratic volatility structure, we combine the POET and kernel-smoothing techniques to estimate the spot volatility matrices for both the latent prices and idiosyncratic errors. Under some mild restrictions, the estimated spot volatility matrices are shown to be uniformly consistent under various matrix norms. We provide Monte-Carlo simulation and empirical studies to examine the numerical performance of the developed estimation methodology.

{\em Keywords}: continuous semimartingale, kernel smoothing, microstructure noise, PCA, spot volatility, time-varying factor models.

Introduction

\setcounter{equation}{0}

In high-frequency financial econometrics, the so-called realised volatility has been commonly used to measure the integrated volatility of asset returns over a fixed time window BS02, BS04, ABDL03, S05, AJ14. However, this results in a question of how to choose the time window, in particular when the financial market is volatile. In practice, it is often important to recover the actual spot/instantaneous volatility structure, which plays an important role in various applications such as testing price jumps LM07 and estimating stochastic volatility models KK16, BR18. There have been many studies about spot volatility estimation. For the case of a single asset without market microstructure noise, FW08 and K10 use a nonparametric kernel smoothing method to estimate the spot volatility function and derive its in-fill asymptotic properties. For the more general high-frequency data setting with microstructure noise, ZB14 propose a local version of the two-scale realised volatility ZMA05 to estimate the spot volatility, whereas KK16 combine classic kernel smoothing with the pre-averaging method JLMPV09, CKP10.

Nowadays, practitioners often have to work with high-frequency financial data collected for a large number of assets. The aforementioned spot volatility estimation methods developed for a single or finite number of assets do not generally work well in the high-dimensional and high-frequency data setting. Under a uniform sparsity condition, BLLW23 estimate high-dimensional spot volatility matrices when the number of assets is ultra large, and derive the uniform convergence properties via the joint in-fill and increasing dimensionality asymptotics. However, the sparsity assumption imposed on large volatility matrices is too restrictive, since the price processes are often highly correlated between a large number of assets (in particular those from the same sector). It is well known that there may exist co-movements between these highly-correlated asset prices, and these co-movements may be captured by some latent risk factors. Hence, to relax the restrictive sparsity assumption and estimate meaningful volatility structures, the following continuous-time factor model is often employed for a $p$-dimensional vector of asset prices:

where $\Lambda$ is a $p\times k$ matrix of constant factor loadings, $F_t$ and $ U_t$ are $k$-dimensional and $p$-dimensional continuous semimartingales, respectively (see Section (ref) for the definition). Model ((ref)) is the approximate factor model, which has been extensively studied for low-frequency data CR83, BN02. By imposing a sparse structural assumption on the volatility of $ U_t$, it follows from ((ref)) that $X_t$ has a low-rank plus sparse volatility structure, which is often called conditional sparsity FLM13. FFX16 estimate the large (integrated) volatility matrix of $X_t$ when the factors $F_t$ are observable; whereas AX17 use the principal component analysis (PCA) to estimate the factor model ((ref)) and further construct the large volatility matrix estimation for $X_t$ when the factors are latent. P19 estimates the factor model allowing jumps in the latent factor process and derives the convergence rates and limit distribution theory for the PCA estimated factors and loadings. DLX19 estimate the conditionally sparse large covariance matrices and their inverse for the asynchronous high-frequency data which may be contaminated by the microstructure noise, combine the pre-averaging and generalised shrinkage in the estimation procedure and cover three different scenarios for the factor model specification. Other recent developments on estimation and testing of model ((ref)) can be found in KL18 and SX22.

The continuous-time factor model ((ref)) is essentially static with constant factor loadings. This model assumption may be insufficient when the main interest lies in the spot volatility matrix estimation. In particular, it becomes invalid when there are smooth structural changes or breaks in the process governing the large data series. This motivates the following time-varying factor model for continuous-time processes:

where $\Lambda(t)$ is a $p\times k$ matrix of time-varying factor loading processes and the other components are the same as those in ((ref)). Model ((ref)) covers model ((ref)) as a special case when $\Lambda(t)=\Lambda$. It can be regarded as a natural extension of the time-varying factor model from the low-frequency data setting MHvS11, SW17 to the high-frequency data setting. The aim of this paper is to estimate the large spot volatility matrix for $X_t$ based on the time-varying factor model ((ref)). K18 generalises FLM13's POET (Principal Orthogonal complEment Thresholding) method to estimate the large spot volatility matrix of $X_t$ defined in ((ref)) and its decomposition into the systematic and idiosyncratic volatility components for noise-free and synchronised high-frequency data. In practice, it is not uncommon that the large high-frequency data are non-synchronised and may be contaminated by the market microstructure noise. A direct application of K18's method in the latter setting would result in biased volatility estimation. Hence, a substantial extension is required to address this problem and consistently estimate large spot volatility matrices. For the noisy high-frequency data but with observed (and possibly noise contaminated) factors, two recent papers by BLLW23 and CMZ23 estimate the spot volatility structure of $X_t$ under ((ref)).

We assume the asset prices are observed with an additive noise structure to be defined in Section (ref), where the microstructure noise is allowed to be temporally correlated with nonlinear heteroskedasticity, asymptotically vanishing and dependent on the latent prices. The empirical studies in JLZ17 and LL22a reveal non-trivial temporal dependence for the microstructure noise of individual asset prices. The assumption of nonlinear heteroskedasticity allows the microstructure noise to have prominent intraday patterns such as the well-known U-shape. Some authors have argued that the microstructure noise is small in magnitude, and so, as in KL08, we allow the microstructure noise to be shrinking as the data frequency increases. Such a “small noise" model structure is also adopted by DX21 and LL22b. Furthermore, we allow the noise to be correlated with the latent prices but do not impose any explicit correlation structure KL08. The existence of endogenous noise may be due to rounding effects, price stickiness and asymmetric information HL06.

In practice different assets are traded at different frequencies and so their transaction prices are observed at times that are not synchronised and can be quite differently spaced due to liquidity levels varying over assets. This often leads to volatility matrix estimation bias and possibly enhances the so-called Epps effect E79. In order to jointly tackle the noise and asynchronicity issues, we extend the conventional pre-averaging method and introduce a kernel-weighted version, i.e., we take a weighted average of the observed prices over a local neighborhood of $t$ by kernel smoothing to obtain an approximation of the latent price at time $t$. This kernel-weighted pre-averaging has been used by KK16 in spot volatility estimation for a single asset, and recently has been adopted by BLLW23 in large spot volatility matrix estimation for synchronised high-frequency data. We show that the the first-step kernel filter of the asset prices is uniformly consistent with the approximation order determined by the bandwidth, the dimension and the sample size, extending the uniform consistency derived by KK16 and BLLW23 to a much more general setting.

To estimate the latent structure in the time-varying factor model ((ref)), it is natural to apply a local version of PCA which relies on the local expansion of the smooth time-varying factor loading functions MHvS11, SW17, K18, WPLL21. Since the latent prices are unavailable, we have to apply the local PCA to the estimated prices obtained via pre-averaging. This extension is non-trivial and leads to an extra challenge in proving the in-fill asymptotics for the local PCA estimates of common factors and factor loadings. The mean squared convergence rates for the estimated factors (taking the first-order difference) and the uniform convergence rates for the estimated time-varying factor loadings are slower than those in K18 who considers synchronised noise-free data. This is reasonable, since approximating latent prices results in extra estimation errors that consequently slow down the convergence rates for local PCA.

It follows from ((ref)) that the spot volatility matrix of $X_t$ is decomposed into the common and idiosyncratic spot volatility matrices. The latter is assumed to satisfy a uniform sparsity condition, which is similar to that used by CXW13, CL16, CLL19 and BLLW23, and results in the low-rank plus sparse matrix structure. Accordingly, we combine the classic POET with kernel smoothing to estimate the spot volatility matrices for both the latent prices and idiosyncratic errors. In particular, a generalised shrinkage technique is applied to the idiosyncratic spot volatility matrix estimation. We derive the uniform convergence property for the estimated idiosyncratic spot volatility matrix in the (elementwise) max and spectral norms, and obtain the rates comparable to those derived in CMZ20 and BLLW23. Due to the spiked volatility structure, we derive the uniform convergence rates for the estimated spot volatility matrix of $X_t$ not only in the max norm but also in the relative error measurement FLM13, WPLL21.

The finite-sample Monte-Carlo simulation studies show that the developed large spot volatility estimation method outperforms the naive estimation (without applying shrinkage to the estimated spot idiosyncratic volatility matrix). An empirical application to the one-min intraday log-price of S&P 500 index constituents reveals significant time-varying patterns of the spot volatility and covariance, and demonstrates rationality of the low-rank plus sparse spot volatility structure. In particular, the estimated factor number varies over time with fewer factors during the period of market collapse and more factors when the market is stable.

The rest of the paper is organised as follows. Section (ref) introduces the additive microstructure noise model, low-rank plus sparse spot volatility matrix structure and estimation technique for large spot volatility matrices. Section (ref) gives some regularity conditions and the in-fill asymptotic properties for the developed estimates together with some remarks. Section (ref) conducts the Monte-Carlo simulation study and Section (ref) reports the empirical application. Section (ref) concludes the paper. Proofs of the main theoretical results and some technical lemmas are available in the Supplementary Material LLZ2024. Throughout the paper, we let $\Vert\cdot\Vert$ be the Euclidean norm of a vector; and for a $p\times p$ matrix $A=(A_{i_1i_2})_{p\times p}$, we let $\Vert A\Vert_s$ and $\Vert A\Vert_F$ be the matrix spectral norm and Frobenius norm, $\vert A\vert_1=\sum_{i_1=1}^p\sum_{i_2=1}^p |A_{i_1i_2}|$, $\Vert A\Vert_1=\max_{1\leq i_2\leq p}\sum_{i_1=1}^p |A_{i_1i_2}|$, $\Vert A\Vert_{\infty,q}=\max_{1\leq i_1\leq p}\sum_{i_2=1}^p |A_{i_1i_2}|^q$ and $\Vert A\Vert_{\max}=\max_{1\leq i_1\leq p}\max_{1\leq i_2\leq p} |A_{i_1i_2}|$.

Model and methodology

\setcounter{equation}{0}

In this section, we first introduce an additive noise structure for contaminated high-frequency data and the low-rank plus sparse spot volatility matrix structure, followed by description of the estimation methodology.

Microstructure noise

For the $i$-th asset, suppose that the asset prices are observed with the additive noise structure:

where $X_{i,t}$ is the $i$-th component of $X_t$, $t_1^i,\cdots,t_{n_i}^i$ are the data collection time points (which may be non-equidistant) for the $i$-th asset, and

is the noise with $\chi_i(\cdot)$ and $\varepsilon_{i,j}^\ast$ satisfying Assumption (ref) in Section (ref) below. The general structure in ((ref)) shows that the microstructure noise not only has a nonlinear heteroskedastic structure $\chi_i(\cdot)$, but also is asymptotically vanishing as the sampling frequency increases if $0<\beta_i<1/2$. In addition, we allow $\varepsilon_{i,j}^\ast$ to be correlated over $i$ and $j$, and dependent on the latent prices $X_{i,t_j^i}$. Throughout the paper, we let $t_{0}^i\equiv0$ and $t_{n_i}^i=T_i$.

The latent vector process $X_t=(X_{1,t},\cdots,X_{p,t})^{^\intercal}$ satisfies the time-varying factor model structure ((ref)), where $F_t$ and $ U_t$ are $k$-dimensional and $p$-dimensional Brownian semimartingales, respectively, solving

and

$\mu_t^F$ and $\mu_t^U$ are vectors of drift with dimensions $k$ and $p$, respectively, $\sigma_t^F$ and $\sigma_t^U$ are $k\times k$ and $p\times p$ matrices of spot volatilities, $W_t^F$ and $W_t^U$ are $k$-dimensional and $p$-dimensional standard Brownian motions. The drifts and spot volatilities are progressively measurable processes. As in K18 and P19, we assume that $\{W_t^F: t\geq0\}$ and $\{W_t^U: t\geq0\}$ are independent standard Brownian motions. This assumption may be relaxed to allow for weak correlation between the factor and idiosyncratic error processes. Without loss of generality, we may further assume that $\sigma_t^F=I_k$, a $k\times k$ identity matrix; otherwise, we can re-define the factor loading matrix as $\Lambda(t) \sigma_t^F$ in ((ref)) and the drift as $(\sigma_t^F)^{-1}\mu_t^F$ in ((ref)).

Large spot volatility matrices with conditional sparsity

Our main interest lies in estimating the factor-based spot volatility structure of the latent process $X_t$. For $0<\tau<T$ with $T$ being fixed, we let $\Sigma_X(\tau)$, $\Sigma_F(\tau)$ and $\Sigma_U(\tau)$ denote the spot volatility matrices for $X_t$, $F_t$ and $ U_t$, respectively. From ((ref)), ((ref)) and ((ref)), assuming $F_t$ and $U_t$ are independent, we readily have that

indicating that $\Sigma_X(\tau)$ can be decomposed into the common and idiosyncratic spot volatility matrices: $\Sigma_C(\tau)$ and $\Sigma_U(\tau)$. Throughout the paper, we assume the following uniform sparsity condition: $\left\{\Sigma_U(t):\ 0\leq t\leq T\right\}\in \mathcal{S}(q,\varpi_p)$ which is defined by

where “$\succ0$" denotes positive definiteness, $0\le q<1$ and $\Psi$ is a positive random variable satisfying ${\mathsf E}(\Psi)\leq C_\Psi<\infty$. This is similar to the sparsity assumption used by CXW13, CL16, CLL19 and BLLW23 and is a natural extension of the approximate sparsity BL08.

With the sparsity restriction on $\Sigma_U(\cdot)$, we obtain a low-rank plus sparse (or conditionally sparse) structure on $\Sigma_X(\cdot)$. In the low-frequency data setting, FLM13 introduce the POET method to estimate the large covariance matrix of $X_t$ which satisfies the discrete-time approximate factor model; and WPLL21 propose a local POET via kernel smoothing to estimate the large dynamic covariance matrix of $X_t$ which satisfies the state-varying or time-varying factor model. Extending their methods to estimate $\Sigma_X(\cdot)$ is non-trivial since $X_t$ is latent and $Y_{i,t}$ is non-synchronised. K18 generalises the POET method to estimate the spot volatilities $\Sigma_X(\cdot)$, $\Sigma_C(\cdot)$ and $\Sigma_U(\cdot)$ as well as their integrated versions for noise-free and synchronised high-frequency data. His method is not applicable to our high-frequency data setting and would result in biased volatility estimation due to presence of the microstructure noise and asynchronicity in model ((ref)).

Estimation of the large spot volatility matrix

Partly motivated by KK16 and BLLW23, we next introduce a kernel-weighted pre-averaging method to jointly tackle the microstructure noise and asynchronicity issues. The proposed technique is a local extension of the conventional pre-averaging which is introduced by JLMPV09 and CKP10 to estimate the integrated volatility for a single asset and has been further extended by KWZ16 and DLX19 to estimate large matrices of integrated volatility.

Start with locally averaging the high-frequency observations $Y_{i,t_j^i}$ via a kernel filter:

where $L(\cdot)$ is a kernel function, $b$ is a bandwidth and $L_b(\cdot)=b^{-1}L(\cdot/b)$. The $\widetilde{X}_{i,t}$ can be seen as a fitted value for the latent $X_{i,t}$, and Proposition (ref) in Section (ref) below establishes the uniform consistency property for $\widetilde{X}_{i,t}$. We adopt the modified kernel weights $(t_j^i-t_{j-1}^i)L_b(t_j^i-t)$ in ((ref)) to address the issue of irregular/asynchronous sampling times. For the special case of equally-spaced time points in high-frequency data collection, i.e., $t_j^i-t_{j-1}^i\equiv\Delta$, $j=1,\cdots,n_i$, the kernel filter in ((ref)) can be simplified to \[\widetilde{X}_{i,t}=\Delta\sum_{j=1}^{n_i} L_b(t_j^i-t)Y_{i,t_j^i},\] see KK16.

Let $\widetilde X_t=(\widetilde{X}_{1,t},\cdots, \widetilde{X}_{p,t})^{^\intercal}$ and $t_j=j\Delta_\circ$, $j=0,1,\cdots,N$, with $\Delta_\circ$ being a user-specified time difference and $N=\lfloor T/\Delta_\circ\rfloor$. We construct the kernel-weighted realised volatility matrix:

where \[\Delta \widetilde X_{j}=\widetilde X_{t_j}-\widetilde X_{t_{j-1}}=\left(\Delta \widetilde X_{1,j}, \cdots, \Delta \widetilde X_{p,j}\right)^{^\intercal}.\] With eigen-analysis on $\widetilde\Sigma_X(\tau)$, we obtain $(\widetilde{\lambda}_{l}(\tau),\widetilde{\eta}_{l}(\tau))$, $l=1,\cdots,p$, as pairs of eigenvalues and normalised eigenvectors, where the eigenvalues are arranged in a descending order. From the matrix spectral decomposition, assuming the factor number is known a priori, we re-write

where $\widetilde\Sigma_C(\tau)$ and $\widetilde\Sigma_U(\tau)$ are the estimated spot volatility matrices at time $\tau$ for the common and idiosyncratic error components, respectively. By the construction in ((ref)), $\widetilde\Sigma_C(\tau)$ is a low-rank covariance matrix. With the sparsity restriction on $\Sigma_U(\cdot)$, we may further apply a generalised shrinkage. Let $s_\rho(\cdot)$ be a shrinkage function satisfying that (i) $\vert s_\rho(u)\vert \leq \vert u\vert$ for $u\in{\cal R}$; (ii) $s_\rho(u)=0$ if $\vert u\vert \leq \rho$; and (iii) $\vert s_\rho(u)-u\vert\leq \rho$, where $\rho$ is a user-specified tuning parameter controlling the level of shrinkage. Letting $\widetilde{\sigma}_{U,i_1i_2}(\tau)$ be the $(i_1,i_2)$-entry of $\widetilde\Sigma_U(\tau)$, we define

where $\rho(\cdot)$ is a time-varying tuning parameter. Consequently, we obtain the kernel POET estimate of $\Sigma_X(\tau)$:

It follows from Theorem 1 in FLM13 and Proposition 1 in WPLL21 that the kernel POET estimate defined in ((ref)) is equivalent to the local PCA method to be introduced shortly. By ((ref)) with some smoothness condition on $\Lambda(\cdot)$ (see Assumption (ref)(iii) in Section (ref)), we may obtain the following local approximation of the time-varying factor model using the increments of the stochastic processes:

where $\Delta X_{j}=X_{t_j}-X_{t_{j-1}}=(\Delta X_{1,j}, \cdots, \Delta X_{p,j})^{^\intercal}$, $\Delta F_j$ and $\Delta U_j$ are defined analogously. As in ((ref)), we replace $\Delta X_j$ by $\Delta \widetilde X_j$ in the following local PCA. Write $\Lambda(\tau)=[\Lambda_{1}(\tau),\cdots,\Lambda_{p}(\tau)]^{^\intercal}$ with $\Lambda_{i}(\tau)=[\Lambda_{i,1}(\tau),\cdots,\Lambda_{i,k}(\tau)]^{^\intercal}$,

Define the kernel-weighted least squares objective function:

where $\overline\Lambda$ and $\overline F_j$ are generic notation for the factor loading matrix and common factor in the local PCA estimation procedure, and $\overline F(\tau)$ is defined similarly to $\Delta F(\tau)$ but with $\Delta F_j$ replaced by $\overline F_{j}$.

Consider the following identification condition:

which are commonly used in PCA estimation of the factor models BN02,FLM13, CMZ20. With the identification condition ((ref)), we replace $\overline\Lambda$ in ((ref)) by $\Delta\widetilde X(\tau)\overline F(\tau)$. Then, the kernel-weighted objective function in ((ref)) becomes \[ \mathsf{trace}\left\{\Delta\widetilde X(\tau)^{^\intercal}\Delta\widetilde X(\tau)\right\}-\mathsf{trace}\left\{\overline F(\tau)^{^\intercal}\Delta\widetilde X(\tau)^{^\intercal}\Delta\widetilde X(\tau)\overline F(\tau)\right\}, \] indicating that minimising ((ref)) subject to the restriction ((ref)) is equivalent to maximising the trace of $\overline F(\tau)^{^\intercal}\Delta\widetilde X(\tau)^{^\intercal}\Delta\widetilde X(\tau)\overline F(\tau)$ subject to $\overline F(\tau)^{^\intercal}\overline F(\tau)=I_k$. Hence, we conduct the eigen-analysis on the $N\times N$ matrix $\Delta\widetilde X (\tau)^{^\intercal}\Delta\widetilde X(\tau)$ and obtain the local PCA estimate of $\Delta F(\tau)$: \[\Delta\widetilde F(\tau)=\left[\Delta\widetilde F_1(\tau),\cdots,\Delta\widetilde F_N(\tau)\right]^{^\intercal}\] which is a matrix consisting of the eigenvectors corresponding to the $k$ largest eigenvalues. Furthermore, the time-varying factor loading matrix $\Lambda(\tau)$ is estimated as \[\widetilde\Lambda(\tau)=\Delta\widetilde X(\tau)\Delta\widetilde F(\tau)=\left[\widetilde\Lambda_1(\tau),\cdots,\widetilde\Lambda_p(\tau)\right]^{^\intercal}.\] The kernel-weighted residuals $\Delta U_j K_h^{1/2}(t_j-\tau)$ in the local approximation ((ref)) are then approximated by \[\widetilde U_j(\tau)=\left[\widetilde{U}_{1,j}(\tau),\cdots,\widetilde{U}_{p,j}(\tau)\right]^{^\intercal}=\Delta\widetilde X_j K_h^{1/2}(t_j-\tau)-\widetilde\Lambda(\tau)\Delta\widetilde F_j(\tau),\] which are subsequently used to estimate $\Sigma_U(\cdot)$. However, the conventional sample covariance matrix using $\widetilde U_j(\tau)$:

usually performs poorly when the number of assets $p$ is ultra large. To address this problem, we again apply the generalised shrinkage and estimate $\Sigma_U(\tau)$ by

Finally, the estimate of $\Sigma_X(\tau)$ is obtained via

It is worth comparing our model assumption and methodology with those in CMZ20 before concluding this section. Although our main model framework is similar to that in CMZ20, we impose more general assumptions on the microstructure noises, allowing them to be nonlinear heteroskedastic and asymptotically vanishing. In particular, the noises can be endogenous but no explicit correlation structure (between noises and latent prices) is required. The estimation methodology in CMZ20 is mainly built on the smoothed two-scale realised volatility introduced in CMZ19, which combines the pre-averaging and two-scale realised volatility. This is further combined with local POET to estimate spot volatility matrices in the high dimension. In contrast, our proposed methodology is virtually simpler, using modified kernel weights in the first step of pre-averaging and kernel POET (or local PCA) in the second step.

Main theoretical results

\setcounter{equation}{0}

In this section, we give some regularity conditions with remarks and then present the theoretical properties for the developed large spot volatility matrix estimates via the in-fill asymptotics.

Regularity conditions

Write

\setcounter{remark}{0}

Asymptotic theorems

We next present the in-fill asymptotic properties of the proposed estimates, letting $t_j^i-t_{j-1}^i\rightarrow0$ and thus $n_i\rightarrow\infty$ for each asset. Proposition (ref) below gives the uniform approximation order of the latent price estimates by kernel-weighted pre-averaging.

\setcounter{prop}{0}

We next derive the uniform convergence property of the estimated idiosyncratic spot volatility matrix $\widehat\Sigma_U(\cdot)$ in both the max and spectral norms. Since $\widetilde\Sigma_U(\cdot)$ is equivalent to $\widehat\Sigma_U(\cdot)$, the following theorem continues to hold for $\widetilde\Sigma_U(\cdot)$.

\setcounter{theorem}{0}

We next explore the uniform convergence property of $\widehat\Sigma_X(\cdot)$. Due to the time-varying factor model structure ((ref)) with the condition ((ref)), the largest $k$ eigenvalues are spiked, diverging at a rate of $p$. Consequently, the spot volatility matrix $\Sigma_X(\cdot)$ cannot be consistently estimated in the absolute term. Motivated by FLM13, we measure the matrix estimate $\widehat\Sigma_X(\tau)$ in the following (time-varying) relative error: \[\left\Vert \widehat\Sigma_X(\tau)-\Sigma_X(\tau)\right\Vert_{\Sigma_X(\tau)}=\frac{1}{\sqrt{p}}\left\Vert \Sigma_X^{-1/2}(\tau)\widehat\Sigma_X(\tau) \Sigma_X^{-1/2}(\tau)-I_{p} \right\Vert_F.\]

Monte-Carlo simulation studies

\setcounter{equation}{0}

Data generating process

We generate the latent price process $X_t$ from a drift-free time-varying factor model:

where ${\Lambda}(t)$ is a $p\times k$ time-varying factor loading process, ${\sigma}^U_t$ is a $p\times p$ (diagonal) idiosyncratic volatility process, ${W}_t^F$ is a $k$-dimensional standard Brownian motion and ${W}_t^U$ is a $p$-dimensional Brownian motion with covariance matrix $\Sigma_\rho$. The number of latent factors is $k=3$. The time-varying factor loading matrix ${\Lambda}(t)=\left[\Lambda_{i,l}(t)\right]_{p\times k}$ is generated from \[ \Lambda^2_{i,l}(t)=b_{i,l}\left[a_{i,l}-\Lambda^2_{i,l}(t)\right] d t+\sigma_{i,l}^{0} \Lambda_{i,l}(t) d W_{i,t}^F, \quad i=1,\cdots,p,\ \ l=1, 2, 3, \] where $a_{i,1}=0.01+i/p, a_{i,2}=0.0115+i/p, a_{i,3}=0.0105+i/p, b_{i,1}=0.006+i/(100 p), b_{i,2}=0.007+i/(100 p), b_{i,3}=0.008+i/(100 p), \sigma_{i,1}^{0}=0.3+i/(5p), \sigma_{i,2}^{0}=\sigma_{i,3}^{0}=0.4+i/(5 p)$, and $W_{i,t}^F$ is the $i$-th element of $W_t^F$. Define $\sigma^U_t = \mathsf{diag}(\sigma^U_{1,t},\cdots,\sigma^U_{p,t})$ with \[ (\sigma_{i,t}^{U})^2=\left[0.00053+i /(100p)\right]\left[0.0017+i / p- (\sigma_{i,t}^{U})^2\right] d t+\left[0.0013+i/(10 p)\right] \sigma_{i,t}^{U} d W_{i,t}^U, \] where $W_{i,t}^U$ is the $i$-th element of $W_t^U$. Choose $\Sigma_\rho=(\rho_{ij})_{p\times p}$ as a correlation matrix satisfying the banded structure GA2017:

where $\mathsf {I}(\cdot)$ is an indicator function and $\rho \sim U(0,0.5)$. Consequently, the generated spot volatility structure has the following low-rank plus sparse decomposition:

The number of assets is set to be $p=100,300,500$ and the simulation is repeated for $100$ times. The sampling frequency $\Delta$ is set to be $1/(6.5\times 60 \times 6)$, which is equivalent to sampling in every $10$ seconds, setting $T$ as one trading day and assuming there are $6.5$ trading hours for each trading day.

In practice, different asset prices are often non-synchronized and contaminated by the micro-structure noise. We generate the contaminated data via model ((ref)). Similar to FK19, we generate $\varepsilon_{i,t}$ independently (over $i$ and $t$) from $\sigma_{\epsilon}\cdot {\mathsf N} (0, \sigma_{ii}(t))$, where $\sigma_{ii}(t)$ is the $i$-th diagonal element of the true spot volatility matrix $\Sigma_X(\cdot)$ at time point $t$ and $\sigma_{\epsilon}$ is the noise-to-signal ratio. The noise is thus correlated with the latent prices. We consider $\sigma_\epsilon=0.05, 0.1$ and $0.2$. To simulate the asynchronous data, we use $p$ independent Poisson processes to generate the true data collection time points $t_j^i$ BN2011. Specifically, for the $i$-th asset, we control $\{t_j^i: j=1,\cdots,n_i\}$ by a Poisson process with mean $\lambda^i \sim U(1,3)$, indicating that asset prices are randomly collected every 10 to 30 seconds and we expect to have $2340/\lambda^i$ observations on average for the $i$-th asset.

Spot volatility matrix estimates

Throughout the simulation studies, we only use the kernel POET estimation method since it is equivalent to the local PCA. To assess the impacts of micro-structure noise and asynchronicity, we consider the following two spot volatility matrix estimators.

We employ the following four shrinkage techniques in the spot idiosyncratic volatility matrix estimation: Smoothly Clipped Absolute Deviation (SCAD), adaptive-lasso (A-lasso), soft thresholding (Soft) and hard thresholding (Hard). In addition, to demonstrate effectiveness of the generalised shrinkage, we also consider the naive estimator as a benchmark that does not apply shrinkage to the spot idiosyncratic volatility matrix estimation.

For each $\tau$, we use the eigenvalue-ratio criterion AH2013 to determine the number of latent factors. Let $\overline{\Sigma}_{X}^{(r)}(\cdot)$ and $\overline{\Sigma}_{U}^{(r)}(\cdot)$ be generic notation for estimates of ${\Sigma}_{X}^{(r)}(\cdot)$ and ${\Sigma}_{U}^{(r)}(\cdot)$, respectively, for the $r$-th replication. To measure the distance between the true spot volatility matrices and the estimated ones, we compute the following two measurements: the mean spectral norm deviation ($\mathsf{MSN}_X$) and the mean relative norm deviation ($\mathsf{MRN}_X$) defined by

where $\| \overline{\Sigma}_{X}^{(r)}(\tau) - {\Sigma}_{X}^{(r)}(\tau) \|_{\Sigma_X^{(r)}(\tau)}$ is defined as in Theorem (ref) for the $r$-th replication, and $\tau_j, j=1,2,\cdots,10$, are equidistant time points. To compare the performance of various shrinkage approaches in estimating idiosyncratic spot volatility matrices, we compute the following mean spectral norm for the idiosyncratic spot volatility matrices: \[ \mathsf{MSN}_U =\frac{1}{100} \sum_{r=1}^{100} \left( \frac{1}{10} \sum_{j=1}^{10} \left \| \overline{\Sigma}_{U}^{(r)}(\tau_j)-\Sigma_U^{(r)}(\tau_j) \right \|_{s} \right). \]

Choice of the tuning parameters

It is well known that the nonparametric kernel-based estimation is sensitive to the bandwidth selection. In the developed estimation procedure, we need to carefully select two bandwidths: $b$ and $h$, which are involved in the kernel-weighted pre-averaging ((ref)) and realised volatility matrix ((ref)), respectively.

As recommended by KK16, we apply the cross-validation to select the optimal bandwidth in the local pre-averaging for each asset, i.e., for the $i$-th asset, obtain \[ b^i_\mathrm{opt}= \operatorname*{arg\,min}_{b>0} \mathsf{CV}_i(b),\ \ \mathsf{CV}_i(b)=\sum_{j=1}^{n_i} \left(Y_{i,t^i_j} - \widetilde{X}_{i,-t_i^j}\right)^{2}\mathsf{I}\left\{T_{l}^i \leq t^i_{j} \leq T_{u}^i\right\}, \] where the truncation mechanism is used to avoid the kernel estimation boundary effect, $T_{l}^i=0.05\times T_i$, $T_{u}^i=0.95\times T_i$, and $\widetilde{X}_{i,-t_i^j}$ is the leave-one-out kernel filter defined as in ((ref)) but removing the observation at $t_i^j$.

In the low-dimensional data setting, we may select the optimal bandwidth entry by entry in the realised volatility matrix estimation. For example, as in K10 and KK16, for $1\leq i_1,i_2\leq p$, we obtain \[ h^{i_1i_2}_\mathrm{opt}= \operatorname*{arg\,min}_{h>0} \mathsf{CV}_{i_1i_2}(h),\ \ \mathsf{CV}_{i_1i_2}(h)=\sum_{\tau_j} \left[ \left(\Delta\widetilde X_{\tau_j} \Delta\widetilde X_{\tau_j}^\prime / \Delta \right)_{i_1i_2} -\widetilde{\sigma}_{X,i_1i_2,-\tau_j}\right]^{2}, \] where $\widetilde{\sigma}_{X,i_1i_2,-\tau_j}$ is the ($i_1,i_2$)-entry of leave-one-out kernel-weighted realised co-volatility estimate. However, this bandwidth selection criterion is rather time consuming when the dimension $p$ is large. Hence, we adopt the following criterion: \[ h_\mathrm{opt}= \operatorname*{arg\,min}_{h>0} \mathsf{CV}(h),\ \ \mathsf{CV}(h)=\sum_{i=1}^{p}\sum_{\tau_j}\left[\left(\Delta\widetilde X_{\tau_j} \Delta\widetilde X_{\tau_j}^\prime / \Delta \right)_{ii}-\widetilde{\sigma}_{X,ii,-\tau_j}\right]^{2}, \] where we only evaluate the quadratic loss for the diagonal elements. To further speed up the computation, as recommended by KK16, we may repeat the bandwidth selection processes over a small number of replications (say, $5$) and then fix the optimal bandwidth as the average of the selected bandwidth values.

When applying the generalised shrinkage in the kernel POET, we need to choose an appropriate tuning parameter to control the shrinkage level. As recommended by FLM13 and WPLL21 and discussed in Remark (ref)(iv), we set \[ \rho_{ij}(\tau)=c_\rho\left[{\widetilde{\sigma}_{U,ii}(\tau){\widetilde{\sigma}_{U,jj}(\tau)}}\right]^{1/2}, \ \ i,j=1,\cdots,p \] where $c_\rho$ is the minimum positive number guaranteeing positive definiteness of the estimated spot idiosyncratic volatility matrix. Specifically, $c_\rho$ is chosen such that \[ c_\rho=\inf \left\{c>0:\ \lambda_{\min}\left\{\left(\overline\Sigma_U^c(\tau)\right)\right\}>0,\ \ \overline\Sigma_U^c(\tau)=\left[\overline{\sigma}_{U,ij}^c(\tau)\right]_{p\times p} \right\}, \] where $\overline{\sigma}_{U,ij}^c(\tau)$ is defined as in ((ref)) using $\rho_{ij}(\tau)=c[{\widetilde{\sigma}_{U,ii}(\tau){\widetilde{\sigma}_{U,jj}(\tau)}}]^{1/2}$.

Simulation results

Tables (ref) and (ref) report the $\mathsf{MSN}_{U}$ values of the spot idiosyncratic volatility matrix estimation for the noise-contaminated synchronised and asynchronised data, respectively. In general, the use of shrinkage results in more accurate estimation than the naive one without shrinkage. The improvement by using the shrinkage becomes more significant as the dimension $p$ increases. The SCAD and adaptive-lasso perform somewhat better than the other two shrinkage methods when the high-frequency data are synchronised (see Table (ref)), and the adaptive-lasso outperforms the other shrinkage methods when the data are asynchronised (see Table (ref)). It is unsurprising that the $\mathsf{MSN}_{U}$ values increase as the dimension or the noise-to-signal ratio $\sigma_\epsilon$ increases.

Tables (ref)--(ref) present the estimation results for $\Sigma_X$. Tables (ref) and (ref) report the $\mathsf{MSN}_{X}$ measurements whereas Tables (ref) and (ref) report the $\mathsf{MRN}_{X}$ measurements. As in Tables (ref) and (ref), the spot volatility matrix estimation with the generalised shrinkage outperforms the naive estimation without shrinkage. The difference between the shrinkage estimation and the naive one is more significant in terms of $\mathsf{MRN}_{X}$. Both $\mathsf{MSN}_{X}$ and $\mathsf{MRN}_{X}$ increase as the dimension $p$ grows. In particular, this divergence pattern is more significant for $\mathsf{MSN}_{X}$, which is not uncommon as the true spot volatility matrix is spiked with the largest $k$ eigenvalues diverging at the same rate of $p$. Hence, as recommended in FLM13, it is more sensible to use the $\mathsf{MRN}_{X}$ to access the factor-based spot volatility matrix estimation. Among the four shrinkage techniques, the SCAD and adaptive-lasso shrinkage slightly outperform the soft and hard thresholding methods, and the hard thresholding often has the poorest performance (although it is still superior to the naive estimator). Due to the data asynchronisation, the $\mathsf{MSN}_{X}$ and $\mathsf{MRN}_{X}$ values in Tables (ref) and (ref) are larger than those provided in Tables (ref) and (ref). The patterns of the estimation performance and comparison are similar over different noise-to-signal ratios. In addition, both the $\mathsf{MSN}_{X}$ and $\mathsf{MRN}_{X}$ measurements increase as $\sigma_\epsilon$ increases from $0.05$ to $0.2$, which is expected as “larger" noise results in “larger" approximation errors of the first-step pre-averaging.

An empirical application

\setcounter{equation}{0}

We next apply the proposed approach to the one-min intraday log-price of S$\&$P 500 index constituents. The log-prices are extracted from Bloomberg, ranging from 29 March to 30 June 2021. There are 66 trading days and 505 stocks in total. The 66 trading days are divided into 11 equal-length time intervals, with 6 trading days in each time interval. In the preliminary data analysis, we remove the illiquid stocks that were not traded for more than half of one-minute trading periods on each trading day, resulting in $353$ stocks in total. We also exclude stock prices that are outside the trading hours of 9:30 am (EST) to 4 pm (EST). A recent paper by LCL2023 estimates the latent factor structure for the same data set, allowing the existence of microstructure noises. In this section, we aim to further explore the spot volatility structure and study its time-varying pattern. For each equal-length time interval, we first apply the kernel-weighted pre-averaging method in ((ref)) to the (noisy) log-price and then construct the kernel-weighted spot volatility matrix estimates in ((ref)). Due to a possible low-rank plus sparse structure, we further implement the kernel POET estimation as in ((ref)). Specifically, we estimate $36$ evenly separated spot volatility matrices in each time interval (i.e., 6 spot volatility matrix estimates for each trading day). As recommended in the simulation, we adopt the adaptive-lasso in the generalised shrinkage with the tuning parameter chosen according to the criteria discussed in Section (ref).

The left panel of Figure (ref) depicts the log-price of the POOL.OQ stock over the sampling period whereas the right panel depicts the fitted log-price via the kernel-weighted pre-averaging. We use the Gaussian kernel function in ((ref)) and determine the optimal bandwidth via the cross-validation discussed in Section (ref). It follows from Figure (ref) that the fitted log-prices are generally close to the original ones which may be noise contaminated, capturing most of the dynamic patterns. Due to the kernel-weighted smoothing, the fitted prices are slightly smoother than the observed ones.

With the kernel POET estimates in ((ref)), we are able to collect the estimated spot volatilities (the diagonal entries of the estimated spot volatility matrix) for all the stocks and the estimated spot covariances (the off-diagonal entries of the estimated spot volatility matrix) between every pair of stocks across all the trading days. We further average the estimated volatilities and covariances within each trading day over the sampling period, and rank them (in decreasing order) according to variances of all the realised volatility and covariance over time. The estimated volatilities and covariances of the 20%, 40%, 60% and 80% quantile stocks (or stock pairs) across all the trading days are plotted in Figures (ref) and (ref), respectively. It follows from the plots that both the spot volatilities and covariances exhibit significant dynamic changes over time, demonstrating that it is imperative to recover the spot volatility structure over time for high-frequency financial data.

Figure (ref) plots the heat maps for the estimated spot correlation matrices on two typical trading days: 19 May and 7 June. The heat maps in the upper panel are obtained for 7 June, which represents a typical period of stable volatility with nine estimated common factors. In contrast, those in the lower panel are for 19 May, which corresponds to the period during the crypto-currency market collapse P21 with only two estimated common factors. We note that the upper-left heat map (with nine estimated factors) is slightly denser than the lower-left one (with two estimated factors), which may be partly explained by the magnitude of the market disruption. During the period of the crypto-currency market collapse, its disruption is so large that it nearly dominates the entire stock market. Consequently, the two estimated factors can explain majority of the market volatility. In contrast, the nine estimated factors during the stable period may only explain a limited portion of the market volatility. Therefore, it is reasonable to expect that the spot idiosyncratic correlation matrix tends to be more sparse (after removing the latent factors) when the market collapses. In fact, by applying the shrinkage, the sparse structure of both the heat maps in the right panel is apparent, in which majority of the off-diagonal entries are either zero or very close to zero. This supports the low-rank plus sparse assumption on the spot volatility structure.

We next further study the time-varying pattern of the factor number estimation and its connection to the market volatility. The left panel in Figure (ref) plots the averages of the estimated factor numbers over the six pre-determined time intervals in each trading day, and the right panel depicts the realised volatility of the S&P 500 index intraday prices across trading days. It is well known that the S&P 500 index is a market value-weighted index and is thus regarded as an indicator of the entire stock market in US. A virtual comparison between the two plots in Figure (ref) reveals that the averaged factor number tends to be low when the realised volatility of S&P 500 index is high, whereas the estimated factor number is relatively high when the market realised volatility is low. This pattern is not uncommon in the financial market as the extreme fluctuation in the financial market is often closely correlated with some unexpected financial events or economic disruption. Consequently, a small number of factors tend to dominate the market during the period of market collapse. For example, in May 2021, the crypto-currency market collapsed for the first time, and this shock rapidly spread to the entire stock market. This adequately explains the peak of the realised volatility and the low estimated factor number in the midst of May which can be observed in Figure (ref). In contrast, when the financial market is relatively stable, more factors tend to affect the market simultaneously, which is exactly what occurred at the end of May and the beginning of June in 2021.

Conclusion

In this paper we propose a new nonparametric methodology and theory for estimating the low-rank plus sparse spot volatility structure of the noise-contaminated and asynchronous high-frequency data with large dimension. The microstructure noises are allowed to be auto-correlated, nonlinear heteroskedastic, asymptotic vanishing and endogenous, and the latent prices satisfy the time-varying continuous-time factor model. Imposing the sparsity restriction on the spot idiosyncratic volatility matrix, we combine the kernel-weighted POET (or local PCA) and pre-averaging techniques to develop the main estimation method. Under some regularity conditions, we derive the uniform convergence property for the estimated spot volatility matrix in various matrix norms and obtain some explicit convergence rates for different scenarios. The Monte-Carlo simulation study shows that the developed estimation method performs well in finite samples. In particular, the empirical application to the S&P 500 index reveals time-varying patterns of the spot volatility matrix and latent factor number and confirms rationality of assuming the low-rank plus sparse structure.

Acknowledgement

The first author's research was partly supported by the Leverhulme Research Fellowship (RF-2023-396) and BA/Leverhulme Small Research Grant (SRG1920/100603).

Supplementary material

The supplement contains the proofs of the main theoretical results (in Appendix A) and some technical lemmas with proofs (in Appendix B).