Nonparametric Estimation of Large Spot Volatility Matrices for High-Frequency Financial Data

Nonparametric Estimation of Large Spot Volatility Matrices for High-Frequency Financial Data

\centerline{\bf Abstract}

In this paper, we consider estimating spot/instantaneous volatility matrices of high-frequency data collected for a large number of assets. We first combine classic nonparametric kernel-based smoothing with a generalised shrinkage technique in the matrix estimation for noise-free data under a uniform sparsity assumption, a natural extension of the approximate sparsity commonly used in the literature. The uniform consistency property is derived for the proposed spot volatility matrix estimator with convergence rates comparable to the optimal minimax one. For the high-frequency data contaminated by microstructure noise, we introduce a localised pre-averaging estimation method that reduces the effective magnitude of the noise. We then use the estimation tool developed in the noise-free scenario, and derive the uniform convergence rates for the developed spot volatility matrix estimator. We further combine the kernel smoothing with the shrinkage technique to estimate the time-varying volatility matrix of the high-dimensional noise vector. In addition, we consider large spot volatility matrix estimation in time-varying factor models with observable risk factors and derive the uniform convergence property. We provide numerical studies including simulation and empirical application to examine the performance of the proposed estimation methods in finite samples.

\noindentKeywords: \ Brownian semi-martingale, Factor model, Kernel smoothing, Microstructure noise, Sparsity, Spot volatility matrix, Uniform consistency.

Introduction

\setcounter{equation}{0}

Modelling high-frequency financial data is one of the most important topics in financial economics and has received increasing attention in recent decades. Continuous-time econometric models such as the It\^{o} semimartingale are often employed in the high-frequency data analysis. One of the main components in these models is the volatility function or matrix. In the low-dimensional setting (with a single or a small number of assets), the realised volatility is often used to estimate the integrated volatility over a fixed time period AB98, BS02, BS04, ABDL03. In practice, it is not uncommon that the high-frequency financial data are contaminated by the market microstructure noise, which leads to biased realised volatility if the noise is ignored. Hence, various modification techniques such as the two-scale, pre-averaging and realised kernel have been introduced to account for the microstructure noise and produce consistent volatility estimation ZMA05, BHLS08, KL08, JLMPV09, PV09, CKP10, PHL16. S05, ABD10 and AJ14 provide comprehensive reviews for estimating volatility with high-frequency financial data under various settings.

In practical applications, financial economists often have to deal with the situation that there are a large amount of high-frequency financial data collected for a large number of assets. A key issue is to estimate the large volatility structure for these assets, which has applications in various areas such as the optimal portfolio choice and risk management. Partly motivated by developments in large covariance matrix estimation for low-frequency data in the statistical literature, WZ10, TWZ13 and KWZ16 estimate the large volatility matrix under an approximate sparsity assumption BL08; ZL11 and XZ18 study large volatility matrix estimation using the large-dimensional random matrix theory BS10; and LF18 propose a nonparametric eigenvalue-regularised integrated covariance matrix for high-dimensional asset returns. Given that there often exists co-movement between a large number of assets and the co-movement is driven by some risk factors which can be either observable or latent, FFX16, AX17, DLX19 extend the methodologies developed by FLM11, FLM13 to estimate the large volatility matrix by imposing a continuous-time factor model structure on the high-dimensional and high-frequency financial data, and AX19 study the principal component analysis of high-frequency data and derive the asymptotic distribution for the realised eigenvalues, eigenvectors and principal components.

The estimation methodologies in the aforementioned literature often rely on the realised volatility (or covariance) matrices, measuring the integrated volatility structure over a fixed time interval. In practice, it is often interesting to further explore the actual spot/instantaneous volatility structure and its dynamic change over a certain time interval, which is a particularly important measurement for the financial assets when the market is in a volatile period (say, the global financial crisis or COVID-19 outbreak). For a single financial asset, FW08 and K10 introduce a kernel-based nonparametric method to estimate the spot volatility function and establish its asymptotic properties including the point-wise and global asymptotic distribution theory and uniform consistency. For the noise-contaminated high-frequency data, ZB14 combine the two-scale realised volatility with the kernel-weighted technique to estimate the spot volatility, whereas KK16 propose a kernel-weighted pre-averaging spot volatility estimation method. Other nonparametric spot volatility estimation methods can be found in FFL07 and FL20. It seems straightforward to extend this local nonparametric method to estimate the spot volatility matrix for a small number of assets. However, a further extension to the setting with vast financial assets is non-trivial. There is virtually no work on estimating the large spot volatility matrix except K18, which considers estimating large spot volatility matrices and their integrated versions under the continuous-time factor model structure for noise-free high-frequency data.

The main methodological and theoretical contributions of this paper are summarised as follows.

We argue that all three of the scenarios we consider above may be practically relevant. Microstructure noise is considered important in the very highest frequency of data, whereas researchers working with five minute data, say, often ignore the noise. For the lower frequency of data there is a lot comovement in returns and the factor model is designed to capture that comovement, whereas at the ultra high frequency, comovement is less of an issue; indeed, under the so-called Epps effect this comovement shrinks to zero with sampling frequency.

