1 Introduction

The covariance matrix of asset returns is a central component, which is required in several contexts in financial economics, such as portfolio composition, pricing of financial instruments, or risk management (e.g., Andersen et al. (2003)). In the past few decades, estimation of quadratic variation (QV) from high-frequency data has been intensively studied in the field of econometrics. The standard estimator of QV, \(\Sigma \), say over a window [0, 1], is the realized variance (RV) (e.g., Andersen and Bollerslev (1998); Barndorff-Nielsen and Shephard (2002, 2004); Jacod (1994)). Given a decreasing sequence \((\Delta _{n})_{n \ge 1}\) with \(\Delta _{n} \rightarrow 0\) and equidistant observations \((Y_{i \Delta _{n}})_{i = 0}^{ \lfloor \Delta _{n}^{-1} \rfloor }\) of a d-dimensional semimartingale \((Y_{t})_{t \in [0,1]}\), the RV is defined as

$$\begin{aligned} \widehat{\Sigma }_{n} = \sum _{k=1}^{ \lfloor \Delta _{n}^{-1} \rfloor } ( \Delta _{k}^{n} Y) ( \Delta _{k}^{n} Y)^{ \top }, \end{aligned}$$
(1.1)

where \(\Delta _{k}^{n} Y = Y_{k \Delta _{n}} - Y_{(k-1) \Delta _{n}}\) is the increment, and \(^{ \top }\) denotes the transpose operator. The asymptotic properties of \(\widehat{\Sigma }_n\) are generally known, when the dimension d is fixed (e.g., Barndorff-Nielsen et al. (2006); Diop et al. (2013); Heiny and Podolskij (2020); Jacod (1994, 2008)). In particular, for any semimartingale \(\widehat{\Sigma }_n\) is by definition a consistent estimator of \(\Sigma \).

When modeling the dynamics of a large financial market over a relatively short time horizon, one is often faced with an ill-posed inference problem, where the number of assets d is comparable to (or even exceeding) the number of high-frequency data \(\lfloor \Delta _{n}^{-1} \rfloor \) available to compute \(\widehat{ \Sigma }_{n}\). This renders RV inherently inaccurate, and it also implies that d cannot be treated as fixed in the asymptotic analysis. To recover an estimate of \(\Sigma \) in such a high-dimensional setting, the covariance matrix has to display a more parsimonious structure.

Currently, the literature on high-dimensional high-frequency data, mostly based on a continuous semimartingale model for the security price process, is rather scarce. Wang and Zou (2010) investigate optimal convergence rates for estimation of QV under sparsity constraints on its entries. Zheng and Li (2011) apply random matrix theory to estimate the spectral eigenvalue distribution of QV in a one-factor diffusion setting. Related work based on random matrix theory and eigenvalue cleaning, or on approximating the high-dimensional system by a low-dimensional factor model, can be found in, e.g., Cai et al. (2020); Hautsch et al. (2012); Lunde et al. (2016). Aït-Sahalia and Xiu (2017) study a high-dimensional factor model, where the eigenvalues of QV are assumed to be dominated by common components (spiked eigenvalue setting), and estimate the number of common factors under sparsity constraints on the idiosyncratic components.

In the context of joint modeling of the cross-section of asset returns, if a factor model is adopted, one should expect a small number of eigenvalues of \(\Sigma \) to be dominant (see also Sect.  4). In other settings, such as when the market is complete and some entries of \(Y_{t}\) reflect prices on derivative contracts (spanned by the existing assets and therefore redundant), one should even expect the rank of the volatility to be strictly smaller than the dimension. In such an instance, the rank of the covariance matrix is smaller than d, meaning that the intrinsic dimension of the system can be represented in terms of a driving Brownian motion of lower dimension. In general, identifying a low rank can, for instance, help with the economic interpretation of the model, or it may lighten the computational load related to its estimation or simulation. We refer to several recent studies on identification and estimation of the intrinsic dimension of the driving Brownian motion and eigenvalue analysis of QV in Aït-Sahalia and Xiu (2017, 2019); Fissler and Podolskij (2017); Jacod et al. (2008); Jacod and Podolskij (2013, 2018), see also Kong (2017, 2020) in the high-dimensional setting. A related contribution of Pelger (2019) incorporates a finite-activity jump process in the model.

In this paper, we develop a regularized estimator of \(\Sigma \) called the penalized realized variance (PRV). We draw from the literature of shrinkage estimation in linear regression (e.g., Hoerl and Kennard (1970); Tibshirani (1996)). In particular, we propose the following LASSO-type estimator of \(\Sigma \) based on the RV in (1.1):

$$\begin{aligned} \widehat{\Sigma }_{n}^{ \lambda } = {{\,\mathrm{arg\,min}\,}}_{A \in \mathbb {R}^{d \times d}} \left( \Vert \widehat{\Sigma }_n - A \Vert _{2}^{2} + \lambda \Vert A \Vert _{1} \right) . \end{aligned}$$
(1.2)

Here, \(\lambda \ge 0\) is a tuning parameter that controls the degree of penalization, while \(\Vert \, \cdot \, \Vert _{1}\) and \(\Vert \, \cdot \, \Vert _{2}\) denote the nuclear and Frobenius matrix norm. It has been shown in various contexts that estimators based on nuclear norm regularization possess good statistical properties, such as having optimal rates of convergence and being of low rank. For instance, Koltchinskii et al. (2011) study such properties in the trace regression model and, in particular, in the matrix completion problem, where one attempts to fill out missing entries of a partially observed matrix. Negahban and Wainwright (2011) adopt a rather general observation model and, among other things, cover estimation in near low rank multivariate regression models and vector autoregressive processes. Further references in this direction include Argyriou et al. (2008); Bach (2008); Candès and Recht (2009); Recht et al. (2010). While previous papers on nuclear norm penalization are solely dedicated to discrete-time models, the current work is—to the best of our knowledge—the first to study this problem in a continuous-time Itô setting. This implies a number of subtle technical challenges, since the associated high-frequency data are dependent and potentially non-stationary.

We provide a complete non-asymptotic theory for the PRV estimator in (1.2) by deriving bounds for the Frobenius norm of its error. We show that it is minimax optimal up to a logarithmic factor. The derivation relies on advanced results from stochastic analysis combined with a Bernstein-type inequalities presented in Theorems 2.3 and 2.4, which constitute some of the major theoretical contributions of the paper. The latter builds upon recent literature on concentration inequalities for matrix-valued random variables, but it is more general as the variables are allowed to be both dependent and unbounded. We discuss the necessary conditions on d and \(\Delta _n\) to ensure consistency results. We further show that in a “local-to-deficient” rank setting, where some eigenvalues are possibly decreasing as the dimension of the system increases, the penalized estimator in (1.2) identifies the number of non-negligible eigenvalues of \(\Sigma \) with high probability.

The estimator depends on a tuning parameter \(\lambda \). We exploit the subsampling approach of Christensen et al. (2017) to choose a data-driven amount of shrinkage in the implementation. We also provide a related theoretical analysis for the local volatility, which is more likely to exhibit a lower rank than QV (the latter is only expected to have a near low rank in many settings).

In the empirical analysis, we look at the 30 members of the Dow Jones Industrial Average index to demonstrate that our results are consistent with a standard Fama–French three-factor structure, and one may expect a lower number of factors during crisis.

The paper proceeds as follows. In Sect. 2, we establish concentration inequalities for the estimation error \(\Vert \widehat{\Sigma }^\lambda _n - \Sigma \Vert _{2}\) and show minimax optimality up to a logarithmic factor. In Sect.  2.3, we present sufficient conditions to ensure that the rank of \(\widehat{\Sigma }^\lambda _n\) coincides with the number of non-negligible eigenvalues of \(\Sigma \) with high probability. In Sect. 2.4, we show how the tuning parameter \(\lambda \) can be selected with a data-driven approach. In Sect. 3, the associated non-asymptotic theory for estimation of the instantaneous variance is developed. In Sect. 4, we design a simulation study to demonstrate the ability of our estimator to identify a low-dimensional factor model through \({{\,\textrm{rank}\,}}( \widehat{ \Sigma }_{n}^{ \lambda })\). In Sect. 5, we implement the estimator on empirical high-frequency data from the 30 members of the Dow Jones Industrial Average index. The proofs are included in the appendix.

1.1 Notation

