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.
94,534 characters · 17 sections · 92 citation commands
Realized Regularized Regressions
\setlist{noitemsep} \onehalfspacing
In recent decades, the empirical asset pricing literature has discovered hundreds of risk factors intended to explain variation in expected returns, i.e., the so-called “factor zoo” cochrane2011presidential,harvey2016and. The increased availability of high-frequency intraday data for both individual assets and aggregate factors (see, e.g., aleti2023high) has further enhanced opportunities for more granular analysis of time variation in factor loadings, or betas, by allowing much shorter estimation windows ait2020high,ait2025continuous.\footnote{Some recent studies apply high-frequency factors to return prediction aleti2025intraday, stochastic discount factor (SDF) estimation aleti2025news, and volatility forecasting cinquetti2025volatility.} Such advances have substantially increased the complexity of the model selection problem from both theoretical and practical perspectives. In this paper, we develop a rigorous methodology for model selection and factor identification in high-frequency regressions with many candidate covariates, while simultaneously retaining valid estimation and inference for the time-varying loadings of the relevant factors.
The standard approach to model selection in static regression models typically relies on penalized estimation, where a regularization term, such as the least absolute shrinkage and selection operator (LASSO; tibshirani1996regression), is added to the objective function. However, in regression models with time-varying coefficients, model selection is considerably less straightforward. A simple strategy is to apply penalized regression separately within each local estimation window, but this can lead to inconsistent outcomes in which the set of selected regressors changes across (possibly overlapping) estimation windows. This issue is especially acute for high-frequency regressions, where the window length usually shrinks to zero under infill asymptotics.
To address this issue, we propose a new approach to continuous-time nonparametric regression that is better suited to “global” model selection over a fixed time interval. Our approach approximates each coefficient path by polynomial splines, and therefore reduces the time-varying beta estimation to the estimation of scalar coefficients in their spline basis expansions. With time variation absorbed by the spline basis functions, it allows us to fit the model once over the entire observation interval, rather than repeatedly re-estimating local regressions. The spline basis expansion also induces a natural group structure: each regressor corresponds to a group of spline coefficients. This enables group-wise model selection to identify the relevant regressors (or equivalently, the groups with nonzero expansion coefficients).
When the number of regressors is fixed, we establish consistency of the time-varying beta estimator, and derive a feasible asymptotic distribution for the integrated beta estimator under infill asymptotics. In high-dimensional settings where the number of irrelevant regressors is allowed to diverge, we show that the penalized spline-based estimator with a truncated $\ell_{1}$-penalty (TLP; shen2012likelihood) attains the “oracle property,” i.e., it selects the correct set of relevant regressors with probability approaching one and estimates their coefficients as if the true model were known a priori fan2001variable.
We investigate the finite-sample performance of the proposed approach through extensive Monte Carlo experiments based on simulated continuous-time factor models with time-varying betas and jumps. In low-dimensional cases, our spline-based estimator performs as well as the benchmark estimator of ait2020high. As the cross-section of candidate factors grows, our estimator remains robust. Our penalized estimation procedure delivers accurate model selection with limited false discoveries, while preserving stable finite-sample estimation performance for the relevant factor exposures.
In the empirical analysis, we illustrate our approach by applying it to the high-frequency “factor zoo” of aleti2023high. We find persistently sparse factor structures in every month across a cross-section of individual stocks and industry portfolios. The market factor emerges as the most stable driver of intraday comovement. A case study of risk decomposition further suggests that a factor set constructed from the most frequently selected factors in our sample achieves greater explanatory power for integrated variance than the standard six-factor benchmark.
Our paper contributes to three strands of the literature. First and foremost, it contributes to the financial econometrics literature on nonparametric regression with high-frequency data. Early studies in this area develop the realized beta estimator as the ratio of realized covariance to realized variance barndorff2004econometric,andersen2005framework,andersen2006realized. mykland2009inference propose the integrated beta estimator by aggregating local coefficient estimates. ait2020high further develop this approach with local ordinary least squares (OLS) and establish feasible limit theorems in the presence of jumps. More recent studies incorporate local regularization to accommodate high-dimensional covariates chen2024realized,shin2025robust,kim2026high. Our contribution to this literature is twofold. First, we propose a spline-based estimator for time-varying betas and establish the associated asymptotic theory. Second, we introduce a method for simultaneous coefficient estimation and model selection in high dimensions which, to the best of our knowledge, is the first to address this issue formally in the context of high-frequency regressions. The closest related work is chen2026high, who develop high-dimensional coefficient tests for the null hypothesis that an additional block of regressors provides no incremental explanatory power. Other notable contributions in this field include, among others: (i) separate inference on factor structure, betas, and risk premia for continuous and jump components todorov2010jumps,bollerslev2016roughing,li2017jump,li2019jump,ait2025continuous, and (ii) inference on time variation in betas patton2012does,reiss2015nonparametric,kong2018testing,kalnina2023inference, and their cross-sectional distribution andersen2021recalcitrant,andersen2023intraday.
Second, we contribute to the rapidly growing literature on high-dimensional problems and statistical learning in finance. As modern financial markets feature a vast set of assets, characteristics, and candidate factors, the penalized regression methods have been widely applied to (i) covariance/correlation matrix estimation and forecasting brownlees2018realized,christensen2023high,bollerslev2025forecasting, (ii) return prediction with large predictor sets chinco2019sparse,gu2020empirical,ait2025predictable,aleti2025intraday, (iii) factor/characteristic selection and SDF estimation for the cross-section of expected returns feng2020taming,freyberger2020dissecting,kozak2020shrinking,chen2025cross, and (iv) implementable portfolio formation and asset allocation ao2019approaching,ding2024statistical, among others.
Third, our work contributes to the statistics literature on varying-coefficient models. Varying-coefficient models extend classical linear regression by allowing coefficients to vary over time or with other variables cleveland1991local,hastie1993varying,fan1999statistical, which provide a flexible and interpretable way to capture dynamic covariate effects, and have been widely adopted in longitudinal studies hoover1998nonparametric,huang2002varying,huang2004polynomial. Work on model selection in varying-coefficient models has closely tracked advances in penalized regression. For example, wang2008variable and wang2009shrinkage adopt the smoothly clipped absolute deviation (SCAD; fan2001variable) and the adaptive LASSO zou2006adaptive, respectively; wei2011variable combine spline-based estimation with the group LASSO yuan2006model in high-dimensional settings; and xue2012variable employ the nonconvex TLP to avoid reliance on a consistent initial estimator, which can be problematic for the SCAD and LASSO approaches, as the irrepresentable conditions for the initial LASSO estimator are often implausible in high-dimensional settings zhao2006model,bach2008consistency. Our contribution to this literature is to establish an infill asymptotic theory for spline-based estimation and model selection in a nonstandard high-frequency setting. Moreover, our asymptotic framework allows for substantially weaker assumptions on the time-varying coefficients. Whereas the classical varying-coefficient literature typically requires smoothness conditions on the coefficient functions, we allow the coefficients to evolve as sample paths of It{\^o} semimartingales with infinite-activity jumps, and hence to be substantially rougher.
The rest of this paper is organized as follows: In (ref), we present the continuous-time regression model, and introduce the required notation and assumptions. In (ref), we propose the estimation procedure and establish its theoretical properties. Monte Carlo experiments and empirical applications are reported in (ref), respectively. (ref) concludes. Proofs and supplementary results can be found in the \hyperref[AP:Proofs]{Appendix}.
We follow the framework of ait2020high to specify a continuous-time multiple regression model and state the assumptions. Consider the following nonparametric regression model,
where $Y=(Y_{t})_{t\geq0}$ is the response process, $X=(X_{1,t},\dots,X_{p,t})^{\top}_{t\geq0}$ is a $p$-dimensional covariate process, and $Z=(Z_{t})_{t\geq0}$ is the residual process. We denote by $X^{c}$ the continuous component of $X$, and by $\Delta X_{t}=X_{t}-X_{t-}$ the jump of $X$ at time $t$. $\beta=(\beta_{1,t},\dots,\beta_{p,t})_{t\geq0}^{\top}$ and $\beta^{J}=(\beta_{1,t}^{J},\dots,\beta_{p,t}^{J})_{t\geq0}^{\top}$ collect the time-varying coefficients with respect to the continuous and discontinuous parts of $X$, respectively.
We assume that both $X$ and $Z$ are It{\^o} semimartingales defined on a filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq0},\mathbb{P})$:
where $X_{0}\in\mathbb{R}^{p}$ and $Z_{0}\in\mathbb{R}$ are $\mathcal{F}_{0}$-measurable, $W$ and $\widetilde{W}$ denote a $p'$-dimensional and a one-dimensional Brownian motion, respectively, the spot volatility $\sigma$ takes values on $\mathbb{R}^{p}\otimes\mathbb{R}^{p'}$, $\mu$ (resp. $\widetilde{\mu}$) is a Poisson random measure with a compensator $\nu$ (resp. $\widetilde{\nu}$) of the form $\nu(dt,dx)=dt\otimes\lambda(dx)$ for some $\sigma$-finite measure $\lambda$ on $\mathbb{R}^{p}$ (resp. $\widetilde{\lambda}$ on $\mathbb{R}$). The spot covariance of $X$, denoted as $c=\sigma\sigma^{\top}$, follows another It{\^o} semimartingale defined on $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq0},\mathbb{P})$:
where $c_{0}\in\mathbb{R}^{p\times p}$ is $\mathcal{F}_{0}$-measurable, and $\sigma'$ takes values on $\mathbb{R}^{p\times p}\otimes\mathbb{R}^{p'}$. We use the same Brownian motion $W$ and Poisson random measure $\mu$ for both $X$ and $c$ without loss of generality, and this specification accommodates leverage effects as well as co-jumps in price and volatility.
We assume that the processes defined above satisfy the following regularity conditions:
The above conditions are standard in the high-frequency literature. Conditions (i) and (ii) ensure the stochastic integrals are well defined and permit the standard localization arguments. Condition (iii) implies that $c$ is uniformly nondegenerate and locally bounded on compact intervals, which guarantees invertibility and rules out pathological cases of vanishing or explosive volatility. Conditions (v) and (vi) allow for both finite- and infinite-activity jumps, while restricting them to be of finite variation.
To account for the time-varying coefficient, we further assume that $\beta$ follows an It{\^o} semimartingale with finite-variation jumps:
Similar to (ref), we impose the following regularity conditions:
Analogous to the classical regression framework, we require an exogeneity assumption over the fixed time interval $[0,T]$. This assumption is specified separately for the continuous and jump components: the former follows barndorff2004econometric and mykland2006anova,mykland2009inference, whereas the latter is in line with li2017jump.
Both $Y$ and $X$ are observed at discrete times $i\Delta_{n}$ for $0\leq i\leq n=\lfloor T/\Delta_{n}\rfloor$. The increments of a generic $p$-dimensional process $A$ are denoted by
For the asymptotic theory below we work in an infill framework, i.e., $n\to\infty$ as $\Delta_{n}\to0$, while $T$ is fixed. Based on the discrete observations of both $Y$ and $X$, our goal is to estimate the time-varying coefficients $\beta$ for the continuous component of $X$. We define the integrated beta ($I\beta$) over $[0,T]$ as
When the dimension $p$ is fixed, $I\beta_{T}$ can be estimated by aggregating local estimates of $\beta$ over $[0,T]$ mykland2009inference,ait2020high. In the fixed-$p$ setting, the local OLS estimator of ait2020high is semiparametrically efficient jacod2013quarticity,renault2017efficient, li2021efficient. However, our main focus is the case where $p \to \infty$. In particular, we assume that the dimension of $X$ increases but the corresponding $\beta$ remains sparse, i.e., a large number of covariates are redundant. In such a high-dimensional setting, the local OLS estimation becomes ill-posed: the local design matrix is nearly singular or not of full column rank, so the local OLS estimator is either not uniquely defined or numerically unstable. For this reason, we seek an alternative estimation procedure whose performance is comparable to that of ait2020high when $p$ is fixed, but which also remains valid as $p$ increases and enables automatic model selection.
One could accommodate the large $p$-case by introducing regularization into each local window. For example, chen2024realized employ the thresholding technique of bickel2008covariance to estimate large spot covariance matrices and to stabilize local regressions in the presence of many regressors; see also, e.g., shin2025robust and kim2026high. Similar ideas could be used to perform model selection in a purely local sense, but such an approach treats model selection as intrinsically timewise: the active set of regressors is allowed to fluctuate across windows, which obscures the economic interpretation of the resulting factor structure.
These considerations motivate our work. The question we consider in this paper is whether one can simultaneously (i) retain the local estimation of $\beta$ (and the estimation of $I\beta$), and (ii) perform global model selection over the entire fixed time interval $[0,T]$.
In this section, we present the proposed estimation method and develop its theoretical properties. (ref) outlines the proposed estimator, with technical details deferred to the subsequent subsections. (ref) derives the asymptotic rate of the spline approximation error. (ref) establishes consistency and derives a feasible asymptotic distribution for our estimators. (ref) extends the framework to high-dimensional settings with group-wise regularization and establishes model selection consistency. (ref) discusses the nonconvex optimization algorithm and the practical choice of tuning parameters.
Our estimation method is based on approximating the coefficient process $\beta=(\beta_{1},\dots,\beta_{p})^{\top}$ on $[0,T]$ by a linear combination of $K_n$ deterministic polynomial spline basis functions with random weights, where $K_n \to \infty$ as $\Delta_n \to 0$:
where $(B_{t}^{(1)},\dots, B_{t}^{(K_{n})})$ collects the spline basis functions of $t$, and $\gamma^{(k)}=(\gamma_{1}^{(k)},\dots,\gamma_{p}^{(k)})^{\top}$ are scalar coefficients in the expansion, which are $\mathcal{F}_{T}$-measurable random variables. Pathwise, for each fixed $\omega\in\Omega$, the expansion coefficients are uniquely determined by the realized path $t\mapsto\beta_t(\omega)$ on $[0,T]$. Specifically, we employ a B-spline basis (see Appendix (ref) for details). This leads to the following approximation for the regression model in (ref):
We collect the spline coefficients in the vector $\bm{\gamma}=(\gamma_{1}^{(1)},\ldots,\gamma_{1}^{(K_n)},\ldots,\gamma_{p}^{(1)},\ldots,\gamma_{p}^{(K_n)})^{\top}\in\mathbb{R}^{pK_n}$. In the absence of jumps on $[0,T]$, we can discretize (ref) in terms of high-frequency increments as
which is a static regression model with dependent variable $\Delta_{i}^{n}Y $ and regressors $B_{(i-1)\Delta_{n}}^{(k)}\Delta_{j,i}^{n}X$ for $j=1,\dots,p$ and $k=1,\dots,K_{n}$. When $p$ is fixed, the spline coefficients in $\bm{\gamma}$ can be estimated by OLS.
In practice, jumps in $X$ or $Z$ may occur on $[0,T]$. To mitigate their impact, we employ the truncation technique of mancini2009non for both $\Delta^n_i Y$ and $\Delta_i^n X$. We then define $\widehat{\bm{\gamma}}$ as the estimator of $\bm{\gamma}$ obtained by minimizing the following truncated least squares criterion:
where $u_{n}>0$ is a truncation threshold satisfying $u_{n}\asymp\Delta_{n}^{\varpi}$ for some $0<\varpi<1/2$. The same truncation threshold $u_{n}>0$ is used for both $\Delta_{i}^{n}Y$ and each component of $\Delta_{i}^{n}X$ for ease of notation, while these thresholds can be different.\footnote{We use the norm truncation $\Vert\Delta_{i}^{n}X\Vert\leq u_{n}$ for notational convenience, as is standard in high-frequency regression literature reiss2015nonparametric,ait2020high. When $p$ is fixed, it is asymptotically equivalent to component-wise truncation. In the high-dimensional analysis in (ref), we instead adopt the component-wise notation $\vert\Delta_{j,i}^{n}X\vert\leq u_{n}$. }
Given $\widehat{\bm{\gamma}}$, the corresponding estimator of $\beta_{t}=(\beta_{1,t},\dots,\beta_{p,t})^{\top}$ takes the form:
The resulting estimator of $I\beta$ is defined as the Riemann sum of $\widehat{\beta}$ on $i\Delta_{n}$ mykland2009inference:
Our proposed estimators $\widehat{\bm{\gamma}}$ and $\widehat{\beta}$ admit closed-form representations in matrix notation. We collect truncated returns of $Y$ and $X$ into the vector $\mathbf{Y}$ and matrix $\mathbf{X}$, respectively:
Define a block-diagonal matrix that collects all B-spline basis functions at time $t$ as
and the design matrix $\mathbf{R}=(\mathbf{R}_{1},\dots,\mathbf{R}_{n})^{\top}\in\mathbb{R}^{n\times pK_{n}}$ with
Then, our estimator $\widehat{\bm{\gamma}}$ takes the following form:
and the spline-based estimator of $\beta_{t}=(\beta_{1,t},\dots,\beta_{p,t})^{\top}$ can be expressed as
We begin by defining a spline approximation to the time-varying coefficient process $\beta=(\beta_{1},\dots,\beta_{p})^{\top}$ on $[0,T]$, where each component is approximated by a polynomial spline of degree $d$. A polynomial spline of degree $d$ on a knot sequence $0=v_{0}<v_{1}<\dots<v_{N_{n}+1}=T$ is a function that coincides with a polynomial of degree $d$ on each subinterval $[v_{i-1},v_{i})$ for $1\le i\le N_{n}$ and on $[v_{N_{n}},v_{N_{n}+1}]$, and that has $d-1$ continuous derivatives when $d\ge 1$. In particular, we consider a quasi-uniform sequence of partitions satisfying
The collection of such spline functions, for a given degree and knot sequence, forms a scalar linear space $\mathcal{G}_{n}$ with a B-spline basis $\{B^{(k)}_{\cdot}\}_{k=1}^{K_{n}}$ of dimension $K_{n}=N_{n}+d+1$; see de1978practical and schumaker2007spline for further details. We then define the vector-valued spline space
Rather than adopting an arbitrary spline approximation of $\beta$, we work with a uniquely defined approximation under a Hilbert-space metric that is intrinsic to the continuous-time regression model. We equip the space of $\mathbb{R}^p$-valued measurable functions on $[0,T]$ with the weighted inner product
where the spot covariance process $c$ is defined in (ref), and write $\Vert f\Vert_{L^{2}(c)}=\langle f,f \rangle_{c}^{1/2}$. We denote by $L^{2}(c)=\{f:\Vert f\Vert_{L^{2}(c)}<\infty\}$ the corresponding Hilbert space. For each $\omega\in\Omega$, let $\Pi_{n}\beta(\omega)$ denote the $L^{2}(c)$-orthogonal projection of the path $t\mapsto \beta_t(\omega)$ onto the spline space $\mathbb{G}_{n}$, i.e., the unique $g\in\mathbb{G}_{n}$ that minimizes
Since, for each $n$, $\mathbb{G}_n$ is a finite-dimensional subspace of $L^2(c)$, the minimizer exists and is unique. In other words, there exists a unique coefficient vector $\bm{\gamma}(\omega)\in\mathbb{R}^{pK_n}$ such that $\Pi_{n}\beta_{t}(\omega) = \mathbf{B}_{t}\bm{\gamma}(\omega)$ for all $0\leq t\leq T$. Consequently, any coefficient path admits the decomposition $\beta_{t}(\omega)=\mathbf{B}_{t}\bm{\gamma}(\omega)+e_{t}(\omega)$, where the spline approximation error $e$ is defined pathwise, and therefore $\widehat{\beta}_{t}$ in (ref) estimates the projected coefficient path $\Pi_{n}\beta_{t}(\omega)$ with the standard OLS estimator $\widehat{\bm{\gamma}}$.
When $\beta$ follows a discontinuous It{\^o} semimartingale in (ref), the spline approximation error satisfies the following asymptotic order:
We establish consistency of our spline-based estimator and derive the corresponding limit theorem for $\widehat{I\beta}$ when $p$ is fixed. Following the literature on varying-coefficient models (see, e.g., hoover1998nonparametric,huang2002varying,huang2004polynomial), we say that $\widehat{\beta}$ is a consistent estimator of $\beta$ if $\Vert\widehat{\beta}-\beta\Vert_{L^{2}}\overset{\mathbb{P}}{\longrightarrow}0$. We denote by $\asymp$ the same order of magnitude, i.e., $f\asymp g$ indicates that $K|g|\leq|f|\leq K'|g|$ for some constants $K,K'>0$.
(ref) shows that the consistency of $\widehat{\beta}$ holds with the standard choice of truncation threshold from mancini2009non under the mild condition $\Delta_{n}K_{n}\log K_{n}\to0$. The identifiability of $\widehat{\bm{\gamma}}$, and hence of $\widehat{\beta}$, follows from the nonsingularity of the empirical Gram matrix $\mathbf{R}^{\top}\mathbf{R}$, which is implied by that of its theoretical counterpart huang2003local, and, in turn, is ensured by the boundedness of B-spline basis functions de1978practical.
We now present a central limit theorem for $\widehat{I\beta}$. We write $X^{n}\xrightarrow{\mathcal{L}-s}X$ to denote stable convergence in law, i.e., $(X^{n},Y)\overset{\mathcal{L}}{\longrightarrow}(X,Y)$ for any $\mathcal{F}$-measurable process $Y$. We denote by $\mathcal{MN}$ a mixed normal distribution, i.e., a conditionally normal distribution with random $\mathcal{F}$-conditional variance.
The above stable central limit theorem is established by decomposing the normalized estimation error into a leading martingale term and a collection of remainder terms. The latter arise from spline approximation and discretization errors in the construction of $\widehat{I\beta}$, as well as from jumps in the covariate and residual processes. The leading term admits a truncated realized-covariance representation and converges stably to a mixed normal limit. This result requires stronger conditions than Theorem (ref) on the truncation threshold $u_n$, as is standard in the asymptotic theory of truncated realized covariance (see, e.g., Theorem 13.2.1, jacod2012discretization). The condition $\Delta_{n}K_{n}\log K_{n}\to0$ follows from (ref) and is analogous to standard conditions in the spline-based varying-coefficient model literature (see, e.g., Theorem 3, huang2004polynomial). The additional condition $K_{n}\sqrt{\Delta_{n}}/\log K_{n}\to\infty$ ensures that the bias induced by spline approximation is asymptotically negligible and thus has no impact on the central limit theorem.
For feasible implementation of the asymptotic distribution in (ref), we can estimate $\Sigma_{T}^{\beta}$ by a jackknife White's heteroskedasticity-consistent estimator mackinnon1985some:
In $\widehat{\Sigma}_{T}^{\beta}$ we use leave-one-out residuals to avoid reusing the same increments for the estimation of both coefficients and covariance matrix, which restores the conditional orthogonality between the coefficient estimate $\widehat{\beta}_{(i-1)\Delta_{n}}^{(-i)}$ and the $i$-th interval increments needed for the covariance matrix estimation consistency.\footnote{For consistency of the standard sandwich estimator based on the full-sample $\widehat{\beta}$, one needs additional uniform control of negligible leverage across all $i$, which becomes restrictive when the regression dimension $pK_n$ diverges; see, e.g., kline2020leave, for related discussions. } For computational convenience, the leave-one-out residuals can be replaced by a block cross-fitted analogue in the spirit of chernozhukov2018double, which requires only a finite number of refits. We leave this extension to future work.
In this section, we assume that the $p$-dimensional coefficient process takes the following form:
where $q\ge1$ is the number of relevant regressors and $0<\Vert\beta_{j}\Vert_{L^{2}}<\infty$ for $j=1,\dots,q$. That is, the first $q$ regressors are relevant with nonzero coefficients, while the remaining $p-q$ regressors are redundant and have coefficients that are almost surely zero for all $t\in[0,T]$.\footnote{We assume $\beta^{J}$ has the same sparsity structure as $\beta$, i.e., $\beta_{t}^{J}=(\beta_{1,t}^{J},\dots,\beta_{q,t}^{J}, 0,\dots,0)^{\top}$. } In line with the model selection literature, we consider a high-dimensional regime in which $p=p_{n}\to\infty$, while the number of relevant regressors $q$ remains fixed.
Our goal is to simultaneously estimate the coefficients of relevant regressors and perform the model selection. For each irrelevant regressor, the associated B-spline coefficient vector $\gamma_j=(\gamma_j^{(1)},\dots,\gamma_j^{(K_n)})^\top$ is identically zero. Consequently, consistent model selection can be achieved by imposing a group penalty on $\gamma_{j}$ for $j=1,\dots,p$.
We consider an estimator $\widehat{\bm{\gamma}}^{*}$ defined as a global minimizer of
where $\rho_{n}(x)$ is the truncated $\ell_{1}$-penalty (TLP) function with a thresholding parameter $\tau_n>0$ shen2012likelihood, and $\mathbf{W}_{j}=(\mathbf{R}_{j}^{*})^{\top}\mathbf{R}_{j}^{*}$ is the empirical Gram matrix of block $j$, with $\mathbf{R}_{j}^{*}$ the $n\times K_{n}$ submatrix of $\mathbf{R}$ corresponding to the $j$-th regressor:
The objective function $Q_{n}^{*}(\bm{\gamma})$ adopts a general form for group selection yuan2006model,huang2012selective, while we select the nonconvex TLP to obtain model selection consistency, which cannot be achieved with standard group LASSO without irrepresentable conditions bach2008consistency. The TLP function is flat beyond the threshold $\tau_n$, so it penalizes small coefficients strongly while leaving sufficiently large true signals unpenalized.
Let $\widehat{\bm{\gamma}}^{\circ}$ be the oracle estimator, i.e., the restricted minimizer of $Q_{n}(\bm{\gamma})$ subject to $\gamma_{j}\equiv0$ for $q<j\leq p$. Equivalently, $\widehat{\bm{\gamma}}^{\circ}$ is the least-squares spline estimator for the model that includes only the $q$ relevant regressors fan2001variable. (ref) shows that $\widehat{\bm{\gamma}}^{\circ}$ is a local minimizer of the penalized objective $Q_{n}^{*}(\bm{\gamma})$, and (ref) provides sufficient conditions such that the global minimizer $\widehat{\bm{\gamma}}^{*}$ coincides with $\widehat{\bm{\gamma}}^{\circ}$, which establishes the model selection consistency.
(ref) is established by verifying the standard Karush-Kuhn-Tucker (KKT) local optimality conditions. The key requirement is that, uniformly over all irrelevant regressors, the effective penalty level $\lambda_{n}/\tau_{n}$ dominates the blockwise gradient of least-squares loss—a standard calibration for both convex and nonconvex penalties for model selection. The first two terms in the denominator of (ref) follow from a Bernstein-type maximal inequality for high-dimensional martingale difference arrays (as an analogue to Lemma A.1, van2008high), combined with a uniform control of the predictable compensators in the third term. Altogether, the condition ensures that irrelevant blocks are shrunk to zero while the active blocks fall into the flat region of TLP and are therefore left unpenalized, such that the local KKT solution coincides with the oracle estimator $\widehat{\bm{\gamma}}^{\circ}$.
(ref) follows from the fact that any misspecified model—whether it omits relevant groups, includes redundant ones, or both—attains a strictly larger penalized objective than the oracle model with probability approaching one. The additional conditions are calibrated to control two types of errors in opposite directions. On the one hand, letting $\lambda_{n}\to0$ ensures the penalty does not over-reward parsimony, and therefore the increase in squared loss from omitting relevant groups dominates any penalty reduction. On the other hand, if instead the model includes redundant groups, it may achieve a decrease in squared loss through overfitting, but the second condition in (ref) guarantees that the penalty is sufficiently strong to outweigh that spurious improvement. In short, the conditions balance fidelity and parsimony: omissions are ruled out because loss inflation dominates when $\lambda_{n}$ is small, and overfitted models are also ruled out because the penalty does not vanish too quickly. Moreover, the second condition in (ref) implies $\Delta_{n}(K_{n}+\log p)\to0$, so the dimension $p$ can diverge only at a rate slower than exponential relative to $n$. Extending the result to accommodate more rapidly growing model complexity is left for future work.
In practice, the nonconvex objective $Q_{n}^{*}(\bm{\gamma})$ is minimized by the difference-of-convex (DC) algorithm in the spirit of shen2012likelihood and xue2012variable. Specifically, the TLP function can be written as $\rho_{n}(x)=\rho_{1,n}(x)-\rho_{2,n}(x)$, where
This allows us to decompose $Q_{n}^{*}(\bm{\gamma})=Q_{1,n}^{*}(\bm{\gamma})-Q_{2,n}^{*}(\bm{\gamma})$, where
Both $Q_{1,n}^{*}(\bm{\gamma})$ and $Q_{2,n}^{*}(\bm{\gamma})$ are convex functions and, in particular, $Q_{1,n}^{*}(\bm{\gamma})$ is a standard group LASSO objective with effective penalty level $\lambda_{n}/\tau_{n}$. We initialize the algorithm at $\widehat{\bm{\gamma}}^{(0)}=0$, since the DC procedure does not require a consistent initial estimator. At iteration $m\geq1$, we keep $Q_{1,n}^{*}(\bm{\gamma})$ unchanged but approximate $Q_{2,n}^{*}(\bm{\gamma})$ by its affine minorant at $\widehat{\bm{\gamma}}^{(m-1)}$. Because $\rho_{2,n}(x)$ is piecewise linear with slope $0$ when $0\leq x<\tau_{n}$ and slope $1/\tau_{n}$ when $x>\tau_{n}$, this leads to the convex optimization for $\widehat{\bm{\gamma}}^{(m)}$:
so that blocks with $\sqrt{(\widehat{\gamma}_{j}^{(m-1)})^{\top}\mathbf{W}_{j}\widehat{\gamma}_{j}^{(m-1)}}\leq\tau_{n}$ are penalized with weight $\lambda_{n}/\tau_{n}$, whereas blocks with $\sqrt{(\widehat{\gamma}_{j}^{(m-1)})^{\top}\mathbf{W}_{j}\widehat{\gamma}_{j}^{(m-1)}}>\tau_{n}$ remain unpenalized. Therefore, the DC algorithm generates a sequence $(\widehat{\bm{\gamma}}^{(m)})_{m\ge0}$, where each update $\widehat{\bm{\gamma}}^{(m)}$ is from a standard group LASSO with group-specific tuning parameters $(\lambda_{j,n}^{(m)})_{j=1}^{p}$. In our implementation, this group LASSO is solved by the standard proximal gradient method.
We consider a data-driven approach to selecting the tuning parameters. For the threshold $\tau_{n}$, the standard KKT conditions suggest that it should be small enough to distinguish whether $\sqrt{\gamma_{j}^{\top}\mathbf{W}_{j}\gamma_{j}}$ for each $j$ is essentially zero or not. Intuitively, $\gamma_{j}^{\top}\mathbf{W}_{j}\gamma_{j}$ measures the contribution of regressor $j$ to the integrated variance of $Y$. In our implementation, we set $\tau_{n}$ to be a small multiple, e.g., $\alpha_{\tau}=0.05$ or 0.01, of the square root of the median realized volatility (MedRV, andersen2012jump) of $Y$.
We employ cubic B-splines throughout the simulation and empirical applications in this paper. Conditional on the above choice of $\tau_{n}$, we select the number of B-spline basis functions $K_{n}$ and the TLP parameter $\lambda_{n}$ (or, equivalently, the effective penalty level $\lambda_{n}/\tau_{n}$) by $K$-fold cross-validation over a reasonable grid of candidate values. The cross-validation criterion is the out-of-sample mean squared prediction error for the truncated returns of $Y$. Furthermore, to improve interpretability in finite samples, we can apply the “one-standard-error” rule with cross-validation, i.e., select the largest effective penalty whose cross-validation criterion is no more than one standard error above the minimum to obtain the most parsimonious model (see, e.g., Chapter 7.10, hastie2009elements). We implement this rule in our empirical applications in (ref).
This section contains a Monte Carlo study to evaluate both the estimation and model selection performance of our estimators, which correspond to the theoretical results developed in (ref).
We assess the finite‐sample performance of our estimator following a simulation design similar to that in ait2020high. We consider a continuous-time factor model in (ref) for the log-price $Y_{t}$ of a single asset, with $q=3$ relevant factors out of $p$ factors in total. The covariate process $X_t = (X_{1,t},\dots,X_{p,t})^{\top}$ evolves as
where $b_{j}$ is the drift of factor $j$, $(W_{1,t},\dots,W_{p,t})^{\top}$ is a $p$-dimensional standard Brownian motion, the symmetric matrix with ones on the diagonal and off-diagonal entries $\rho_{jk}$ controls the dependence across factors, $J_{j,t}$ is the jump size of factor $j$ at time $t$, and $N_t$ is a Poisson process with annualized intensity $\lambda$. The factor volatilities follow the Cox-Ingersoll-Ross (CIR) dynamics:
where $\kappa'_j$ is the mean-reversion speed, $\alpha'_j$ is the long-run volatility level, $v'_j$ is the volatility-of-volatility parameter, and $(W'_{j,t})_{j=1}^{p}$ are standard Brownian motions independent of $(W_t, N_t)$. For each $j = 1,\dots,p$, the factor jump sizes $J_{j,t}$ are i.i.d. across jump times and follow a double-exponential distribution:
where $\omega_{j}$ is the mixture probability. The jump sizes in the volatility process, $J'_{j,t}$, are exponentially distributed with mean $g'_{j}$. The Poisson process $N_t$ is common to all factors and their volatilities, so jumps occur simultaneously across components.
The idiosyncratic component $Z_t$ is simulated as a jump-diffusion,
where $\widetilde{W}_{t}$ is a standard Brownian motion independent of $(W_t, W'_t, N_t)$, and the idiosyncratic jump sizes $\widetilde{J}_{t}$ follow another double-exponential distribution:
while $\widetilde{N}_{t}$ is a Poisson process with the same annualized intensity $\lambda$ and is independent of $N_{t}$.
Finally, we simulate the time-varying coefficients $(\beta_{1,t},\beta_{2,t},\beta_{3,t},0,\dots,0)^{\top}$ following an Ornstein-Uhlenbeck process:
where $(B_{1,t},B_{2,t},B_{3,t})^{\top}$ is a three-dimensional standard Brownian motion independent of all previously defined sources of randomness. The response process $Y_{t}$ is then obtained as
The parameters are set as follows:
We simulate 5-minute observations for both $Y$ and $X$ over a month (21 days) for each replication. (ref) illustrates a simulated path of $Y$ together with the three relevant factors, with all initial values set to zero. Most of our parameter choices are calibrated to be consistent with ait2020high in their benchmark case with $p=q=3$. The only deviation is that we set relatively high long-run means for the coefficients of relevant factors in order to avoid a “weak factor” configuration, which would make true signals difficult to distinguish from noise and thus complicate the finite-sample assessment of model selection performance. Under this calibration, the three relevant factors contribute approximately 18.7%, 7.6% and 3.7%, respectively, to the quadratic variation of $Y$.
Following the tuning parameter selection in (ref), for each Monte Carlo experiment, we use the first 50 simulated paths to perform 5-fold cross-validation, and then fix the parameters at the median selected values for all replications. In (ref), we report the finite-sample performance of the $\widehat{I\beta}$ estimator for the relevant factors and the model selection results, for $p$ varying from 3 to 500, which covers both the three-factor benchmark case of ait2020high and a high-dimensional “factor zoo” setting.\footnote{We also consider designs with more pronounced cross-factor dependence, which is common in high-dimensional settings; see Appendix (ref) for further discussion.}
(ref) reports the sample bias, standard deviation and root mean square error (RMSE) of the $\widehat{I\beta}$ estimators for relevant factors. Panel A shows that, in the three-factor benchmark case $p=q=3$, the unpenalized spline-based estimator delivers very accurate estimates for $I\beta$, with small bias and low RMSE, and performs comparably to the estimators of ait2020high. The penalized spline-based estimator yields nearly the same results as the unpenalized one, since all three factors are truly relevant and, with our data-driven truncation threshold $\tau_{n}$, the TLP is effectively inactive in this case. Overall, in this low-dimensional case the spline-based estimator matches the performance of the local OLS benchmark, in line with our asymptotic results for fixed $p=q$.
As the number of candidate factors increases from $p=3$ to 100, the RMSEs of the unpenalized spline-based estimator remain stable and rise only slightly across Panels A–D. Introducing the penalty has very little impact on the estimation error of $\widehat{I\beta}$, which indicates that the TLP regularizes the contribution of irrelevant factors without shrinking the true coefficients of the relevant ones, in line with the main intuition behind TLP and the oracle property in (ref). In particular, for $p=100$, the ait2020high estimator uses a local window size smaller than $p$, which renders the local design matrix ill-conditioned and makes estimation infeasible. By contrast, our spline-based estimator is obtained from a single “global” regression of all truncated increments of $Y$ on those of $X$ over $[0,T]$, and thus it remains numerically stable and more robust than the local estimation procedure when $p$ is large.
The high-dimensional case with $p=500$ (Panel E) further highlights the role of regularization. In this case, the unpenalized spline-based estimator breaks down: the biases of $\widehat{I\beta}$ exceed those in the $p=100$ case by more than two orders of magnitude, and the RMSEs are of the same order, which indicates that OLS is no longer reliable when the number of nuisance factors is very large such that the effective number of coefficients in $\bm{\gamma}$, $pK_{n}$, exceeds the sample size $n$. In sharp contrast, the penalized spline-based estimators retain RMSEs of about 0.03 and small biases, comparable to those observed for $p=10$, 50 and 100. In conclusion, the Monte Carlo evidence in (ref) shows that our spline-based estimator matches the efficiency of the local OLS method in small dimensions but is substantially more robust when the cross-section of candidate factors becomes large, and that the specific nonconvex penalty is crucial to maintaining good finite-sample performance when $p$ is very large.
In this section, we examine the finite-sample model selection performance of the penalized spline-based estimator. Consistent with (ref), “selection” refers to a factor whose corresponding block of B-spline coefficient estimates $\widehat{\gamma}_{j}=(\widehat{\gamma}_{j}^{(1)},\dots,\widehat{\gamma}_{j}^{(K_n)})^\top$ is nonzero.
To illustrate the impact of the tuning parameters, we first consider the case $p=10$ and $q=3$. (ref) reports the true discovery rate (TDR), defined as the fraction of truly relevant factors that are correctly selected, and the false discovery rate (FDR), defined as the fraction of irrelevant factors that are incorrectly included in the active set, over a two-dimensional grid of tuning parameters, i.e., the number of B-spline basis functions $K_{n}$ and the effective penalty level $\lambda_{n}/\tau_{n}$. For small values of $\lambda_{n}/\tau_{n}$, the penalization is weak and both the TDR and FDR are high. As $\lambda_{n}/\tau_{n}$ increases, the FDR decreases sharply and the procedure becomes increasingly conservative, until, for very large penalties, the TDR starts to deteriorate because some relevant factors are also shrunk out. Hence, the “good” region of the grid is characterized by combinations of $(\lambda_{n}/\tau_{n},K_{n})$ for which the TDR is essentially one and the FDR is close to zero. For this example, the values of $K_{n}=4$ and $\lambda_{n}/\tau_{n}=0.07$ selected by 5-fold cross-validation fall into this region.
(ref) reports the empirical selection frequencies of each factor in the case with $p=10$ and $q=3$, with two choices of the truncation constant $\alpha_{\tau}=0.05$ and 0.01, and corresponding cross-validated $K_{n}$ and $\lambda_{n}/\tau_{n}$. The three relevant factors are selected in essentially all replications, while the irrelevant factors are selected only rarely when $\alpha_{\tau}=0.05$ and nearly never when $\alpha_{\tau}=0.01$. This pattern is fully in line with our model selection consistency in (ref).
(ref) summarizes the model selection results across all cases considered in the previous section, which reports the average selection frequencies for relevant and irrelevant factors, and the frequencies of exactly correct model specification, i.e., all three relevant factors selected and all irrelevant ones excluded. For all cases and both values of $\alpha_{\tau}$, the (near-)unity selection frequencies for the relevant factors and the near-zero frequencies for the irrelevant ones, jointly with the estimation results in (ref), indicate that the proposed penalized spline-based estimator simultaneously delivers accurate coefficient estimation and reliable recovery of sparse factor structure, which provides a direct finite-sample illustration of the oracle property established in (ref).
In this section, we apply the methodology developed in this paper to a large number of high-frequency factors and anomalies. We start with a standard six-factor analysis for some representative stock, and compare our high-frequency estimates with their low-frequency counterparts to illustrate the benefits of exploiting granular intraday information. Then we examine whether a relatively small number of high-frequency factors from the “factor zoo” can explain comovement across a large cross-section of liquid U.S. stocks and industry portfolios, where we implement simultaneous coefficient estimation and model selection by employing our penalized spline-based estimation procedure.
We use the high-frequency factor zoo dataset of aleti2023high. This minute-level dataset comprises 272 high-frequency portfolios over 1996--2020: (i) the five fama2015five (FF5) factors plus momentum, constructed following the original Fama-French methodology; (ii) 218 characteristic-sorted factor portfolios based on the predictor libraries of chen2022open and jensen2023there; and (iii) 48 Fama-French industry portfolios fama1997industry. For our empirical analysis, we consider the period from January 2016 to December 2020, and adopt a 5-minute sampling frequency for all 224 factors, individual stocks, and industry portfolios.
For the individual stocks, we start with S&P 100 constituents and retain only firms that are consistently included in the index throughout January 2016 to December 2020. As is standard in empirical research with the Trade and Quote (TAQ) data,\footnote{We use the SAS code from holden2014liquidity to extract all tick-by-tick transaction records matched with relevant ask/bid quotes from the daily TAQ dataset of WRDS. } we use the standard data filters as in barndorff2009realized to eliminate data errors, remove all transactions in the original record that are later corrected, canceled or otherwise invalidated. In addition, we remove all trading days with an early market closure, and restrict our sample to transactions between 9:30:00--16:00:00 Eastern Time (ET). More importantly, we keep only stocks whose 5-minute returns are almost always non-zero (fewer than 5% zero returns in each year). After these screens, our sample retains 57 stocks, which span 9 sectors and 21 Fama-French industry groups.\footnote{Appendix (ref) reports the list of individual stocks retained in our sample, together with their sector and industry group assignments, and the fraction of zero 5-minute returns in each year. Sector classifications are based on the GICS sector code from WRDS Compustat. The Fama-French 48 industry group for each firm is assigned based on its four-digit SIC code from Compustat and the SIC-to-industry definitions from Kenneth R. French's data library. } In addition, we include the full set of 48 Fama-French industry portfolios from aleti2023high, for a total of 105 test assets.
We estimate the betas in a standard six-factor model (FF5 plus momentum) for Apple Inc. (AAPL) with three different approaches: (i) our unpenalized spline-based estimator, (ii) the local OLS estimator of ait2020high, and (iii) a low-frequency OLS estimator based on daily data. Following the Monte Carlo simulations in (ref), we select the number of B-spline basis functions $K_{n}$ via 5-fold cross-validation. For the AKX estimator, we set the local window length to $k_n = 78$.
(ref) reports monthly beta estimates for each factor. For the two high-frequency methods, the plotted series are monthly $\widehat{I\beta}$'s computed from truncated 5-minute returns of both AAPL and the six factors. For the daily-data benchmark, the OLS beta is estimated from daily returns within each month and is therefore constant over that month. We find that the spline-based estimator and the local OLS estimator deliver nearly identical beta estimates across all six factors. Moreover, both high-frequency estimators produce markedly smoother and more stable beta dynamics than the low-frequency benchmark, which exhibits substantially larger month-to-month variation and pronounced spikes. These differences are consistent with the role of jumps. Since large price moves identified as jumps are truncated and do not contribute to $\widehat{I\beta}$, the high-frequency beta estimators primarily reflect the continuous covariation between the asset and the factors. In contrast, the daily OLS estimator is directly exposed to large daily returns when jumps occur, which can induce sizable distortions in monthly beta estimates. Overall, (ref) indicates that high-frequency beta estimators are more robust and deliver more stable beta dynamics, which is consistent with the empirical evidence in ait2020high.
We now apply our penalized spline-based estimator to the high-frequency factor zoo. The goal is to investigate whether the continuous component of each selected stock and industry portfolio can be explained by a relatively small set of factors with time-varying exposures. As summarized in (ref), we work with the 224 high-frequency factors of aleti2023high as candidate covariates, and study a cross-section of liquid S&P 100 stocks together with the Fama-French 48 industry portfolios as test assets. Consistent with the Monte Carlo applications in (ref), we select the tuning parameters, i.e., the number of B-spline basis functions $K_{n}$ and the effective penalty level $\lambda_{n}/\tau_{n}$, via 5-fold cross-validation using the one-standard-error rule. We set the TLP cutoff level $\tau_{n}=0.01\sqrt{\mathrm{MedRV}_{T}}$, where $\mathrm{MedRV}_{T}$ is the MedRV of the test asset computed over the corresponding estimation window $[0,T]$.
(ref) reports the average number of selected factors for a set of representative stocks, i.e., Apple (AAPL), Microsoft (MSFT), JPMorgan Chase (JPM), Johnson & Johnson (JNJ), Exxon Mobil (XOM), and their corresponding Fama-French industry portfolios, for each year from 2016 to 2020. The results suggest that the factor representations fitted by our penalized estimation procedure are generally parsimonious. For most of the representative stocks, the average number of selected factors falls between 15 and 25, whereas the corresponding industry portfolios tend to select fewer factors, with averages around 10. This difference is consistent with the additional diversification at the industry level, which attenuates idiosyncratic variation and concentrates the continuous comovement on a smaller set of common drivers. Furthermore, the results suggest time variation in sparsity: For AAPL and MSFT, the average number of selected factors declines over 2016--2020, while JPM, JNJ, and XOM exhibit an increase over the last two years of the sample. The most pronounced change occurs for the Petroleum and Natural Gas industry portfolio, for which our estimator selects around 30 factors on average in 2020, making it one of the industry portfolios with the largest selected sets in that year. The full results for all selected S&P 100 stocks and industry portfolios are reported in Appendix (ref).
Based on factor selections across the full cross-section of individual stocks and industry portfolios, (ref) lists the twenty factors that are selected most frequently over the period from January 2016 to December 2020. The reported selection rate is the fraction of asset-month pairs for which the factor is selected by our penalized estimation procedure. To facilitate interpretation, we also report broad economic labels from chen2022open and the theme classification from jensen2023there for these frequently selected factors. To assess whether the dominant channels of intraday comovement differ between firm-level and industry-level returns, (ref) reports the same rankings separately for the selected S&P 100 stocks and for the industry portfolios.
Both (ref) show that the market factor is selected for the vast majority of asset-month pairs, which indicates that the market remains the most stable driver of continuous intraday comovement in our sample. Beyond the market factor, however, selection rates decline rapidly: the next-ranked factors are selected in fewer than $15\%$ of asset-months. This pattern suggests that there is no single non-market factor that plays a pervasive role across assets and time. Instead, the non-market component is accounted for by a set of factors that are heterogeneous and episodic, and is therefore well captured by a sparse representation estimated from more granular information within each month. The results also indicate that the most frequently selected factors cluster into a small number of economically coherent themes, including trading activity and liquidity, intangibles and innovation, investment, and balance-sheet-related characteristics. Moreover, the relative importance of these themes differs across firm-level and industry-level analyses. In particular, the stock-level rankings place more weight on trading activity and liquidity-related factors, together with valuation and investment-related measures, whereas the industry-level rankings tilt more toward broader risk-sensitivity and lead-lag-type predictors.
We further conduct a case study of total risk decomposition using two alternative factor sets: the standard six factors (FF5 plus momentum) and the twenty most frequently selected factors from (ref). The analysis parallels the decomposition in ait2020high (Fig. 6) with six factors to quantify how much of AAPL’s total intraday variation is attributable to systematic covariation with a given factor set versus the residual (idiosyncratic) component.
(ref) illustrates the monthly risk decomposition for AAPL from January 2016 to December 2020. Specifically, the (annualized) quadratic variation of AAPL is estimated by realized volatility constructed from 5-minute returns within each month, and the integrated variance is estimated by the truncated realized volatility of mancini2009non with the same truncation threshold used in our regression analysis. For each month and for each factor set, we estimate AAPL’s time-varying betas using our unpenalized spline-based regression on truncated 5-minute increments. Given the fitted beta paths, we calculate the residual increments as the difference between the truncated increments of AAPL and their corresponding predicted values. The unexplained integrated variance is then obtained as the annualized sum of squared residual increments, and the monthly $R^2$ is defined as the fraction of integrated variance explained by the regression model.
Relative to the six-factor benchmark, the specification based on the twenty most frequently selected factors provides a systematically tighter fit. The unexplained integrated variance is smaller, and the corresponding $R^2$ is higher in nearly every month. The results indicate that our cross-sectionally selected factors provide a more effective sparse representation for AAPL’s continuous intraday variation than the standard six-factor specification over this sample period.
This paper develops a high-frequency penalized regression framework for It\^{o} semimartingales with time-varying coefficients. By approximating coefficient paths with polynomial splines, the approach provides a simple “global” alternative to local regressions over rolling windows, and makes high-dimensional model selection feasible with a group-wise $\ell_{1}$-penalty. We establish a comprehensive asymptotic theory for the spline-based estimator, and show that the penalized estimator attains the oracle property when there exists a large number of redundant factors. Monte Carlo experiments corroborate that the proposed procedure simultaneously delivers reliable coefficient estimation and model selection in finite samples. An empirical application to a large collection of high-frequency factors further indicates that continuous intraday comovement is well captured by a sparse and economically coherent set of factors. The selection results underscore the important role of the market factor and reveal a sparse theme structure: The frequently selected factors cluster into a small number of thematic classifications, such as liquidity, investment and valuation, and these themes vary systematically across stock- and industry-level analyses.
\numberwithin{assumption}{section} \numberwithin{equation}{section} \numberwithin{table}{section} \numberwithin{figure}{section} \numberwithin{lemma}{section} \numberwithin{remark}{section} \numberwithin{proposition}{section} \numberwithin{corollary}{section} \numberwithin{definition}{section}