The rest of the paper is organised as follows. In Section (ref), we estimate the large spot volatility matrix in the noise-free high-frequency data setting and give the uniform consistency property. In Section (ref) , we extend the methodology and theory to the noise-contaminated data setting and further estimate the time-varying noise volatility matrix. Section (ref) considers the large spot volatility matrix with systematic factors. Section (ref) reports the simulation studies and Section (ref) provides an empirical application. Section (ref) concludes the paper. Proofs of the main theoretical results are available in Appendix A. The supplementary document contains proofs of some technical lemmas and propositions and discussions on the spot precision matrix estimation and the asynchronicity issue. Throughout the paper, we let $\Vert\cdot\Vert_2$ be the Euclidean norm of a vector; and for a $d\times d$ matrix ${\mathbf{A}}=(A_{ij})_{d\times d}$, we let $\Vert { \mathbf{A}}\Vert$ and $\Vert {\mathbf A}\Vert_F$ be the matrix spectral norm and Frobenius norm, $\vert{\mathbf{A}} \vert_1=\sum_{i=1}^d\sum_{i=1}^d |A_{ij}|$, $\Vert{\mathbf{A}} \Vert_1=\max_{1\leq j\leq d}\sum_{i=1}^d |A_{ij}|$, $\Vert{\mathbf{A}} \Vert_{\infty,q}=\max_{1\leq i\leq d}\sum_{j=1}^d |A_{ij}|^q$ and $\Vert { \mathbf{A}}\Vert_{\max}=\max_{1\leq i\leq d}\max_{1\leq j\leq d} |A_{ij}|$.

Estimation with noise-free data

\setcounter{equation}{0}

Suppose that ${\mathbf{X}}_t=\left(X_{1,t},\cdots,X_{p,t}\right)^{^\intercal} $ is a $p$-variate Brownian semi-martingale solving the following stochastic differential equation:

where ${\mathbf{W}}_t=\left(W_{1,t},\cdots,W_{p,t}\right)^{^\intercal}$ is a $p$-dimensional standard Brownian motion, ${\boldsymbol{\mu}} _t=(\mu_{1,t},\cdots,\mu_{p,t})^{^\intercal}$ is a $p$-dimensional drift vector, and ${\boldsymbol{\sigma}}_t=\left(\sigma_{ij,t}\right)_{p\times p}$ is a $p\times p$ matrix. The spot volatility matrix of ${\mathbf{X}}_t$ is defined as

Our main interest lies in estimating ${\boldsymbol{\Sigma}}_t$ when the size $p$ is large. As in CXW13 and CL16, we assume that the true spot volatility matrix satisfies the following uniform sparsity condition: $ \left\{{\boldsymbol{\Sigma}}_t:\ 0\leq t\leq T\right\}\in \mathcal{S} (q,\varpi(p), T)$, where

where $0\le q<1$, $\varpi(p)$ is larger than a positive constant, $T$ is a fixed positive number and $\Lambda$ is a positive random variable satisfying $\mathsf{E}[\Lambda]\leq C_\Lambda<\infty$. This is a natural extension of the approximate sparsity assumption BL08. Section (ref) below will relax this assumption and consider estimating large spot volatility matrices with systematic factors. The asset prices are assumed to be collected over a fixed time interval $[0,T]$ at $0,\Delta,2\Delta,\cdots,n\Delta$, where $\Delta$ is the sampling frequency and $n=\lfloor T/\Delta\rfloor$ with $\lfloor\cdot\rfloor$ denoting the floor function. In the main text, we focus on the case of equidistant time points in the high-frequency data collection. The asynchronicity issue will be discussed in Appendix C.2 of the supplement.

For each $1\leq i,j\leq p$, we estimate the spot co-volatility $\Sigma_{ij,t} $ by