This paragraph introduces some notation used throughout the paper. For a vector or a matrix A the transpose of A is denoted by \(A^{ \top }\). We heavily employ the fact that any real matrix \(A \in \mathbb {R}^{m \times q}\) admits a singular value decomposition (SVD) of the form:

$$\begin{aligned} A = \sum _{k=1}^{m \wedge q} s_{k} u_{k} v_{k}^{ \top } \end{aligned}$$
(1.3)

where \(s_{1} \ge \cdots \ge s_{m \wedge q}\) are the singular values of A, and \(\{u_1,\dots , u_{m\wedge q}\}\) and \(\{v_1,\dots , v_{m\wedge q}\}\) are orthonormal vectors. If \(m = q\) and A is symmetric and positive semidefinite, its singular values coincide with its eigenvalues and the SVD is the orthogonal decomposition (meaning that one can take \(v_{k} = u_{k}\)). Sometimes we write \(s_{k}(A)\), \(v_{k}(A)\) and \(u_{k}(A)\) to explicitly indicate the dependence on the matrix A. The rank of A is denoted by \({{\,\textrm{rank}\,}}(A)\), whereas

$$\begin{aligned} {{\,\textrm{rank}\,}}(A; \varepsilon )= \sum _{k=1}^{m\wedge q} \mathbbm {1}_{[ \varepsilon , \infty )}(s_{k}) \end{aligned}$$
(1.4)

is the number of singular values of A exceeding a certain threshold \(\varepsilon \in [0, \infty )\). Note that \(\varepsilon \mapsto {{\,\textrm{rank}\,}}(A; \varepsilon )\) is non-increasing on \([0,\infty )\), and càdlàg and piecewise constant on \((0,\infty )\). Furthermore, \({{\,\textrm{rank}\,}}(A;0) = m\wedge q\) and \(\lim _{ \varepsilon \downarrow 0} {{\,\textrm{rank}\,}}(A;\varepsilon ) = {{\,\textrm{rank}\,}}(A)\). We denote by \(\Vert A \Vert _{p}\) the Schatten p-norm of \(A \in \mathbb {R}^{m \times q}\), i.e.,

$$\begin{aligned} \Vert A \Vert _{p} = \left( \sum _{k=1}^{m \wedge q} s_{k}^p\right) ^{1/p} \end{aligned}$$
(1.5)

for \(p \in (0,\infty )\). Moreover, \(\Vert A \Vert _{ \infty }\) denotes the maximal singular value of A, i.e., \(\Vert A\Vert _{\infty }= s_{1}\). In particular, in this paper we work with \(p=1\), \(p=2\), and \(p = \infty \) corresponding to the nuclear, Frobenius, and spectral norm. The Frobenius norm is also induced by the inner product \(\langle A,B \rangle = {{\,\textrm{tr}\,}}(A^{ \top } B)\) and we have the trace duality

$$\begin{aligned} \langle A, B \rangle \le \Vert A \Vert _{1} \Vert B \Vert _\infty . \end{aligned}$$
(1.6)

For a linear subspace \(S \subseteq \mathbb {R}^m\), \(S^{\perp }\) denotes the linear space orthogonal to S and \(P_{S}\) stands for the projection onto S.

Given two sequences \((a_{n})_{n \ge 1}\) and \((b_{n})_{n \ge 1}\) in \((0, \infty )\), we write \(a_{n} = o(b_{n})\) if \(\lim _{n \rightarrow \infty } a_{n}/b_{n} = 0\) and \(a_{n} = O (b_{n})\) if \(\limsup _{n \rightarrow \infty } a_{n}/b_{n}\) is bounded. Furthermore, if both \(a_{n} = O(b_{n})\) and \(b_{n} = O(a_{n})\), we write \(a_{n} \asymp b_{n}\). Finally, to avoid trivial cases we assume throughout the paper that the dimension d of the process \((Y_t)_{t\in [0,1]}\) is at least two.

2 Theoretical framework

We suppose a d-dimensional stochastic process \(Y_{t} = (Y_{1,t}, \dots , Y_{d,t})^{ \top }\) is defined on a filtered probability space \((\Omega , \mathcal {F}, ( \mathcal {F}_{t})_{t \in [0,1]}, \mathbb {P})\). In our context, \(Y_{t}\) is a time t random vector of log-prices of financial assets, which are traded in an arbitrage-free, frictionless market and therefore semimartingale (e.g., Delbaen and Schachermayer (1994)). In particular, we assume \((Y_t)_{t \in [0,1]}\) is a continuous Itô semimartingale, which can be represented by the stochastic differential equation:

$$\begin{aligned} \textrm{d}Y_{t} = \mu _{t} \textrm{d}t + \sigma _{t} \textrm{d} W_{t}, \qquad t \in [0,1], \end{aligned}$$
(2.1)

where \((W_{t})_{t \in [0,1]}\) is an adapted d-dimensional standard Brownian motion, \(( \mu _{t})_{t \in [0,1]}\) is a progressively measurable and locally bounded drift with values in \(\mathbb {R}^{d}\), and \((\sigma _{t})_{t \in [0,1]}\) is a predictable and locally bounded stochastic volatility with values in \(\mathbb {R}^{d \times d}\). In the above setting the QV (or integrated variance) is \(\Sigma = \int _{0}^{1} c_{t} \textrm{d}t\), where \(c_t \equiv \sigma _{t} \sigma _{t}^{ \top }\). We remark that the length \(T = 1\) of the estimation window [0, T] is fixed for expositional convenience and is completely without loss of generality, since a general T can always be reduced to the unit interval via a suitable time change.

Note that we exclude a jump component from (2.1). It should be possible to extend our analysis and allow for at least finite-activity jump processes. To do so, we can replace the RV with its truncated version, following the ideas of Mancini (2009). We leave the detailed analysis of more general jump processes for future work. Apart from these restrictions, our model is essentially nonparametric, as it allows for almost arbitrary forms of random drift, volatility, and leverage.

2.1 Upper bounds for the PRV

The goal of this section is to derive a sharp upper bound on the estimation error \(\Vert \widehat{ \Sigma }_{n}^{ \lambda } - \Sigma \Vert _{2}\) that apply with high probability, where the PRV estimator \(\widehat{ \Sigma }_{n}^{ \lambda }\) has been defined in (1.2). We first show that the PRV can alternatively be found by soft-thresholding of the eigenvalues in the orthogonal decomposition of \(\widehat{ \Sigma }_{n}\).

Proposition 2.1

Let \(\widehat{ \Sigma }_{n} = \sum _{k=1}^{d} s_{k} u_{k}u_{k}^{ \top }\) be the orthogonal decomposition of \(\widehat{ \Sigma }_{n}\). Then, the unique solution \(\widehat{ \Sigma }_{n}^{ \lambda }\) to (1.2) is given by

$$\begin{aligned} \widehat{ \Sigma }_{n}^{ \lambda } = \sum _{k=1}^{d} \max \left\{ s_{k} - \frac{ \lambda }{2},0 \right\} u_{k} u_{k}^{ \top }. \end{aligned}$$
(2.2)

Proposition 2.1 is standard (see Koltchinskii et al. (2011)), but we state it explicitly for completeness. The interpretation of the PRV representation in (2.2) is that all “significant” eigenvalues of \(\widehat{ \Sigma }_{n}\) are shrunk by \(\lambda /2\), while the smallest ones are set to 0. Hence, we only retain the principal eigenvalues in the orthogonal decomposition of \(\widehat{ \Sigma }_{n}\). The PRV can therefore account for a (near) low rank constraint. In the next proposition, we present a general non-probabilistic oracle inequality for the performance of \(\widehat{ \Sigma }_{n}^{ \lambda }\).

Proposition 2.2

Assume that \(2 \Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{ \infty } \le \lambda \). Then, it holds that

$$\begin{aligned} \Vert \widehat{ \Sigma }_{n}^{ \lambda } - \Sigma \Vert _{2}^{2} \le \inf _{A \in \mathbb {R}^{d \times d}} \left( \Vert A - \Sigma \Vert _{2}^{2} + \min \left\{ 2 \lambda \Vert A \Vert _{1}, 3 \lambda ^{2} {{\,\textrm{rank}\,}}(A) \right\} \right) . \end{aligned}$$
(2.3)

In particular, \(\Vert \widehat{ \Sigma }_{n}^{ \lambda } - \Sigma \Vert _{2}^{2} \le \min \left\{ 2 \lambda \Vert \Sigma \Vert _{1}, 3 \lambda ^{2} {{\,\textrm{rank}\,}}(\Sigma ) \right\} \).

