Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
68,130 characters · 11 sections · 16 citation commands
High-dimensional estimation of quadratic variation based on penalized realized variance
\thispagestyle{empty}
\setcounter{page}{1}
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-bollerslev-diebold-labys:03a). 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-bollerslev:98a, barndorff-nielsen-shephard:02a, barndorff-nielsen-shephard:04a,jacod:94a). Given a decreasing sequence $(\Delta_{n})_{n \geq 1}$ with $\Delta_{n} \to 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
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-graversen-jacod-podolskij-shephard:06a, diop-jacod-todorov:13a, heiny-podolskij:20a, jacod:94a, jacod:08a). 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. \citet*{wang-zou:10a} investigate optimal convergence rates for estimation of QV under sparsity constraints on its entries. \citet*{zheng-li:11a} 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-hu-li-zheng:20a, hautsch-kyj-oomen:12a, lunde-shephard-sheppard:16a. \citet*{ait-sahalia-xiu:17a} 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 Section (ref)). 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 ait-sahalia-xiu:17a, ait-sahalia-xiu:19b, fissler-podolskij:17a,jacod-lejay-talay:08a,jacod-podolskij:13a, jacod-podolskij:18a, see also \citet*{kong:17a,kong:20a} in the high-dimensional setting. A related contribution of \citet*{pelger:19a} 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-kennard:70a, tibshirani:96a). In particular, we propose the following LASSO-type estimator of $\Sigma$ based on the RV in (ref):
Here, $\lambda \geq 0$ is a tuning parameter that controls the degree of penalization, while $\lVert \: \cdot \: \rVert_{1}$ and $\lVert \: \cdot \: \rVert_{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, \citet*{koltchinskii-lounici-tsybakov:11a} 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. \citet*{negahban-wainwright:11a} 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-evgeniou-pontil:08a, bach:08a, candes-recht:09a, recht-fazel-parrilo:10a. 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\^{o} 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 (ref) 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 (ref) and (ref), 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 (ref) 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 \citet*{christensen-podolskij-thamrongrat-veliyev:17a} 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 Section (ref), we establish concentration inequalities for the estimation error $\lVert \widehat{\Sigma}^\lambda_n - \Sigma \rVert_{2}$ and show minimax optimality up to a logarithmic factor. In Section (ref), 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 Section (ref), we show how the tuning parameter $\lambda$ can be selected with a data-driven approach. In Section (ref), the associated non-asymptotic theory for estimation of the instantaneous variance is developed. In Section (ref), we design a simulation study to demonstrate the ability of our estimator to identify a low-dimensional factor model through $\operatorname{rank} ( \widehat{ \Sigma}_{n}^{ \lambda})$. In Section (ref), 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.
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:
where $s_{1} \geq \cdots \geq 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 $\operatorname{rank}(A)$, whereas
is the number of singular values of $A$ exceeding a certain threshold $\varepsilon \in [0, \infty)$. Note that $\varepsilon \mapsto \operatorname{rank} (A; \varepsilon)$ is non-increasing on $[0,\infty)$, and c\`{a}dl\`{a}g and piecewise constant on $(0,\infty)$. Furthermore, $\operatorname{rank} (A;0) = m\wedge q$ and $\lim_{ \varepsilon \downarrow 0} \operatorname{rank} (A;\varepsilon) = \operatorname{rank}(A)$. We denote by $\lVert A \rVert_{p}$ the Schatten $p$-norm of $A \in \mathbb{R}^{m \times q}$, i.e.,
for $p \in (0,\infty)$. Moreover, $\lVert A \rVert_{ \infty}$ denotes the maximal singular value of $A$, i.e., $\lVert A\rVert_{\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 = \operatorname{tr} (A^{ \top} B)$ and we have the trace duality
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 \geq 1}$ and $(b_{n})_{n \geq 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.
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-schachermayer:94a). In particular, we assume $(Y_t)_{t \in [0,1]}$ is a continuous It\^{o} semimartingale, which can be represented by the stochastic differential equation:
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} \mathrm{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 (ref). 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 \citet*{mancini:09a}. 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.
The goal of this section is to derive a sharp upper bound on the estimation error $\lVert \widehat{ \Sigma}_{n}^{ \lambda} - \Sigma \rVert_{2}$ that apply with high probability, where the PRV estimator $\widehat{ \Sigma}_{n}^{ \lambda}$ has been defined in (ref). 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 (ref) is standard (see koltchinskii-lounici-tsybakov:11a), but we state it explicitly for completeness. The interpretation of the PRV representation in (ref) 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}$.
The statement of Proposition (ref) is also a general algebraic result that is standard for optimization problems under nuclear penalization (e.g., koltchinskii-lounici-tsybakov:11a). It says that an oracle inequality for $\widehat{ \Sigma}_{n}^{ \lambda}$ is available if we can control the stochastic error $\lVert \widehat{ \Sigma}_{n} - \Sigma \rVert_{ \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.
Mathematically speaking, Assumption (ref) 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 $\lVert \widehat{ \Sigma}_{n} - \Sigma \rVert_{ \infty}$ in Theorem (ref) 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 Section 4.4.1 in \citet*{jacod-protter:12a}. In most financial applications, however, Assumption (ref) 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., (ref)). For example, if the components of $\mu$ and $c$ are uniformly bounded, we readily deduce that
To study the concentration of $\lVert \widehat{ \Sigma}_{n} - \Sigma \rVert_{ \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:17a, tropp:11a, tropp:12a, tropp:15a 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.
Assume that Theorem (ref) applies and that $C_{1}, \dots, C_{n}$ in (ref) can also be chosen to be deterministic. Then, with $\nu = \sum_{k=1}^{n} C_{k}$,
so $\lVert\sum_{k=1}^n X_k \rVert_\infty$ has a sub-exponential tail. The next theorem derives a concentration inequality for $\lVert \widehat{ \Sigma}_{n} - \Sigma \rVert_{ \infty}$. We remark that, although we have Theorem (ref) at our disposal, the derivation of Theorem (ref) requires a number of non-standard inequalities in order to prove that $R$ in (ref) can be chosen sufficiently small. For further details, see the discussion in relation to its proof in the supplementary file.
We now combine Proposition (ref) and Theorem (ref) to deliver a result on the concentration of $\lVert \widehat{ \Sigma}_{n} - \Sigma \rVert_{2}$, which is the main statement of this section.
If the drift term is non-dominant, in the sense that $\nu_{ \mu} \leq \nu_{c, \infty}$, it follows that the regularization parameter $\lambda$ should meet
for an absolute constant $\gamma$. To get a concentration probability as large as possible without impairing the rate implied by (ref), we should choose $\tau \asymp \log(d)$. Moreover, the first term of the maximum in (ref) is largest for $1/ \Delta_{n} \geq ( \log(d)+ \tau) \nu_{c,2}/ \nu_{c, \infty}$. In light of these observations, the following corollary is an immediate consequence of Theorem (ref), so we exclude the proof.
It follows from (ref) 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:15a, we can also often assume that $\lVert c_{s} \rVert_{ \infty}$ and, hence, $\nu_{c, \infty}$ can be chosen independently of $d$. When this is the case, the rate implied by (ref) is $\operatorname{rank}( \Sigma) \Delta_{n}d \log(d)$, and since $\operatorname{rank}( \Sigma) \leq d$, Corollary (ref) 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 $\operatorname{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 $\lVert \widehat{ \Sigma}_{n} - \Sigma \rVert_{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}$.
We denote by $\mathcal{C}_{r}$, for a given non-zero integer $r \leq d$, the subclass of $d \times d$ symmetric positive semidefinite matrices $\mathcal{S}_{+}^{d}$ whose effective rank is bounded by $r$:
where
Compared to the rank, the effective rank is a more stable measure for the intrinsic dimension of a covariance matrix (see, e.g., vershynin:10a). In the following we argue that $\widehat{ \Sigma}_{n}^{ \lambda}$, with $\lambda$ of the form (ref), is a minimax optimal estimator of $\Sigma$ in the parameters $\nu_{c,2}$, $\nu_{c, \infty}$, and $\operatorname{rank}( \Sigma)$ over the parametric class of continuous It\^{o} processes generated by (ref) 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 (ref) with no drift and constant volatility $c_{s} \equiv A$. In this setting, we can choose $\nu_{c,2} = \operatorname{tr}(A)$ and $\nu_{c,\infty} = \lVert A \rVert_{\infty}$, which by Corollary (ref) means that
for an absolute constant $\overline{ \gamma}$ and $\Delta_{n}^{-1} \geq 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 \citet*{lounici:14a}. It shows that, up to a logarithmic factor, no estimator can do better than (ref).
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 Section (ref), we begin by stating a general non-probabilistic inequality in Theorem (ref). In the formulation of this result, we recall that $\operatorname{rank}(A; \varepsilon)$ is the number of singular values of $A$ exceeding $\varepsilon \in [0, \infty)$.
With this result in hand we can rely on the exponential inequality for the quantity $\lVert \widehat{ \Sigma}_{n} - \Sigma \rVert_{ \infty}$ established in Theorem (ref) 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$.
In the setting of Corollary (ref), it follows that with high probability an eigenvalue $s$ of $\Sigma$ affects $\operatorname{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, $\operatorname{rank}( \widehat{ \Sigma}_{n}^{ \lambda}) = \operatorname{rank}( \Sigma)$ with a probability tending to one as $n \to \infty$.
The following example illustrates a stylized setting, where many eigenvalues of $\Sigma$ are negligible.
In view of Theorem (ref) and Corollary (ref), it follows that $\widehat{ \Sigma}_{n}^{ \lambda}$ can be of low rank and accurate, given optimal tuning of $\lambda$. However, as evident from (ref), $\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 (ref) with $\widehat{ \Sigma}_{n}$ replaced by $\widehat{ \Sigma}_{n} - \Sigma$). By Proposition (ref), this means $\lambda = 2 \lVert \widehat{ \Sigma}_{n} - \Sigma \rVert_{ \infty}$.
The above is unobservable. To circumvent this problem and facilitate the calculation of our estimator in Section (ref) and (ref), we adopt a data-driven shrinkage selection by exploiting the subsampling technique of \citet*{christensen-podolskij-thamrongrat-veliyev:17a}, see also \citet*{politis-romano-wolf:99a} and kalnina:11a. 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 $l$th subsample consists of the increments $\bigl( \Delta_{(k - 1)L + l}^{n} Y \bigr)_{k = 1, \ldots, n/L}$ for $l = 1, \dots, L$. We denote the associated subsampled RV by
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 \citet*{christensen-podolskij-thamrongrat-veliyev:17a} 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 \citet*[][Theorem 1]{barndorff-nielsen-shephard:04a}.
Hence, in each subsample $\lambda = 2 \lVert \widehat{\Sigma}_{n} - \Sigma \rVert_{ \infty}$ can be approximated by
As an estimator of $\lambda$, we therefore take the sample average:
As noted in Section (ref), 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 $\operatorname{rank} (\Sigma) = d$. On the other hand, there are many situations where it is reasonable to expect that $\operatorname{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 $\operatorname{rank} (c_t)$.
To estimate $c_{t}$, we follow \citet*{jacod-protter:12a} and apply a localized realized variance, which is defined over a block of $\lfloor h_{n} / \Delta_{n} \rfloor \geq 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:
The corresponding penalized version $\widehat{ \Sigma}_{n}^{ \lambda} (t;t+h_{n})$ is computed as
or by using Proposition (ref) 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
Also, we write
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 \colon [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}, \lVert \: \cdot \: \rVert_{2})$ is defined as:
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)$.
The bound on the estimation error of $\hat{c}_{n}^{ \lambda}(t)$ builds on Corollary (ref), but Theorem (ref) 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 (ref). For example, if $(c_{t})_{t \in[0,1]}$ is $1/2$-H\"{o}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\"{o}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, (ref) 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., \citet*{jacod-protter:12a}.
The next result, which relies on Corollary (ref), shows that the rank and performance of $\hat{c}_{n}^{ \lambda}(t)$ can be controlled simultaneously.
Consider the setting of Theorem (ref). For the upper bound in (ref) to be useful, we must have that $\underline{ \varepsilon} > 0$ or, equivalently,
Suppose that $\nu_{c,2}$, $\nu_{c,\infty}$, and $\nu_{c,\psi}$ do not depend on $n$. The inequality (ref) 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 $\lVert \hat{c}_{n}^{ \lambda}(t) - c_{t} \rVert_{2}^{2}$. In that case, (ref) 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 $\operatorname{rank}( \Sigma (t- \frac{h_{n}}{2}; t+ \frac{3h_{n}}{2}))$ in the upper bound of (ref) be replaced by $\operatorname{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.
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 (ref), 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 \citet*{ait-sahalia-xiu:19b}:
where $F_{t} \in \mathbb{R}^{r}$ is a vector of systematic risk factors, which has the dynamic:
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} = \operatorname{diag} ( \sigma_{1,t}, \sigma_{2,t}, \sigma_{3,t}) \in \mathbb{R}^{r \times r}$, where $\operatorname{diag} ( \: \cdot\:)$ is a diagonal matrix with coordinates $\:\cdot\:$, is simulated as a \citet*{heston:93a} square-root process:
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 \citet*{ait-sahalia-xiu:19b}, 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$:
The idiosyncratic component $Z_{t} \in \mathbb{R}^{d}$ is given by
where $g_{t} = \operatorname{diag} ( \gamma_{1,t}^2, \dots, \gamma_{d,t}^2)$ with
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.
The above setup implies
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 \citet*{ait-sahalia-xiu:19b}, 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 (ref) in Section (ref). 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 Figure (ref), 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 Section (ref) with $L = 6$ to choose $\lambda^{*}$.\footnote{We also employed $L = 13$ subsamples, but there were no discernible change in the results.}
In Panel A of Figure (ref), we report the distribution of the $\lambda^{ \ast}$ 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.
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.
We look at high-frequency data from the 30 members of the Dow Jones Industrial Average (DJIA) index as of April 6, 2020.\footnote{On April 6, 2020, Raytheon Technologies (RTX) replaced United Technologies (UTX) in the DJIA index following a merger of United Technologies and Raytheon Company. Moreover, on April 2, 2019, Dow replaced DowDuPont following a spin-off from the parent company, which itself was a fusion between Dow Chemical Company and DuPont in 2017. We employ high-frequency data for the preceding member prior to each index update.} The ticker codes of the various firms are listed in Table (ref), along with selected descriptive statistics. As readily seen from Table (ref), 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{We truncate 5-minute log-returns that exceed three local standard deviations, as gauged by the daily bipower variation estimator, in order to get shelter from potential jumps.} 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 (ref) 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 Section (ref).
In Figure (ref), 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 \citet*{ait-sahalia-xiu:19b}, 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.
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 Figure (ref), 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. (ref)), 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.
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.