with \[ K_{h}^\ast(t_k-t)=K_h\left(t_k-t\right)/\left[\Delta\sum_{l=1}^nK_h\left(t_l-t\right)\right], \] where $t_k=k\Delta$, $K_h(u)=h^{-1}K(u/h)$, $K(\cdot)$ is a kernel function, $h$ is a bandwidth shrinking to zero and $\Delta X_{i,k}=X_{i,t_k}-X_{i,t_{k-1}}$. The use of $K_h^\ast(t_k-t)$ rather than $K_h(t_k-t)$ in the estimation ((ref)) is to correct a constant bias when $t$ is close to the boundary points $0$ and $T$. A naive method of estimating the spot volatility matrix ${\boldsymbol{\Sigma}}_{t}$ is to directly use $\widehat\Sigma_{ij,t}$ to form an estimated matrix. However, this estimate often performs poorly in practice when the number of assets is very large (say, $p>n$). To address this issue, a commonly-used technique is to apply a shrinkage function to $\widehat\Sigma_{ij,t}$ when $i\neq j$, forcing very small estimated off-diagonal entries to be zeros. Let $ s_\rho(\cdot)$ denote a shrinkage function satisfying the following three conditions: (i) $\vert s_\rho(u)\vert \leq \vert u\vert$ for $u\in\mathscr{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. With the shrinkage function, we construct the following nonparametric estimator of ${\boldsymbol{\Sigma}}_{t}$:

where $\rho_1(t)$ is a tuning parameter which is allowed to change over $t$ and $I(\cdot)$ denotes the indicator function. Section (ref) discusses the choice of $\rho_1(t)$, ensuring that $\widehat{\boldsymbol\Sigma}_t$ is positive definite in finite samples. Our estimation method of the spot volatility matrix can be seen as a natural extension of the kernel-based large sparse covariance matrix estimation CXW13, CL16,CLL19 from the low-frequency data setting to the high-frequency one. We next give some technical assumptions which are needed to derive the uniform convergence property of $\widehat{\boldsymbol{ \Sigma}}_{t}$.

The following theorem gives the uniform convergence property (in the matrix spectral norm) for the spot volatility matrix estimator $\widehat{\boldsymbol{\Sigma}}_{t}$ under the uniform sparsity assumption.

Estimation with contaminated high-frequency data

\setcounter{equation}{0}

In practice, it is not uncommon that high-frequency financial data are contaminated by the market microstructure noise. The kernel estimation method proposed in Section (ref) would be biased if the noise is ignored in the estimation procedure. Consider the following additive noise structure:

where $t_k=k\Delta$, $k=1,\cdots,n$, ${\mathbf{Z}}_{t}=\left(Z_{1,t}, \cdots,Z_{p,t}\right)^{^\intercal}$ is a vector of observed asset prices at time $t$, and ${\boldsymbol{\xi}}_k=(\xi_{1,k},\cdots,\xi_{p,k})^{^\intercal} $ is a $p$-dimensional vector of noises with nonlinear heteroskedasticity, ${ \boldsymbol{\omega}}(\cdot)=\left[\omega_{ij}(\cdot)\right]_{p\times p}$ is a $p\times p$ matrix of deterministic functions, and ${\boldsymbol{\xi}} _k^\ast=\left(\xi_{1,k}^\ast,\cdots,\xi_{p,k}^\ast\right)^{^\intercal}$ independently follows a $p$-variate identical distribution. The noise structure defined in ((ref)) is similar to the setting considered in KL08 which also contains a nonlinear mean function and allows the existence of endogeneity for a single asset. Throughout this section, we assume that $\{{\boldsymbol{\xi}}_k^\ast\}$ is independent of the Brownian semimartingale $\{{\mathbf{X}}_t\}$.

Estimation of the spot volatility matrix

To account for the microstructure noise and produce consistent volatility matrix estimation, we apply the pre-averaging technique as the realised kernel estimate BHLS08 can be seen as a member of the pre-averaging estimation class whereas the two-scale estimate ZMA05 can be re-written as the realised kernel estimate with the Bartlett-type kernel (up to the first-order approximation). The pre-averaging method has been studied by JLMPV09, PV09 and CKP10 in estimating the integrated volatility for a single asset and is further extended by KWZ16 and DLX19 to the large high-frequency data setting. KK16 use a localised pre-averaging technique to estimate the spot volatility function for a single asset and derive the uniform convergence rate for the developed estimate. A similar technique is also used by XL02 to improve convergence of the nonparametric spectral density estimator for time series with general autocorrelation for low-frequency data.

We first pre-average the observed high-frequency data via a kernel filter, i.e.,

with $L_{b}^\dagger(t_k-\tau)=L_b\left(t_k-\tau\right)/\int_0^TL_b(s-\tau)ds$, where $L_b(u)=b^{-1}L(u/b)$, $L(\cdot)$ is a kernel function and $b$ is a bandwidth. Let $\Delta \widetilde{X}_{i,l}=\widetilde{X}_{i,\tau_l}-\widetilde{X}_{i,\tau_{l-1}}$, where $\widetilde{X}_{i,\tau_l}$ is the $i$-th component of $\widetilde{ \mathbf{X}}_{\tau_l}$ and $\tau_0,\tau_1,\cdots,\tau_N$ are the pseudo-sampling time points in the fixed interval $[0,T]$ with equal distance $\Delta_\ast=T/N$. Replacing $\Delta X_{i,k}$ by $\Delta\widetilde{X }_{i,l}$ in ((ref)), we estimate the spot co-volatility $\Sigma_{ij,t}$ by

where \[ K_{h}^\dagger(\tau_l-t)=K_h\left(\tau_l-t\right)/\left[\Delta_\ast\sum_{k=1}^NK_h\left(\tau_k-t\right)\right]. \] Furthermore, to obtain a stable spot volatility matrix estimate in finite samples when the dimension $p$ is large, as in ((ref)), we apply shrinkage to $\widetilde{\Sigma}_{ij,t}$, $1\leq i\neq j\leq p$, and subsequently construct

where $\rho_2(t)$ is another time-varying shrinkage parameter. We next give some conditions needed to derive the uniform consistency property of $ \widetilde{\boldsymbol{\Sigma}}_t$.

Estimation of the time-varying noise volatility matrix

It is often interesting to further explore the volatility structure of microstructure noise. CHLT21 estimate the constant covariance matrix for high-dimensional noise and derive the optimal convergence rates for the developed estimate. In the present paper, we consider the time-varying noise covariance matrix defined by

It is sensible to assume that $\left\{{\boldsymbol{\Omega}}(t):\ 0\leq t\leq T\right\}$ satisfies the uniform sparsity condition as in ((ref)). For each $1\leq i,j\leq p$, we estimate $\Omega_{ij}(t)$ by the kernel smoothing method:

where $h_1$ is a bandwidth, $\Delta Z_{i,t_{k}}=Z_{i,t_k}-Z_{i,t_{k-1}}$ and $K_{h_1}^\ast(t_k-t)$ is defined similarly to $K_{h}^\ast(t_k-t)$ in ((ref)) but with $h_1$ replacing $h$. As in ((ref)) and ((ref)), we again apply shrinkage to $ \widehat{\Omega}_{ij}(t)$, $1\leq i\neq j\leq p$, and construct

where $\rho_3(t)$ is a time-varying shrinkage parameter. To derive the uniform consistency property of $\widehat{\boldsymbol{\Omega}}(t)$, we need to impose stronger moment condition on ${\boldsymbol{\xi}}_k^\ast$ and smoothness restriction on $\Omega_{ij}(\cdot)$.

In the following theorem, we state the uniform consistency result for $ \widehat{\boldsymbol{\Omega}}(t)$ with convergence rate comparable to that in Theorem (ref).

Estimation with observed factors

\setcounter{equation}{0}

The large spot volatility matrix estimation with the shrinkage technique developed in Sections (ref) and (ref) heavily relies on the uniform sparsity assumption ((ref)). However, the latter may be too restrictive in practice since the price processes of a large number of assets are often driven by some common factors such as the market factors, resulting in strong correlation among assets and failure of the sparsity condition. To address this problem, we next consider the nonparametric time-varying regression at high frequency:

where ${\boldsymbol\beta}(t)=\left[\beta_{1}(t),\cdots,\beta_p(t)\right]^{^\intercal}$ is a $p\times k$ matrix of time-varying betas (or factor loadings), ${\mathbf F}_t$ and ${\mathbf X}_t$ are $k$-variate and $p$-variate continuous semi-martingales defined by

respectively, ${\boldsymbol{\mu}}_t^F$ and ${\boldsymbol{\mu}}_t^X$ are drift vectors, ${\boldsymbol{\sigma}}_t^F=\left(\sigma_{ij,t}^F\right)_{k\times k}$, ${\boldsymbol{\sigma}}_t^X=\left(\sigma_{ij,t}^X\right)_{p\times p}$, ${\mathbf{W}}_t^F$ and ${\mathbf{W}}_t^X$ are $k$-dimensional and $p$-dimensional standard Brownian motions. For the time being, we assume that ${\mathbf Y}_t$ and ${\mathbf F}_t$ are observable and noise free but ${\mathbf X}_t$ is latent. Extension of the methodology and theory to the noise-contaminated high-frequency data will be considered later in this section.

Estimation of the constant betas via the ratio of realised covariance to realised variance is proposed by BS04, and extension to time-varying beta estimation has been studied by MZ06, RTT15 and AKX20, some of which allow jumps in the semi-martingale processes. The main interest of this section lies in estimating the large spot volatility structure ${\boldsymbol\Sigma}_t^Y$ of ${\mathbf Y}_t$. Letting ${\boldsymbol\Sigma}_t^F={\boldsymbol\sigma}_t^F\left({\boldsymbol\sigma}_t^F\right)^{^\intercal}$ and ${\boldsymbol\Sigma}_t^X={\boldsymbol\sigma}_t^X\left({\boldsymbol\sigma}_t^X\right)^{^\intercal}$, and assuming orthogonality between ${\mathbf X}_t$ and ${\mathbf F}_t$, see Assumption (ref)(iii) below, it follows from ((ref)) that

As in FLM11, FLM13, we impose the uniform sparsity restriction on ${\boldsymbol\Sigma}_t^X$ instead of ${\boldsymbol\Sigma}_t^Y$, i.e., $ \left\{{\boldsymbol{\Sigma}}_t^X:\ 0\leq t\leq T\right\}\in \mathcal{S}(q,\varpi(p), T)$. This is a reasonable assumption in practical applications as the asset prices, after removing the influence of systematic factors, are expected to be weakly correlated. FFX16 and DLX19 use a similar framework with constant betas to estimate large integrated volatility matrices.

Suppose that we observe ${\mathbf Y}_t$ and ${\mathbf F}_t$ at regular points: $t_k=k\Delta$, $k=1,\cdots,n$, as in Sections (ref) and (ref). Let ${\boldsymbol\Sigma}_t^{YF}$ be the spot covariance between ${\mathbf Y}_t$ and ${\mathbf F}_t$. We may use the kernel smoothing method as in ((ref)) to estimate ${\boldsymbol\Sigma}_t^{Y}$, ${\boldsymbol\Sigma}_t^{F}$ and ${\boldsymbol\Sigma}_t^{YF}$, i.e.,

where $\Delta{\mathbf Y}_k={\mathbf Y}_{t_k}-{\mathbf Y}_{t_{k-1}}$, $\Delta{\mathbf F}_k={\mathbf F}_{t_k}-{\mathbf F}_{t_{k-1}}$, and $K_h^\ast(t_k-t)$ is defined as in ((ref)). Consequently, the time-varying betas ${\boldsymbol\beta}(t)$ and the spot idiosyncratic volatility matrix ${\boldsymbol\Sigma}_t^X$ are estimated by

and

With the uniform sparsity condition, it is sensible to further apply shrinkage to $\widehat\Sigma_{ij,t}^X$, i.e.,

where $\rho_4(t)$ is a time-varying shrinkage parameter. We finally estimate ${\boldsymbol\Sigma}_t^Y$ as

We need the following assumption to derive the uniform convergence property for $\widehat{\boldsymbol\Sigma}_t^{X,s}$ and $\widehat{\boldsymbol\Sigma}_t^{Y,s}$.

We next present the convergence property of $\widehat{\boldsymbol\Sigma}_t^{X,s}$ and $\widehat{\boldsymbol\Sigma}_t^{Y,s}$ defined in ((ref)) and ((ref)), respectively. Due to the nonparametric factor regression model structure ((ref)), the largest $k$ eigenvalues of ${\boldsymbol\Sigma}_t^Y$ are spiked, diverging at a rate of $p$. Hence, ${\boldsymbol\Sigma}_t^Y$ cannot be consistently estimated in the absolute term. To address this problem, as in FLM11, FLM13, we measure the spiked volatility matrix estimate in the following relative error: \[\left\Vert \widehat{\boldsymbol\Sigma}_t^{Y,s}-{\boldsymbol\Sigma}_t^Y\right\Vert_{{\boldsymbol\Sigma}_t^Y}=\frac{1}{\sqrt{p}}\left\Vert \left({\boldsymbol\Sigma}_t^Y\right)^{-1/2}\left(\widehat{\boldsymbol\Sigma}_t^{Y,s}-{\boldsymbol\Sigma}_t^Y\right) \left({\boldsymbol\Sigma}_t^Y\right)^{-1/2} \right\Vert_F,\] where the normalisation factor $p^{-1/2}$ is used to guarantee that $\left\Vert{\boldsymbol\Sigma}_t^Y\right\Vert_{{\boldsymbol\Sigma}_t^Y}=1$.

We next modify the above methodology and theory to accommodate microstructure noise in the asset prices and factors. Assume that

where ${\boldsymbol{\omega}}_Y(\cdot)$ and ${\boldsymbol{\omega}}_F(\cdot)$ are matrices of deterministic functions similar to ${\boldsymbol{\omega}}(\cdot)$, and $\{{\boldsymbol{\xi}}_{Y,k}^\ast\}$ and $\{{\boldsymbol{\xi}}_{F,k}^\ast\}$ are i.i.d. sequences of random vectors similar to $\{{\boldsymbol{\xi}}_{k}^\ast\}$. Since both ${\mathbf Y}_t$ and ${\mathbf F}_t$ are latent, we need to first adopt the kernel pre-averaging technique proposed in Section (ref) to obtain the approximation of ${\mathbf Y}_t$ and ${\mathbf F}_t$, and then apply the kernel smoothing and generalised shrinkage as in ((ref))--((ref)). This results in a three-stage estimation procedure which we describe as follows.

As shown in Theorem (ref), the existence of microstructure noises slows down the uniform convergence rates. Following the proof of Lemma B.1 in Appendix B, we may show that \[ \max_{0\leq l\leq N}\left\vert \widetilde{\mathbf{Y}}_{\tau_l}-{\mathbf{Y}}_{\tau_l}\right\vert_{\max}+\max_{0\leq l\leq N}\left\vert \widetilde{\mathbf{F}}_{\tau_l}-{\mathbf{F}}_{\tau_l}\right\vert_{\max}=O_P\left(\nu_{\Delta,p,N}\right), \] where $\vert\cdot\vert_{\max}$ denotes the $L_\infty$-norm of a vector, and $\nu_{\Delta,p,N}$ is defined in Assumption (ref)(iii). Modifying Proposition (ref) and the proof of Theorem (ref) in Appendix A, we can prove that ((ref)) and ((ref)) hold but with $\zeta_{\Delta,p}$ replaced by $\zeta_{N,p}^\ast+\nu_{\Delta,p,N}$ defined in Assumption (ref)(iii), i.e.,

Monte-Carlo Study

\setcounter{equation}{0}

In this section, we report the Monte-Carlo simulation studies to assess the numerical performance of the proposed large spot volatility matrix and time-varying noise volatility matrix estimation methods under the sparsity condition and the factor-based spot volatility matrix estimation. Here we only consider the synchronous high-frequency data. Additional simulation results for asynchronous high-frequency data are provided in the supplement.

Simulation for sparse volatility matrix estimation

{\em 5.1.1.\ \ Simulation setup}

We generate the noise-contaminated high-frequency data according to model ((ref)), where ${\boldsymbol{\omega}}(t)$ is taken as the Cholesky decomposition of the noise covariance matrix ${\boldsymbol\Omega}(t)=\left[\Omega _{ij}(t)\right] _{p\times p}$, ${\boldsymbol{\xi }}_{k}^{\ast }=\left(\xi _{1,k}^{\ast },\cdots ,\xi _{p,k}^{\ast }\right)^{^{\intercal }}$ is an independent $p$-dimensional random vector of cross-sectionally independent standard normal random variables, the latent return process $\mathbf{X}_{t}$ of $p$ assets is generated from the following drift-free model:

$\mathbf{W}_{t}^{X}=\left( W_{1,t}^{X},\cdots ,W_{p,t}^{X}\right)^{^{\intercal }}$ is a standard $p$-dimensional Brownian motion, and ${\boldsymbol{\sigma }}_{t}$ is chosen as the Cholesky decomposition of the spot covariance matrix ${\boldsymbol{\Sigma }}_{t}=\left(\Sigma_{ij,t}\right) _{p\times p}$. In the simulation, we consider the volatility matrix estimation over the time interval of a full trading day, and set the sampling interval to be $15$ seconds, i.e., $\Delta =1/(252\times 6.5\times 60\times 4)$, to generate synchronous data. We consider three structures in ${\boldsymbol{\Sigma}}_{t}$ and ${\boldsymbol\Omega}(t)$: “banding", “block-diagonal", and “exponentially decaying". Following WZ10, we generate the diagonal elements of ${\boldsymbol\Sigma}_{t}$ from the following geometric Ornstein-Uhlenbeck model BS02:

where $\mathbf{W}_{t}^\ast=\left(W_{1,t}^\ast,\cdots,W_{p,t}^\ast\right)^{^\intercal}$ is a standard $p$-dimensional Brownian motion independent of $\mathbf{W}_{t}^{X}$, and $\iota_{i}$ is a random number generated uniformly between $-0.62$ and $-0.30$, reflecting the leverage effects. The diagonal elements of ${\boldsymbol\Omega}(t)$ are defined as daily cyclical deterministic functions of time: \[ \Omega _{ii}\left( t\right) =c_{i}\left\{ \frac{1}{2}\left[ \cos \left(2\pi t/T\right) +1\right] \times \left( \overline\omega -\underline\omega\right) +\underline\omega\right\}, \] where $\overline\omega=1$ and $\underline\omega=0.1$ reflect the observation by KL08 that the noise level is high at both the opening and the closing times of a trading day and is low in the middle of the day, and the scalar $c_{i}$ controls the noise ratio for each asset which is chosen to match the highest noise ratio considered by WZ10. As in BS02, BS04, we define a continuous-time stochastic process $\kappa^{\Sigma}_t$ by

where $W_{t}^\diamond$ is a standard univariate Brownian motion independent of $\mathbf{W}_{t}^{X}$ and $\mathbf{W}_{t}^\ast$. Let \[ \kappa^\Omega_t =\frac{\overline\kappa-\underline\kappa}{2}\left[ \cos \left( 2\pi t/T\right) +1\right] +\underline\kappa, \] where $\overline\kappa=0.5$ and $\underline\kappa=-0.5$. We will use $\kappa^{\Sigma}_t$ and $\kappa^{\Omega}_t$ to define the off-diagonal elements in ${\boldsymbol{\Sigma}}_{t}$ and ${\boldsymbol\Omega}(t)$, respectively, which are specified as follows.

It is clear that the sparsity condition is not satisfied when the off-diagonal elements of ${\boldsymbol{\Sigma}}_{t}$ and ${\boldsymbol\Omega}(t)$ are exponentially decaying as in ((ref)). The number of assets $p$ is set as $p=200$ and $500$ and the replication number is $R=200$

{\em 5.1.2.\ \ Volatility matrix estimation}

In the simulation studies, we consider the following volatility matrix estimates.

The choice of tuning parameter in shrinkage is similar to that in DLX19. For example, in the noise-free spot volatility estimate, we set the tuning parameter as $\rho_{ij}(t)=\rho(t)(\widehat{\Sigma}_{ii,t}\widehat{\Sigma}_{jj,t})^{1/2}$ where $\rho(t)$ is chosen as the minimum value among the grid of values on $[0,1]$ such that the shrinkage estimate of the spot volatility matrix is positive definite. To evaluate the estimation performance of $\widehat{\boldsymbol\Sigma }_{t}$, we consider $21$ equidistant time points on $[0,T]$ and compute the following Mean Frobenius Loss (MFL) and Mean Spectral Loss (MSL) over $200$ repetitions:

where $t_{j}$, $j=1,2,\cdots,21$ are the equidistant time points on the interval $[0,T]$, and $\boldsymbol{\widehat{\Sigma}}_{t_{j}}^{(m)}$ and $\boldsymbol{\Sigma }_{t_{j}}^{(m)}$ are respectively the estimated and true spot volatility matrices at $t_{j}$ for the $m$-th repetition. The “MFL" and “MSL" can be similarly defined for $\widetilde{\boldsymbol\Sigma}_{t}$ and $\widehat{\boldsymbol\Omega}(t)$.

{\em 5.1.3.\ \ Simulation results}

Table 1 reports the simulation results when the dimension is $p=200$. The three panels in the table (from top to bottom) report the results where the true volatility matrix structures are banding, block-diagonal, and exponentially decaying, respectively. In each panel, the MFL results are reported on the left, whereas the MSL results are on the right. The first two rows of each panel contain the MFL and MSL results for the spot volatility matrix estimation whereas the third row contains the results for the time-varying noise volatility matrix estimation.

For the noise-free estimate $\widehat{\boldsymbol\Sigma }_{t}$, when the volatility matrix structure is banding, the performance of the four shrinkage estimators are substantially better than that of the naive estimate (without any shrinkage). In particular, the results of the soft thresholding, adaptive LASSO and SCAD are very similar and their MFL and MSL values are approximately one third of those of the naive estimator. Meanwhile, the performance of the hard thresholding is less accurate (despite the much stronger level of shrinking used), but is still much better than the naive estimate. These results show that the shrinkage technique is an effective tool in estimating the sparse volatility matrix. Similar results are obtained for the noise-contaminated estimate $\widetilde{\boldsymbol\Sigma}_{t}$. Unsurprisingly, due to the microstructure noise, the MFL and MSL values of the local pre-averaging estimates are noticeably higher than the corresponding values of the noise-free estimates. We next turn the attention to the time-varying noise volatility matrix estimate $\widehat{\boldsymbol\Omega}(t)$. As in the spot volatility matrix estimation, the naive method again produces the highest MFL and MSL values. The performance of the four shrinkage estimators are similar with the adaptive LASSO and SCAD being slightly better than the hard and soft thresholding. The simulation results for the block-diagonal and exponentially decaying covariance matrix settings, reported in the middle and bottom panels of Table 1, are fairly close to those for the banding setting. Overall, the results in Table 1 show that the shrinkage methods perform well not only in the sparse covariance matrix settings but also in the non-sparse one (i.e., the exponentially decaying setting).

{\scriptsize The selected bandwidths are $h^\ast=90$ for $\widehat{\boldsymbol\Sigma}_{t}$, $h^\ast=90$ and $b^\ast=4$ for $\widetilde{\boldsymbol\Sigma}_{t}$, and $h_1^{\ast}=90$ for $\widehat{\boldsymbol\Omega}(t)$, where $h^{\ast}=h/\Delta$, $b^{\ast}=b/\Delta$, and $h_{1}^{\ast}=h_1/\Delta$.}

The simulation results when the dimension is $p=500$ are reported in Table 2. Overall the results are very similar to those in Table 1, so we omit the detailed discussion and comparison to save the space.

{\scriptsize The selected bandwidths are $h^{\ast }=240$ for $\widehat{\boldsymbol\Sigma}_{t}$, $h^{\ast }=240$, $b^{\ast }=4$ for $\widetilde{\boldsymbol\Sigma}_{t}$ and $h_1^{\ast }=240$ for $\widehat{\boldsymbol\Omega}(t)$, where $h^{\ast}=h/\Delta$, $b^{\ast}=b/\Delta$, and $h_{1}^{\ast}=h_1/\Delta$.}

Simulation for factor-based spot volatility matrix estimation

{\em 5.2.1.\ \ Simulation setup}

We generate ${\mathbf Y}_t$ via ((ref)), where the $p$-dimensional idiosyncratic returns follow the dynamics of $d\mathbf{X}_{t}$ defined in ((ref)). In this simulation, we only consider $p=500$. As in AKX20, we adopt a three-factor model, where the factors $\mathbf{F}_{t}=\left(F_{1,t},F_{2,t},F_{3,t}\right) ^{^\intercal}$ are generated by

The factor volatilities are driven by

where ${\sf E}[ dW_{k,t}^{F}d\widetilde{W}_{k,t}] =\rho_{k}dt$, allowing for potential leverage effects in the factor dynamics. Both $W_{k,t}^{F}$ and $\widetilde{W}_{k,t}$ are standard univariate Brownian motions. In the simulation, we set $\left( \tilde{\kappa}_{1},\tilde{\kappa}_{2},\tilde{ \kappa}_{3}\right) =\left( 3,4,5\right) $, $\left( \tilde{\alpha}_{1},\tilde{ \alpha}_{2},\tilde{\alpha}_{3}\right) =\left( 0.09,0.04,0.06\right) $, $ \left( \tilde{\nu}_{1},\tilde{\nu}_{2},\tilde{\nu}_{3}\right) =\left( 0.3,0.4,0.3\right) $, $\left( \mu _{1}^{F},\mu _{2}^{F},\mu _{3}^{F}\right) =\left( 0.05,0.03,0.02\right) $, $\left( \rho _{1},\rho _{2},\rho _{3}\right) =\left( -0.6,-0.4,-0.25\right) $ and $\left( \rho _{12},\rho _{13},\rho _{23}\right) =\left( 0.05,0.10,0.15\right) $.

We consider the following three cases for generating the time-varying beta processes: $\beta_{i}(t)=\left[\beta_{i,1}(t), \beta_{i,2}(t), \beta_{i,3}(t)\right]^{^\intercal}$, $i=1,\cdots,p$.

{\em 5.2.2.\ \ Simulation Results}

The spot idiosyncratic volatility matrix is estimated via ((ref)). For ease of comparison, we use exactly the same bandwidth as in our first experiment. The results for the noise-free and noise-contaminated spot idiosyncratic volatility matrix estimates $\widehat{\boldsymbol{\Sigma }}_{t}$ and $\widetilde{\boldsymbol{\Sigma }}_{t}$ measured by MFL and MSL are reported in Table 3, which reveal some desirable observations. Firstly, we note that our estimation results in terms of MFL and MSL are almost identical across different types of dynamics of factor loadings, indicating that the developed estimation procedure is robust in finite samples to different assumption of the factor loading dynamics as long as they satisfy our smooth restriction, see Assumption (ref)(iv). Secondly, the MFL and MSL values are similar to those reported in Table 2 which were obtained based on data generating model without common factors. This means that the proposed nonparametric time-varying high frequency regression can effectively remove common factors, resulting in accurate estimation of the spot idiosyncratic volatility matrix.

The factor-based spot volatility matrix of ${\mathbf Y}_t$ is estimated via ((ref)). As discussed in Section (ref), we measure the accuracy of the spiked volatility matrix estimate by the relative error defined above Theorem (ref), i.e., consider the following Mean Relative Loss (MRL): \[ \text{MRL}=\frac{1}{200}\sum_{m=1}^{200}\left( \frac{1}{21} \sum_{j=1}^{21}\left\Vert \widehat{\boldsymbol{\Sigma }}_{t_{j}}^{Y,(m)}- \boldsymbol{\Sigma }_{t_{j}}^{Y,(m)}\right\Vert _{\boldsymbol{\Sigma } _{t_{j}}^{Y,(m)}}\right) . \] The relevant results are reported in Table 4, where $\widehat{\boldsymbol\Sigma }_{t}^{Y}$ and $\widetilde{\boldsymbol\Sigma }_{t}^{Y}$ denote the noise-free and noise-contaminated factor-based spot volatility matrix estimates, respectively. We can see that the performance of the shrinkage estimates is substantially better than that of the naive estimate. Unsurprisingly, due to the presence of microstructure noise, the MRL results of $\widetilde{\boldsymbol\Sigma }_{t}^{Y}$ are much higher than those of $\widehat{\boldsymbol\Sigma }_{t}^{Y}$. As in Table 3, our proposed estimation is robust to different factor loading dynamics.

Empirical Study

\setcounter{equation}{0}

We apply the proposed methods to the intraday returns of the S&P 500 component stocks to demonstrate the effectiveness of our nonparametric spot volatility matrix estimation in revealing time-varying patterns. We consider the 5-minute returns of the S&P 500 stocks collected in September 2008. On September 15 Lehman Brothers filed for bankruptcy, causing shockwaves throughout the global financial system. Hence, it is interesting to examine how the spot volatility structure of the returns evolved during this one-month period. In addition, to demonstrate the effectiveness of our model with observed risk factors in explaining the systemic component of the dependence structure, we also collect the 5-minute returns of twelve factors. The first three factors are constructed in AKX20 as our proxy for the market (MKT), small-minus-big market capitalisation (SMB), and high-minus-low price-earning ratio (HML). The other nine factors are the widely available sector SDPR ETFs, which are intended to tract the following nine largest S&P sectors: Energy (XLE), Materials (XLB), Industrials (XLI), Consumer Discretionary (XLY), Consumer Staples (XLP), Health Care (XLV), Financial (XLF), Information Technology (XLK), Utilities (XLU). We sort our stocks according to their GICS (Global Industry Classification Standard) codes, so that they are grouped by sectors in the above order. Consequently, the correlation (sub)matrix for stocks within each sector corresponds to a block on the diagonal of the full correlation matrix FFX16.

We only use stocks that are included in the S&P 500 index and whose GICS codes are unchanged in September 2008. We also exclude stocks that do not belong to any of the above nine sectors. This leaves us with a total of $p=482$ stocks. All the returns are synchronised via the previous-tick subsampling technique Zh11, and overnight returns are removed because of potential dividends and stock splits. Consequently, we have $1638$ time series observations for each of the $482$ stocks. For the 5-minute returns, we may assume that the potential impact of microstructure noises are negligible. The smoothing parameter in our kernel estimation is chosen as $h=2/252$ (equivalent to 2 trading days)\footnote{We experimented three bandwidth choices, namely, 1 day, 2 days and 3 days. We found that $h=1/252$ ($h=3/252$) produced clearly undersmoothed (oversmoothed) time series of estimated deciles of the cross sectional distribution of the variances and pairwise correlations of our returns, whereas $h=2/252$ seems to be most reasonable. Our qualitative conclusion is unaffected by the choices of $h$ within the range of 1 to 3 days.}.

We start with estimating the spot volatility matrices of the total returns (i.e., the observed returns) without incorporating the observed factors or applying any shrinkage. To visualise the potential time variation of the estimated spot matrices, as in BHMR19, we plot the time series of deciles of the distribution of the estimated spot variances and the pairwise correlations\footnote{The nine decile levels we use in this study are the 10th, 20th,..., and 90th percentiles.}. The patterns of the spot variances and correlations in Figure (ref)(a) and (b) reveal some clear evidence of time variation in our sampling period. We note that the distributions of the variances are relatively narrow and stay low on the first few days of the month. However, close to Lehman Brothers' announcement on the 15th, they start to rise and get wider quite rapidly and reach the peak around the 17th and the 18th. The spot variances at the peak are much higher than those on the earlier days of the month. The distributions return to the earlier level in the following week. In contrast, the distributions of pairwise spot correlations also start to shift up around the same time, but quickly reach the peak on the 16th (only one day after the bankruptcy news), and then dip to a relatively low point around the 19th before returning to the earlier level. Such time-varying features in the dynamics covariance structure are quite interesting and sensible, reflecting the impact of market news. Hence, our proposed spot volatility matrix estimation methodology provides a useful tool for revealing such dynamics.

To examine whether it is appropriate to directly apply shrinkage techniques to the spot volatility matrices of the total returns, following FFX16 we plot in Figure (ref)(a) and (b) their sparsity patterns on the 16th and the 19th of September \footnote{Recall that as in DLX19 our tuning parameter used for each pairwise spot covariance in our shrinkage method is proportional to the product of the spot standard deviations of the returns of that pair of assets. Therefore, the sparsity pattern is effectively determined by the spot correlation matrix.}. The deep blue dots correspond to the locations of pairwise correlations that are at least $0.15$, whereas the white dots correspond to those smaller than $0.15$. Note that the covariance structure of the total returns is very dense on these two days. Therefore, it is not appropriate to directly apply the shrinkage technique as in Sections (ref) and (ref). Meanwhile, although both are quite dense, we can still clearly see their differences. Consistent with our observation from the decile plots of the correlations, we can see that the plot for the 16th is almost completely covered by blue dots, but the plot for the 19th in contrast has significantly more areas covered in white.

We next incorporate the twelve observed factors in the large spot volatility matrix estimation as suggested in Section (ref). In particular, we are interested on estimating the spot idiosyncratic volatility matrix, which is expected to satisfy the sparsity restriction. To save space, we choose to only report results using the SCAD shrinkage due to its satisfactory performance in our simulations. In Figure (ref)(c) and (d), we plot the deciles of the estimated spot idiosyncratic variances and correlations over trading days. In Figure (ref)(c), we observe a significant upward shift of the distribution of the spot variances of the idiosyncratic returns around the time of Lehman Brothers' bankruptcy, indicating that the observed factors may not fully capture the time variation of the spot variances. In contrast, the deciles of the spot correlations in Figure (ref)(d) seem to be quite flat throughout the entire month, suggesting that the systematic factors may explain the time variation in the distribution of the pairwise correlations better than that of the variances.

We finally plot the sparsity patterns of the two estimated spot idiosyncratic volatility matrices on 16 and 19 September in Figure (ref)(c) and (d), respectively. Unlike Figure (ref)(a) and (b), we note that the estimated spot idiosyncratic volatility matrices are highly sparse on the two days. This is consistent with our observation from Figure (ref)(d), confirming that the observed factors can effectively account for the time variation in the spot covariance structure of the returns. Meanwhile, we also note that the two idiosyncratic volatility matrices are clearly not diagonal and still carry some visible time variation. Lastly, it is worth mentioning that the estimated spot idiosyncratic volatility matrices do not exhibit significant correlations within the blocks along the diagonal lines, except for some very limited actions in the lower right corner of the two matrices (the lower right corner corresponds to the XLU sector according to our sorting).

Conclusion

We developed nonparametric estimation methods for large spot volatility matrices under the uniform sparsity assumption. We allowed for microstructure noise and observed common risk factors and employed kernel smoothing and generalised shrinkage. In each scenario we obtained the uniform convergence rates for the large estimated covariance matrices and these reflect the smoothness and sparsity assumptions we made. The simulation results show that the proposed estimation methods work well in finite samples for both the noise-free and noise-contaminated data. The empirical study demonstrated the effectiveness on S&P 500 stocks five minute data. Several issues can be further explored. For example, it is worthwhile to further study the spot precision matrix estimation which is briefly discussed in Appendix C.1 of the supplement and explore its application to optimal portfolio choice.

Acknowledgements

The authors would like to thank a Co-Editor and two reviewers for the constructive comments, which helped to improve the article. The first author's research was partly supported by the BA Talent Development Award (No. TDA21$\backslash$210027). The second author’s research was partly supported by the BA/Leverhulme Small Research Grant funded by the Leverhulme Trust (No. SRG1920/ 100603).