The statement of Proposition 2.2 is also a general algebraic result that is standard for optimization problems under nuclear penalization (e.g., Koltchinskii et al. (2011)). It says that an oracle inequality for \(\widehat{ \Sigma }_{n}^{ \lambda }\) is available if we can control the stochastic error \(\Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{ \infty }\). The latter is absolute key to our investigation.

In order to assess this error, we need to impose the following assumption on the norms of the drift and volatility processes.

Assumption (A): \(\sup _{t\in [0,1]} \Vert \mu _{t} \Vert _{2}^{2} \le \nu _{ \mu }\), \(\sup _{t \in [0,1]} {{\,\textrm{tr}\,}}(c_{t}) \le \nu _{c,2}\), and \(\sup _{t \in [0,1]} \Vert c_{t} \Vert _{ \infty } \le \nu _{c, \infty }\) for some constants \(\nu _{ \mu }\), \(\nu _{c,2}\), and \(\nu _{c, \infty }\) with \(\nu _{c, \infty } \le \nu _{c,2}\).

Mathematically speaking, Assumption (A) appears a bit strong, since it imposes an almost sure bound on the drift and volatility. We can weaken the requirement on the drift to a suitable moment condition without affecting the rate of \(\Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{ \infty }\) in Theorem 2.4 below, but as usual the cost of this is more involved expressions. If the volatility does not meet the boundedness condition, which it does not for most stochastic volatility models, one can resort to the localization technique from Sect. 4.4.1 in Jacod and Protter (2012). In most financial applications, however, Assumption (A) is not too stringent if the drift and volatility do not vary strongly over short time intervals, such as a day.

The constants \(\nu _{ \mu }\), \(\nu _{c,2}\), and \(\nu _{c, \infty }\) may depend on d, but they should generally be as small as possible to get the best rate (see, e.g., (2.13)). For example, if the components of \(\mu \) and c are uniformly bounded, we readily deduce that

$$\begin{aligned} \nu _{ \mu } = O(d) \quad \text {and} \quad \nu _{c,2} = O(d). \end{aligned}$$
(2.4)

To study the concentration of \(\Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{ \infty }\), we present an exponential Bernstein-type inequality, which applies to matrix-valued martingale differences as long as the conditional moments are sufficiently well-behaved. While there are several related concentration inequalities (see, e.g., Minsker (2017); Tropp (2011, 2012, 2015) and references therein), existing results require that summands are either independent or bounded. Thus, to be applicable in our setting we needed to modify these.

Theorem 2.3

Suppose that \((X_{k})_{k=1}^{n}\) is a martingale difference sequence with respect to a filtration \(( \mathcal {G}_{k})_{k=0}^{n}\) and with values in the space of symmetric \(d \times d\) matrices. Assume that, for some predictable sequence \((C_{k})_{k=1}^{n}\) of random variables and constant \(R \in (0, \infty )\),

$$\begin{aligned} \big \Vert \mathbb {E} \left[ X_{k}^{p} \mid \mathcal {G}_{k-1} \right] \big \Vert _{ \infty } \le \frac{p!}{2} R^{p} C_{k}, \quad \text {{for} } k = 1, \dots , n \text { {and} } p = 2, 3, \dots \end{aligned}$$
(2.5)

Then, for all \(x,\nu \in (0, \infty )\), it holds that

$$\begin{aligned} \mathbb {P} \left( \Big \Vert \sum _{k=1}^{n} X_{k} \Big \Vert _{ \infty } > x,\, \sum _{k=1}^{n} C_{k} \le \nu \right) \le 2d \exp \left( - \frac{x^{2}}{2(xR+ \nu R^{2})} \right) . \end{aligned}$$
(2.6)

Assume that Theorem 2.3 applies and that \(C_{1}, \dots , C_{n}\) in (2.5) can also be chosen to be deterministic. Then, with \(\nu = \sum _{k=1}^{n} C_{k}\),

$$\begin{aligned} \mathbb {P} \left( \Big \Vert \sum _{k=1}^{n} X_{k} \Big \Vert _{ \infty } > x \right) \le 2d \exp \left( - \frac{x}{4 \nu R^{2}} \min \{x, \nu R \} \right) , \end{aligned}$$
(2.7)

so \(\Vert \sum _{k=1}^n X_k \Vert _\infty \) has a sub-exponential tail. The next theorem derives a concentration inequality for \(\Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{ \infty }\). We remark that, although we have Theorem 2.3 at our disposal, the derivation of Theorem 2.4 requires a number of non-standard inequalities in order to prove that R in (2.5) can be chosen sufficiently small. For further details, see the discussion in relation to its proof in the supplementary file.

Theorem 2.4

Suppose that Assumption (A) holds. Then, there exists an absolute constant \(\gamma \) such that

$$\begin{aligned} \mathbb {P} \left( \Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{ \infty } > x \right) \le 6d \exp \left( - \frac{x}{ \gamma \nu _{c,2} \nu _{c, \infty } \Delta _{n}} \min \{x , \nu _{c, \infty } \} \right) , \end{aligned}$$
(2.8)

for all \(\displaystyle x \ge \gamma \max \left\{ \frac{ \nu _{c,2} \nu _{\mu } \Delta _{n}}{ \nu _{c, \infty }}, \sqrt{ \nu _{c,2} \nu _{ \mu } \Delta _{n}}, \nu _{c, \infty } \Delta _{n} \right\} \).

We now combine Proposition 2.2 and Theorem 2.4 to deliver a result on the concentration of \(\Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{2}\), which is the main statement of this section.

Theorem 2.5

Suppose that Assumption (A) holds. For a given \(\tau \in (0,\infty )\) consider a regularization parameter \(\lambda \) such that

$$\begin{aligned} \lambda \ge \gamma \max \left\{ \sqrt{ \nu _{c,2} \nu _{c, \infty } \Delta _{n} ( \log (d) + \tau )}, \nu _{c,2} \Delta _{n} ( \log (d) + \tau ), \frac{ \nu _{c,2} \nu _{ \mu } \Delta _{n}}{ \nu _{c, \infty }} \right\} , \end{aligned}$$
(2.9)

where \(\gamma \) is a large absolute constant. Then, with probability at least \(1 - e^{- \tau }\),

$$\begin{aligned} \Vert \widehat{ \Sigma }_{n}^\lambda - \Sigma \Vert _{2}^{2} \le \Vert \Sigma - A \Vert _{2}^{2} + \min \{2 \lambda \Vert A \Vert _{1}, 3 \lambda ^{2} {{\,\textrm{rank}\,}}(A) \}, \end{aligned}$$
(2.10)

for all \(A \in \mathbb {R}^{d \times d}\).

If the drift term is non-dominant, in the sense that \(\nu _{ \mu } \le \nu _{c, \infty }\), it follows that the regularization parameter \(\lambda \) should meet

$$\begin{aligned} \lambda \ge \gamma \max \left\{ \sqrt{ \nu _{c,2} \nu _{c, \infty } \Delta _{n} \left( \log (d)+ \tau \right) }, \nu _{c,2} \Delta _{n}\left( \log (d)+ \tau \right) \right\} , \end{aligned}$$
(2.11)

for an absolute constant \(\gamma \). To get a concentration probability as large as possible without impairing the rate implied by (2.10), we should choose \(\tau \asymp \log (d)\). Moreover, the first term of the maximum in (2.11) is largest for \(1/ \Delta _{n} \ge ( \log (d)+ \tau ) \nu _{c,2}/ \nu _{c, \infty }\). In light of these observations, the following corollary is an immediate consequence of Theorem 2.5, so we exclude the proof.

Corollary 2.6

Suppose that Assumption (A) holds with \(\nu _{ \mu } \le \nu _{c, \infty } \le \nu _{c,2}\). Assume further that \(\Delta _{n}^{-1} \ge 2 \frac{ \nu _{c,2}}{ \nu _{c, \infty }} \log (d)\) and choose the regularization parameter

$$\begin{aligned} \lambda = \gamma \sqrt{ \nu _{c,2} \nu _{c, \infty } \Delta _{n} \log (d)}, \end{aligned}$$
(2.12)

where \(\gamma \) is a sufficiently large absolute constant. Then, it holds that

$$\begin{aligned} \Vert \widehat{ \Sigma }_{n}^{ \lambda } - \Sigma \Vert _{2}^{2} \le 3 \gamma ^{2} \nu _{c,2} \nu _{c, \infty } {{\,\textrm{rank}\,}}( \Sigma ) \Delta _{n} \log (d) \end{aligned}$$
(2.13)

with probability at least \(1-d^{-1}\).

It follows from (2.13) that the estimation error of \(\widehat{ \Sigma }_{n}^{ \lambda }\) is closely related to the size of \(\nu _{c,2}\) and \(\nu _{c, \infty }\), which are both determined from the properties of the volatility process. If c is uniformly bounded, \(\nu _{c,2} = O(d)\). As emphasized by Tropp (2015, Section 1.6.3), we can also often assume that \(\Vert c_{s} \Vert _{ \infty }\) and, hence, \(\nu _{c, \infty }\) can be chosen independently of d. When this is the case, the rate implied by (2.13) is \({{\,\textrm{rank}\,}}( \Sigma ) \Delta _{n}d \log (d)\), and since \({{\,\textrm{rank}\,}}( \Sigma ) \le d\), Corollary 2.6 implies consistency of \(\widehat{ \Sigma }_{n}^{ \lambda }\) when both \(\Delta _{n} \rightarrow 0\) and \(d \rightarrow \infty \) so long as \(d^{2} \log (d) = o( \Delta _{n}^{-1})\). If \({{\,\textrm{rank}\,}}( \Sigma )\) can be further bounded by a constant (that does not depend on d), the growth condition on d improves to \(d \log (d) = o( \Delta _{n}^{-1})\). In contrast, for the estimation error \(\Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{2}^{2} \rightarrow 0\), one cannot expect a better condition than \(d^{2} = o( \Delta _{n}^{-1})\), since this estimation error corresponds to a sum of \(d^{2}\) squared error terms, each of the order \(\Delta _{n}\).

2.2 Minimax optimality

We denote by \(\mathcal {C}_{r}\), for a given non-zero integer \(r \le d\), the subclass of \(d \times d\) symmetric positive semidefinite matrices \(\mathcal {S}_{+}^{d}\) whose effective rank is bounded by r:

$$\begin{aligned} \mathcal {C}_{r} = \left\{ A \in \mathcal {S}_{+}^{d} : r_e (A) \le r \right\} , \end{aligned}$$
(2.14)

where

$$\begin{aligned} r_{e} (A) \equiv \frac{{{\,\textrm{tr}\,}}(A)}{\Vert A \Vert _{ \infty }}. \end{aligned}$$
(2.15)

Compared to the rank, the effective rank is a more stable measure for the intrinsic dimension of a covariance matrix (see, e.g., (Vershynin 2010, Remark 5.53)). In the following we argue that \(\widehat{ \Sigma }_{n}^{ \lambda }\), with \(\lambda \) of the form (2.12), is a minimax optimal estimator of \(\Sigma \) in the parameters \(\nu _{c,2}\), \(\nu _{c, \infty }\), and \({{\,\textrm{rank}\,}}( \Sigma )\) over the parametric class of continuous Itô processes generated by (2.1) with no drift \(\mu _{s} \equiv 0\) and a constant volatility \(\sigma _{s} \equiv \sqrt{A}\) for \(A \in \mathcal {C}_{r}\). To this end, denote by \(\mathbb {P}_{A}\) a probability measure for which \((Y_{t})_{t \in [0,1]}\) is defined as in (2.1) with no drift and constant volatility \(c_{s} \equiv A\). In this setting, we can choose \(\nu _{c,2} = {{\,\textrm{tr}\,}}(A)\) and \(\nu _{c,\infty } = \Vert A \Vert _{\infty }\), which by Corollary 2.6 means that

$$\begin{aligned} \sup _{A \in \mathcal {C}_{r}} \mathbb {P}_{A} \left( \Vert \widehat{ \Sigma }_{n}^{ \lambda } - \Sigma \Vert _{2}^{2} > \overline{ \gamma } \Vert A \Vert _{ \infty }^{2} r_{e}(A) {{\,\textrm{rank}\,}}(A) \Delta _{n} \log (d) \right) \le \frac{1}{d}, \end{aligned}$$
(2.16)

for an absolute constant \(\overline{ \gamma }\) and \(\Delta _{n}^{-1} \ge 2r \log (d)\).

Now, since the log-price increments \(\Delta _{1}^{n} Y, \dots , \Delta _{ \lfloor \Delta _{n}^{-1} \rfloor }^{n} Y\) are i.i.d. Gaussian random vectors under \(\mathbb {P}_{A}\), we can exploit existing results from the literature to assess the performance of \(\widehat{ \Sigma }_{n}^{ \lambda }\). Hence, the following is effectively an immediate consequence of Theorem 2 in Lounici (2014). It shows that, up to a logarithmic factor, no estimator can do better than (2.16).

Theorem 2.7

Let \(\mathcal {C}_r\) be given as in (2.14) and suppose that \(\lfloor \Delta _n^{-1}\rfloor \ge r^2\). Then, there exist absolute constants \(\beta \in (0,1)\) and \(\underline{\gamma }\in (0,\infty )\) such that

$$\begin{aligned} \inf _{\widehat{\Sigma }}\sup _{A\in \mathcal {C}_r} \mathbb {P}_A \left( \Vert \widehat{\Sigma } - \Sigma \Vert _{2}^{2} > \underline{ \gamma } \frac{ \Vert A \Vert _{ \infty }^{2} r_e(A) {{\,\textrm{rank}\,}}(A)}{ \lfloor \Delta _{n}^{-1} \rfloor } \right) \ge \beta , \end{aligned}$$
(2.17)

where the infimum runs over all estimators \(\widehat{\Sigma }\) of \(\Sigma \) based on \(( \Delta _{1}^{n}Y, \dots , \Delta _{ \lfloor \Delta _{n}^{-1} \rfloor }^{n}Y)\).

2.3 Bound on the rank of PRV

In this section, we study the rank of \(\widehat{ \Sigma }_{n}^{ \lambda }\) relative to the number of non-negligible eigenvalues and, in particular, the true rank of \(\Sigma \). In line with Sect.  2.1, we begin by stating a general non-probabilistic inequality in Theorem 2.8. In the formulation of this result, we recall that \({{\,\textrm{rank}\,}}(A; \varepsilon )\) is the number of singular values of A exceeding \(\varepsilon \in [0, \infty )\).

Theorem 2.8

Suppose that \(2 \Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{ \infty } \le \bar{ \lambda }\) for some \(\bar{ \lambda } \in (0, \infty )\). Then

$$\begin{aligned} {{\,\textrm{rank}\,}}( \Sigma ; \lambda ) \le {{\,\textrm{rank}\,}}( \widehat{ \Sigma }_{n}^{ \lambda }) \le {{\,\textrm{rank}\,}}\left( \Sigma ; \frac{1}{2}( \lambda - \bar{ \lambda }) \right) , \end{aligned}$$
(2.18)

for any \(\lambda \in ( \bar{ \lambda }, \infty )\). In particular, if \(\lambda \in ( \bar{ \lambda }, s]\), where s is the smallest non-zero eigenvalue of \(\Sigma \), then \({{\,\textrm{rank}\,}}( \widehat{ \Sigma }_{n}^{ \lambda }) = {{\,\textrm{rank}\,}}( \Sigma )\).

With this result in hand we can rely on the exponential inequality for the quantity \(\Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{ \infty }\) established in Theorem 2.4 to show that, with high probability and in addition to converging to \(\Sigma \) at a fast rate, \(\widehat{ \Sigma }_{n}^{ \lambda }\) automatically has the rank of \(\Sigma \) when we neglect eigenvalues of sufficiently small order. In particular, if \(\Sigma \) has full rank, but many of its eigenvalues are close to zero, \(\widehat{ \Sigma }_{n}^{ \lambda }\) is of low rank and reflects the number of important components (or factors) in \(\Sigma \).

Corollary 2.9

Suppose that Assumption (A) holds with \(\nu _{ \mu } \le \nu _{c, \infty } \le \nu _{c,2}\) and fix \(\delta \in (0,1/2)\). Assume further that \(\Delta _{n}^{-1} \ge 2 \frac{ \nu _{c,2}}{ \nu _{c, \infty }} \log (d)\) and, for a sufficiently large constant \(\gamma \) depending only on \(\delta \), choose the regularization parameter \(\lambda \) as follows:

$$\begin{aligned} \lambda = \gamma \sqrt{ \nu _{c,2} \nu _{c, \infty } \Delta _{n} \log (d)}. \end{aligned}$$
(2.19)

Then, with probability at least \(1-d^{-1}\), it holds that

$$\begin{aligned} \Vert \widehat{ \Sigma }_{n}^{ \lambda } - \Sigma \Vert _{2}^{2} \le 3 \gamma ^{2} \nu _{c,2} \nu _{c, \infty } {{\,\textrm{rank}\,}}( \Sigma ) \Delta _{n} \log (d), \end{aligned}$$
(2.20)

and

$$\begin{aligned} {{\,\textrm{rank}\,}}( \Sigma ; \lambda ) \le {{\,\textrm{rank}\,}}( \widehat{ \Sigma }_{n}^{ \lambda }) \le {{\,\textrm{rank}\,}}( \Sigma ; \delta \lambda ). \end{aligned}$$
(2.21)

If, in addition, \(\lambda \le s\), where s is the smallest non-zero eigenvalue of \(\Sigma \), then both \({{\,\textrm{rank}\,}}( \widehat{ \Sigma }_{n}^{ \lambda }) = {{\,\textrm{rank}\,}}( \Sigma )\) and (2.20) hold with probability at least \(1-d^{-1}\).

In the setting of Corollary 2.9, it follows that with high probability an eigenvalue s of \(\Sigma \) affects \({{\,\textrm{rank}\,}}( \widehat{ \Sigma }_{n}^{ \lambda })\) for large n if \(\lambda = o(s)\), while it does not if \(s = o( \lambda )\). The first condition says that, relative to the level of shrinkage, an eigenvalue is significant (or non-negligible), whereas the second says the opposite. Hence, the notion of negligibility depends on \(\lambda \), which in turn depends on the model through the constants \(\nu _{c,2}\) and \(\nu _{c, \infty }\). However, we know that for \(\widehat{ \Sigma }_{n}^{ \lambda }\) to be a consistent estimator of \(\Sigma \), a necessary condition is that \(\lambda \rightarrow 0\) as \(n \rightarrow \infty \), which implies that for an eigenvalue of \(\Sigma \) to be negligible, it must tend to zero as n increases. In particular, if d is fixed, \({{\,\textrm{rank}\,}}( \widehat{ \Sigma }_{n}^{ \lambda }) = {{\,\textrm{rank}\,}}( \Sigma )\) with a probability tending to one as \(n \rightarrow \infty \).

The following example illustrates a stylized setting, where many eigenvalues of \(\Sigma \) are negligible.

Example 2.10

Let \(r\in \mathbb {N}\) be an absolute constant (i.e., independent of d and n). In line with the factor model studied in Aït-Sahalia and Xiu (2017), suppose that \(c_{t}\) is of the form

$$\begin{aligned} c_{t} = \beta e_{t} \beta ^{ \top } + g_{t}, \end{aligned}$$
(2.22)

where \(e_{t} \in \mathbb {R}^{r \times r}\) and \(g_{t} \in \mathbb {R}^{d \times d}\) are predictable, symmetric, and positive definite processes, which are uniformly bounded such that \(\Vert e_{t} \Vert _{ \infty } \le C_{e}\) and \(\Vert g_{t} \Vert _{ \infty } \le C_{g}\) for some constants \(C_{e}, C_{g} = O(1)\). The matrix \(\beta \in \mathbb {R}^{d \times r}\) of factor loadings is deterministic and constant in time; a common assumption in the literature, which is also supported by Reiss et al. (2015). The form (2.22) of \(c_{t}\) can be motivated by a standard multi-factor model, where the total risk associated with any of the d assets can be decomposed into a systematic and an idiosyncratic component. The systematic component \(\beta e_{t} \beta ^{ \top }\) is composed of loadings on the r underlying priced common factors in the economy, \(\beta \), and the risk of those factors, \(e_{t}\). The idiosyncratic component \(g_{t}\) is asset-specific and can therefore be reduced by diversification.

Suppose that, given an initial set of d individual securities, we start forming a new set of (orthogonalized) portfolios by taking linear combinations of the original assets, as prescribed by the APT model of Ross (1976), and a standard way to implement factor models in practice when constructing sorting portfolios. Then, these new assets are generally diversified enough to assume that \(C_{g} = O(d^{-1})\). To identify the number of driving factors, r, we assume that \(\Vert d^{- \gamma } \beta ^{ \top } \beta - I_{r} \Vert _{ \infty } = o(1)\) (as \(d \rightarrow \infty \)) and \(s_{r}(e_{t}) \ge \varepsilon \) for some absolute constants \(\gamma \in [0,1]\) and \(\varepsilon \in (0, \infty )\). This corresponds to Assumption 4 in Aït-Sahalia and Xiu (2017), but it is slightly more general as the assets are assumed to be diversified portfolios.

Now, given a suitable drift process \(( \mu _{t})_{t \in [0,1]}\), Corollary 2.9 applies with regularization parameter \(\lambda = O \left( d^{ \gamma } \sqrt{ \Delta _{n} \log (d)} \right) \). In particular, due to the fact \(0 \le s_{k} ( \Sigma ) - s_{k} ( \beta E \beta ^{ \top }) \le \Vert \Gamma \Vert _{ \infty }\) with \(E = \int _0^1 e_t \textrm{d} t\) and \(\Gamma = \int _0^1 g_t \textrm{d} t\), it follows from (2.21) that

$$\begin{aligned} {{\,\textrm{rank}\,}}( \beta E \beta ^{ \top }; \lambda ) \le {{\,\textrm{rank}\,}}( \widehat{ \Sigma }_{n}^{ \lambda }) \le {{\,\textrm{rank}\,}}(\beta E \beta ^\top ; \delta \lambda -\Vert \Gamma \Vert _\infty ) \end{aligned}$$
(2.23)

with probability at least \(1 - d^{-1}\). By assumption \(\vert s_{r}( \beta E \beta ^{ \top }) - d^{ \gamma } s_{r}(E) \vert = o(d^{ \gamma })\) and \(s_{r}(E) \ge \varepsilon \) (the former follows from arguments as in the proof of Theorem 1 in Aït-Sahalia and Xiu (2017)), and hence \(d^{ \gamma } = O(s_{r}( \beta E \beta ^{ \top }))\). Combining this with (2.23), we conclude that \({{\,\textrm{rank}\,}}( \widehat{ \Sigma }_{n}^{ \lambda }) = r\) with probability at least \(1 - d^{-1}\) if both \(\log (d) = o( \Delta _{n}^{-1})\) and \(\Delta _{n}^{-1}= o(d^{2+2 \gamma } \log (d))\), and n is large. Or to put it differently, \({{\,\textrm{rank}\,}}( \widehat{ \Sigma }_{n}^{ \lambda })\) is, with a probability tending to one, exactly the number of underlying factors in the model.

We remark that the restrictions imposed above are too weak to ensure that the Frobenius estimation error \(\Vert \widehat{ \Sigma }_{n}^{ \lambda } - \Sigma \Vert _{2}\) tends to zero, meaning that our estimator can be used to detect the underlying factor structure even in very high-dimensional settings. To ensure the estimation error is small, we need \(d^{1+2 \gamma } \log (d) = o( \Delta _{n}^{-1})\) rather than \(\log (d) = o( \Delta _{n}^{-1})\).

2.4 Selection of tuning parameter

In view of Theorem 2.7 and Corollary 2.9, it follows that \(\widehat{ \Sigma }_{n}^{ \lambda }\) can be of low rank and accurate, given optimal tuning of \(\lambda \). However, as evident from (2.19), \(\lambda \) depends on the latent spot variance process \((c_t)_{t \in [0,1]}\) through \(\nu _{c,2}\) and \(\nu _{c, \infty }\) as well as the unknown absolute constant \(\gamma \), whose value can be important in finite samples.

We remark that \(\widehat{ \Sigma }_{n} - \Sigma \) provides an estimate of the null matrix, but it is perturbed by randomness in the data. Hence, a good choice of shrinkage parameter exactly shrinks the eigenvalues of \(\widehat{ \Sigma }_{n} - \Sigma \) to zero (in view of problem (1.2) with \(\widehat{ \Sigma }_{n}\) replaced by \(\widehat{ \Sigma }_{n} - \Sigma \)). By Proposition 2.1, this means \(\lambda = 2 \Vert \widehat{ \Sigma }_{n} - \Sigma \Vert _{ \infty }\).

The above is unobservable. To circumvent this problem and facilitate the calculation of our estimator in Sect. 4 and 5, we adopt a data-driven shrinkage selection by exploiting the subsampling technique of Christensen et al. (2017), see also Politis et al. (1999) and Kalnina (2011). To explain the procedure in short, suppose for notational convenience that \(n = \Delta _{n}^{-1}\). We select an integer L—the number of subsamples—that divides n and assign log-returns successively to each subsample. Hence, the lth subsample consists of the increments \(\left( \Delta _{(k - 1)L + l}^{n} Y \right) _{k = 1, \ldots , n/L}\) for \(l = 1, \dots , L\). We denote the associated subsampled RV by

$$\begin{aligned} \widehat{\Sigma }_{n,l} = \frac{1}{n/L} \sum _{k=1}^{n/L} ( \sqrt{n} \Delta _{(k-1)L+l}^{n} Y) ( \sqrt{n} \Delta _{(k-1)L+l}^{n} Y)^{ \top }, \end{aligned}$$
(2.24)

Note that the random matrices in the sequence \((\widehat{ \Sigma }_{n,l})_{l = 1}^{L}\) are asymptotically conditionally independent, as \(n \rightarrow \infty \). Moreover, it follows from Christensen et al. (2017) that as \(n \rightarrow \infty \), \(L \rightarrow \infty \) and \(n/L \rightarrow \infty \), the half-vectorization of \(\sqrt{n/L}( \widehat{\Sigma }_{n,l} - \widehat{\Sigma }_{n})\) and \(\sqrt{n}(\widehat{\Sigma }_{n} - \Sigma )\) converge in law to a common mixed normal distribution, derived in Barndorff-Nielsen and Shephard (2004, Theorem 1).

Hence, in each subsample \(\lambda = 2 \Vert \widehat{\Sigma }_{n} - \Sigma \Vert _{ \infty }\) can be approximated by

$$\begin{aligned} \lambda _{l} = \frac{2}{ \sqrt{L}} \Vert \widehat{ \Sigma }_{n,l} - \widehat{ \Sigma }_{n} \Vert _{ \infty } \end{aligned}$$
(2.25)

As an estimator of \(\lambda \), we therefore take the sample average:

$$\begin{aligned} \lambda ^{ *} = \frac{1}{L} \sum _{l=1}^{L} \lambda _{l}. \end{aligned}$$
(2.26)

3 Estimation of the local variance

As noted in Sect. 1, the rank of the local variance \(c_{t}\) is often much smaller than the rank of \(\Sigma \). In fact, although \(\Sigma \) may be well-approximated by a matrix of low rank, we should expect that \({{\,\textrm{rank}\,}}(\Sigma ) = d\). On the other hand, there are many situations where it is reasonable to expect that \({{\,\textrm{rank}\,}}(c_{t})\) is small, in which case it provides valuable insight about the complexity of the underlying model. This motivates developing a theory for the spot version of our penalized estimator with particular attention on its ability to identify \({{\,\textrm{rank}\,}}(c_t)\).

To estimate \(c_{t}\), we follow Jacod and Protter (2012) and apply a localized realized variance, which is defined over a block of \(\lfloor h_{n} / \Delta _{n} \rfloor \ge 2\) log-price increments with \(h_{n} \in (0,1)\). The RV computed over the time window \([t, t+h_{n}]\), for \(t \in (0,1-h_{n}]\), is then defined as follows:

$$\begin{aligned} \widehat{ \Sigma }_{n}(t;t+h_n) = \sum _{k = \lfloor t / \Delta _{n} \rfloor +1}^{ \lfloor (t+h_{n}) / \Delta _{n} \rfloor } ( \Delta _{k}^{n} Y)( \Delta _{k}^{n} Y)^{ \top }. \end{aligned}$$
(3.1)

The corresponding penalized version \(\widehat{ \Sigma }_{n}^{ \lambda } (t;t+h_{n})\) is computed as

$$\begin{aligned} \widehat{ \Sigma }_{n}^{ \lambda } (t;t+h_{n}) = {{\,\mathrm{arg\,min}\,}}_{A \in \mathbb {R}^{d \times d}} \left( \Vert \widehat{ \Sigma }_n(t;t+h_{n}) - A \Vert _{2}^{2} + \lambda \Vert A \Vert _{1} \right) , \end{aligned}$$
(3.2)

or by using Proposition 2.1 with \(\widehat{ \Sigma }_{n}\) replaced by \(\widehat{ \Sigma }_{n}(t;t+h_{n})\). Then, the penalized estimator \(\hat{c}_{n}^{ \lambda }(t)\) of \(c_{t}\) is given by

$$\begin{aligned} \hat{c}_{n}^{ \lambda }(t) = \frac{ \widehat{ \Sigma }_{n}^{ \lambda } (t;t+h_{n})}{h_{n}}. \end{aligned}$$
(3.3)

Also, we write

$$\begin{aligned} \Sigma (t_{1};t_{2}) = \int _{t_{1}}^{t_{2}} c_{s} \textrm{d}s, \end{aligned}$$
(3.4)

for the QV over \((t_{1},t_{2}]\) for arbitrary time points \(t_{1},t_{2} \in [0,1]\) with \(t_{1} < t_{2}\).

Recall that for a convex and increasing function \(\psi :[0, \infty ) \rightarrow [0, \infty )\) with \(\psi (0) = 0\), the Orlicz norm of a random variable Z with values in the Hilbert space \(( \mathbb {R}^{d \times d}, \Vert \, \cdot \, \Vert _{2})\) is defined as:

$$\begin{aligned} \Vert Z \Vert _{ \psi } \equiv \inf \left\{ c > 0 : \mathbb {E}[ \psi ( \Vert Z \Vert _{2}/c)] \le 1 \right\} . \end{aligned}$$
(3.5)

This setting includes, in particular, the \(L^{p}( \mathbb {P})\) norms (with \(\psi (s) = s^{p}\)), but also the \(\psi _{1}\) and \(\psi _{2}\) norms for sub-exponential (\(\psi (s) = e^{s}-1\)) and sub-Gaussian (\(\psi (s) = e^{s^{2}}-1\)) random variables. For convenience, we further impose the mild restriction that the range of \(\psi \) is \([0, \infty )\), so it admits an inverse \(\psi ^{-1}\) on \([0, \infty )\).

Theorem 3.1

Suppose that Assumption (A) holds with \(\nu _{ \mu } \le \nu _{c, \infty } \le \nu _{c,2}\). Furthermore, let \(t \in \big [ \frac{h_n}{2},1- \frac{3h_n}{2} \big ]\) and assume \(\sup _{0 \le u<s \le 1} \Vert c_{s} - c_{u} \Vert _{ \psi }/ \sqrt{s-u} \le \nu _{c, \psi }\) and \(h_{n} \ge 2 \Delta _{n} \frac{ \nu _{c,2}}{ \nu _{c, \infty }} \log (d)\). Take the regularization parameter \(\lambda \) as

$$\begin{aligned} \lambda = \gamma \sqrt{h_{n} \Delta _{n} \nu _{c,2} \nu _{c, \infty } \log (d)} \end{aligned}$$
(3.6)

for a sufficiently large absolute constant \(\gamma \). Then,

$$\begin{aligned} \Vert \hat{c}_{n}^{ \lambda }(t) - c_{t} \Vert _{2}^{2} \le \kappa \gamma ^{2} \left( \frac{ \nu _{c,2} \nu _{c, \infty } \Delta _{n} {{\,\textrm{rank}\,}}\left( \Sigma (t- \frac{h_n}{2};t+ \frac{3h_n}{2}) \right) }{h_{n}} + h_{n} \nu _{c, \psi }^{2} \right) \log (d), \end{aligned}$$
(3.7)

with probability at least \(1-d^{-1} - \psi ( \sqrt{ \log (d)})^{-1}\) for some absolute constant \(\kappa \).

The bound on the estimation error of \(\hat{c}_{n}^{ \lambda }(t)\) builds on Corollary 2.6, but Theorem 3.1 further enforces a smoothness condition on \((c_{t})_{t \in [0,1]}\). It entails a trade-off between the smoothness of spot variance and the concentration probability associated with (3.7). For example, if \((c_{t})_{t \in [0,1]}\) is 1/2-Hölder continuous in \(L^{2}( \mathbb {P})\), which corresponds to setting \(\psi (s) = s^{2}\), then we end up with a concentration probability converging to one at rate \(\log (d)^{-1}\), which is slower than for the PRV. On the other hand, if \((c_{t})_{t \in [0,1]}\) is 1/2-Hölder continuous in the sub-Gaussian norm \(\psi (s) = e^{s^{2}}-1\), the concentration probability converges to one at the rate \(d^{-1}\), which is equivalent to the PRV.

Note that with d fixed, (3.7) reveals how to select the length of the estimation window \(h_{n}\) optimally. The upper bound depends on \(\Delta _{n}/h_{n}\) and \(h_{n}\), and for these terms to converge to zero equally fast, we should take \(h_{n} \asymp \sqrt{ \Delta _{n}}\). This is consistent with the literature on spot variance estimation; see, e.g., Jacod and Protter (2012).

The next result, which relies on Corollary 2.9, shows that the rank and performance of \(\hat{c}_{n}^{ \lambda }(t)\) can be controlled simultaneously.

Theorem 3.2

Suppose that Assumption (A) holds with \(\nu _{ \mu } \le \nu _{c, \infty } \le \nu _{c,2}\) and fix \(\delta \in (0,1/4)\). Furthermore, let \(t \in [ \frac{h_n}{2},1 - \frac{3h_n}{2}]\) and assume that \(\sup _{0 \le u < s \le 1} \Vert c_{s} - c_{u} \Vert _{ \psi } / \sqrt{s-u} \le \nu _{c, \psi }\). Suppose further that \(h_{n} \ge 2 \Delta _{n} \frac{ \nu _{c,2}}{ \nu _{c, \infty }} \log (d)\) and consider a regularization parameter \(\lambda \) such that

$$\begin{aligned} \lambda = \gamma \sqrt{h_{n} \Delta _{n} \nu _{c,2} \nu _{c, \infty } \log (d)}, \end{aligned}$$
(3.8)

for a sufficiently large constant \(\gamma \) depending only on \(\delta \). Then, with probability at least \(1-d^{-1}- \psi ( \sqrt{ \log (d)})^{-1}\), it holds that

$$\begin{aligned} \Vert \hat{c}^{ \lambda }_{n}(t) - c_{t} \Vert _{2}^{2} \le \kappa \gamma ^{2} \left( \frac{ \nu _{c,2} \nu _{c, \infty } \Delta _{n} {{\,\textrm{rank}\,}}\left( \Sigma (t- \frac{h_{n}}{2};t+ \frac{3h_{n}}{2}) \right) }{h_{n}} + h_{n} \nu _{c, \psi }^{2} \right) \log (d), \end{aligned}$$
(3.9)

and

$$\begin{aligned} {{\,\textrm{rank}\,}}(c_{t}; \overline{ \varepsilon }) \le {{\,\textrm{rank}\,}}( \hat{c}_{n}^{ \lambda }(t)) \le {{\,\textrm{rank}\,}}(c_{t}; \underline{ \varepsilon }), \end{aligned}$$
(3.10)

where

$$\begin{aligned} \begin{aligned} \underline{ \varepsilon }&\equiv \delta \gamma \max \biggl \{ \sqrt{ \frac{ \nu _{c,2} \nu _{c, \infty } \Delta _{n}}{h_{n}}} - \nu _{c, \psi } \sqrt{h_{n}},0 \biggr \} \sqrt{ \log (d)}, \\ \overline{ \varepsilon }&\equiv \gamma \left( \sqrt{ \frac{ \nu _{c,2} \nu _{c, \infty } \Delta _{n}}{h_{n}}} + \nu _{c, \psi } \sqrt{h_{n}} \right) \sqrt{ \log (d)}, \end{aligned} \end{aligned}$$
(3.11)

and \(\kappa \) is an absolute constant. If, in addition, \(\underline{ \varepsilon } > 0\) and \(\overline{ \varepsilon } \le s_{ {{\,\textrm{rank}\,}}(c_{t})} (c_{t})\), then both (3.9) and \({{\,\textrm{rank}\,}}(\hat{c}_{n}^{ \lambda }(t)) = {{\,\textrm{rank}\,}}(c_{t})\) hold with probability at least \(1-d^{-1}- \psi ( \sqrt{ \log (d)})^{-1}\).

Consider the setting of Theorem 3.2. For the upper bound in (3.10) to be useful, we must have that \(\underline{ \varepsilon } > 0\) or, equivalently,

$$\begin{aligned} \nu _{c, \psi }^{2} h_{n}^{2} < \nu _{c,2} \nu _{c, \infty } \Delta _{n}. \end{aligned}$$
(3.12)

Suppose that \(\nu _{c,2}\), \(\nu _{c,\infty }\), and \(\nu _{c,\psi }\) do not depend on n. The inequality (3.12) can then be achieved in large samples by choosing \(h_{n} = o( \sqrt{ \Delta _{n}})\). However, as pointed out above one should take \(h_{n} \asymp \sqrt{ \Delta _{n}}\) in order to achieve an optimal rate for the estimation error \(\Vert \hat{c}_{n}^{ \lambda }(t) - c_{t} \Vert _{2}^{2}\). In that case, (3.12) translates directly into an upper bound on \(\nu _{c, \psi }\), which concerns the smoothness of \((c_t)_{t \in [0,1]}\).

When can the term \({{\,\textrm{rank}\,}}( \Sigma (t- \frac{h_{n}}{2}; t+ \frac{3h_{n}}{2}))\) in the upper bound of (3.9) be replaced by \({{\,\textrm{rank}\,}}(c_t)\)? This question appears difficult to answer in general, but it is not too difficult to show that the replacement is valid if the volatility process is locally constant.

4 Simulation study

In this section, we do a Monte Carlo analysis to inspect the finite sample properties of the PRV, \(\widehat{\Sigma }_n^\lambda \). The aim here is to demonstrate the ability of our estimator to detect the presence of near deficient rank in a large-dimensional setting. In view of Theorem 3.2, this can be achieved by a proper selection of \(\lambda \). We simulate a \(d = 30\)-dimensional log-price process \(Y_{t}\). The size of the vector corresponds to the number of assets in our empirical investigation. We assume \(Y_{t}\) follows a slightly modified version of the \(r = 3\)-factor model proposed in Aït-Sahalia and Xiu (2019):

$$\begin{aligned} \textrm{d} Y_{t} = \beta \textrm{d} F_{t} + \textrm{d} Z_{t}, \quad t \in [0,1], \end{aligned}$$
(4.1)

where \(F_{t} \in \mathbb {R}^{r}\) is a vector of systematic risk factors, which has the dynamic:

$$\begin{aligned} \textrm{d} F_{t} = \alpha \textrm{d} t + \sigma _{t} L \textrm{d} W_{t}, \end{aligned}$$
(4.2)

and \(W_{t} \in \mathbb {R}^{r}\) is a standard Brownian motion. The drift term \(\alpha \in \mathbb {R}^{r}\) is constant \(\alpha = (0.05,0.03,0.02)^{ \top }\). The random volatility matrix \(\sigma _{t} = {{\,\textrm{diag}\,}}( \sigma _{1,t}, \sigma _{2,t}, \sigma _{3,t}) \in \mathbb {R}^{r \times r}\), where \({{\,\textrm{diag}\,}}( \, \cdot \,)\) is a diagonal matrix with coordinates \(\,\cdot \,\), is simulated as a Heston (1993) square-root process:

$$\begin{aligned} \textrm{d} \sigma _{j,t}^{2} = \kappa _{j} ( \theta _{j} - \sigma _{j,t}^{2}) \textrm{d}t + \eta _{j} \sigma _{j,t} \left( \rho _{j} \textrm{d} W_{j,t} + \sqrt{1 - \rho _{j}^{2}} \textrm{d} \widetilde{W}_{j,t} \right) , \end{aligned}$$
(4.3)

for \(j = 1, \dots , r\), where \(\widetilde{W}_{t} \in \mathbb {R}^{r}\) is an r-dimensional standard Brownian motion independent of \(W_{t}\). As in Aït-Sahalia and Xiu (2019), the parameters are \(\kappa = (3, 4, 5)^{ \top }\), \(\theta = (0.05, 0.04, 0.03)^{ \top }\), \(\eta = (0.3, 0.4, 0.3)^{ \top }\), and \(\rho = (-0.60, -0.40, -0.25)^{ \top }\). The factor correlation is captured by the Cholesky component L:

$$\begin{aligned} \rho = L L^{ \top } = \begin{bmatrix} 1.00 &{} 0.05 &{} 0.10 \\ 0.05 &{} 1.00 &{} 0.15 \\ 0.10 &{} 0.15 &{} 1.00 \end{bmatrix}. \end{aligned}$$
(4.4)

The idiosyncratic component \(Z_{t} \in \mathbb {R}^{d}\) is given by

$$\begin{aligned} \textrm{d} Z_{t} = \sqrt{g_t} \textrm{d} B_{t}, \end{aligned}$$
(4.5)

where \(g_{t} = {{\,\textrm{diag}\,}}( \gamma _{1,t}^2, \dots , \gamma _{d,t}^2)\) with

$$\begin{aligned} \textrm{d} \gamma _{j,t}^{2} = \kappa _{Z} ( \theta _{Z} - \gamma _{j,t}^{2}) \textrm{d}t + \eta _{Z} \gamma _{j,t} \textrm{d} \widetilde{B}_{j,t}, \end{aligned}$$
(4.6)

for \(j = 1, \dots , d\), and \(B_{t} \in \mathbb {R}^{d}\) and \(\widetilde{B}_{t} \in \mathbb {R}^{d}\) are independent d-dimensional standard Brownian motions. Thus, \(g_{t}\) is the instantaneous variance of the unsystematic component. We set \(\kappa _{Z} = 4\), \(\theta _{Z} = 0.25\), and \(\eta _{Z} = 0.06\). Hence, the idiosyncratic error varies independent in the cross-section of assets, but the parameters controlling the dynamics are common.

Fig. 1
figure 1

Relative importance of eigenvalues in simulated model. Note. In Panel A, we plot for daily RV the relative size of the three largest eigenvalues, and also the average of the remaining 27 eigenvalues. In Panel B, we plot a histogram of the sample average of the relative size of each of the eigenvalues across simulations

Full size image

The above setup implies

$$\begin{aligned} c_{t} = \beta \sigma _{t} \rho \sigma _{t} \beta ^{ \top } + g_{t}, \end{aligned}$$
(4.7)

so the spot covariance matrix has full rank at every time point t. However, it has only \(r = 3\) large eigenvalues associated with the systematic factors, while the remaining \(d-r\) associated with the idiosyncratic variance are relatively small. Note that in contrast to Aït-Sahalia and Xiu (2019), we allow for time-varying idiosyncratic variance.

We set \(\Delta _{n} = 1/n\) with \(n = 78\). This corresponds to 5-minute sampling frequency within a 6.5 hour window. Hence, d is relatively large compared to n. We construct 10,000 independent replications of this model. At the beginning of each simulation, we draw the associated factor loadings \(\beta \in \mathbb {R}^{d \times r}\) at random such that the first column (interpreted as the loading on the market factor) are from a uniform distribution on (0.25, 1.75), i.e. \(\beta _{i1} \sim U(0.25,1.75)\), \(i = 1, \dots , d\). The range of the support is consistent with the realized beta values reported in Table 1 in Sect.  5. The remaining columns are generated as \(\beta _{ij} \sim N(0, 0.5^{2})\), \(i = 1, \dots , d\) and \(j = 2\) and 3.

In Panel A of Fig. 1, we plot the relative size of the three largest eigenvalues, extracted from the corresponding RV, together with the average of the remaining 27 eigenvalues. In Panel B, we include a histogram of the relative size of the eigenvalues across simulations. We observe that it is generally challenging to distinguish important factors from idiosyncratic variation, since the relative size of the eigenvalues decays smoothly with the exception of the largest (and perhaps the second largest) eigenvalue.

In each simulation, we employ the subsampling procedure described in Sect. 2.4 with \(L = 6\) to choose \(\lambda ^{*}\).Footnote 1

Fig. 2
figure 2

Properties of the PRV. Note. In Panel A, we report a kernel density estimate of the shrinkage parameter \(\lambda \), expressed in percent of the maximal eigenvalue of the RV matrix. In Panel B, we show the relative frequency histogram of the associated rank of the PRV

Full size image

In Panel A of Fig. 2, we report the distribution of the \(\lambda ^{ *}\) parameter across simulations, which is here expressed in percent of the maximal eigenvalue of RV. Overall, the penalization is relatively stable with a tendency to shrink around one-fourth to one-third of the largest eigenvalue most of the times.

In Panel B, we plot a histogram of the relative frequency of the rank of the PRV. In general, the PRV does a good job at identifying the number of driving factors, given the challenging environment. Although there are variations in the rank across simulations, caused by the various sources of randomness incorporated in our model, the picture is close to the truth. This means we tend to identify about one-to-three driving factors. The average rank estimate is 1.79, so as expected there is a slight tendency to overshrink leading to a small downward bias. Based on these findings, we can confidently apply the PRV in the empirical application.

5 Empirical application

In this section, we apply the PRV to an empirical dataset. The sample period is from January 3, 2007 through May 29, 2020 (3,375 trading days in total) and includes both the financial crisis around 2007 – 2009 and partly the recent events related to Covid-19.

Table 1 Descriptive statistics of high-frequency data
Full size table

We look at high-frequency data from the 30 members of the Dow Jones Industrial Average (DJIA) index as of April 6, 2020.Footnote 2 The ticker codes of the various firms are listed in Table 1, along with selected descriptive statistics. As readily seen from Table 1, most of these companies are very liquid. In order to compute a daily RV matrix estimate we collect the individual 5-minute log-return series spanning the 6.5 hours of trading on U.S. stock exchanges from 9:30am to 4:00pm EST.Footnote 3 Hence, our analysis is based on a relatively large cross-section (i.e., \(d = 30\)) compared to the sampling frequency (i.e., \(n = 78\)). The realized beta in Table 1 is calculated with SPY acting as market portfolio. The latter is an exchange-traded fund that tracks the S &P500 and its evolution is representative of the overall performance of the U.S. stock market. The dispersion of realized beta is broadly consistent with the simulated values generated in Sect. 4.

Fig. 3
figure 3

Proportion of variance explained by each eigenvalue. Note. In Panel A, we plot for each day in the sample the relative size of the three largest eigenvalues of RV, as well as the average of the remaining 27. In Panel B, we plot a histogram of the time series average of the relative size of each eigenvalues

Full size image

In Fig. 3, we depict the factor structure in the time series of RV. In Panel A, we compute on a daily basis the proportion of the total variance (i.e., trace) explained by each eigenvalue, whereas Panel B reports the sample average of this statistic. As consistent with Aït-Sahalia and Xiu (2019), we observe a pronounced dynamic evolution in the contribution of each eigenvalue to RV with notably changes corresponding to times of severe distress in financial markets. The first factor captures the vast majority of aggregate return variation (about 35% on average). It is followed by a few smaller—but still important—factors accountable for an incremental 25% – 30% of the total variance, whereas the last 25 or so eigenvalues are relatively small.

Fig. 4
figure 4

Rank of PRV. Note. In Panel A, we report the relative frequency of the rank of the PRV. In Panel B, we compare the estimated rank to the effective rank, where both are smoothed over a three-month moving average window

Full size image

Next, we turn our attention to the PRV estimator. To select the shrinkage parameter \(\lambda \), we follow the subsampling implementation from the simulation section. The relative frequency histogram of the rank of the PRV is reported in Panel A of Fig.  4, whereas Panel B reports a three-month (90-day) moving average of the rank. The vast majority (around 95%) of the daily rank estimates are between one and three, which is consistent with a standard Fama–French three-factor interpretation of the data. There are relatively few rank estimates at four, and it never exceeds five (with about a dozen of the latter). In Panel B, we observe the rank varies over time in an intuitive fashion, often dropping close to a single-factor representation during times of crisis, where correlations tend to rise. As a comparison, we compute the effective rank of the RV (cf. (2.15)), which does not depend on a shrinkage parameter. There is a high association between the series, which is consistent with a soft-thresholding that eliminates smaller eigenvalues of the RV.

6 Conclusion

In this paper, we develop a novel and powerful penalized realized variance estimator, which is applicable to estimate the quadratic variation of high-dimensional continuous-time semimartingales under a low rank constraint. Our estimator relies on regularization and adapts the principle ideas of the LASSO from regression analysis to the field of high-frequency volatility estimation in a high-dimensional setting. We derive a non-asymptotic analysis of our estimator, including bounds on its estimation error and rank. The estimator is found to be minimax optimal up to a logaritmic factor. We design a completely data-driven procedure for selecting the shrinkage parameter based on a subsampling approach. In a simulation study, the estimator is found to possess good properties. In our empirical application, we confirm a low rank environment that is consistent with a three-factor structure in the cross-section of log-returns from the large-cap segment of the U.S. stock market.