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.
69,005 characters · 14 sections · 61 citation commands
Machine Learning Panel Data Regressions with Heavy-tailed Dependent Data: Theory and Application
{\it Keywords:} High-dimensional panels, large $N$ and $T$ panels, mixed-frequency data, sparse-group LASSO, fat tails. \\ \thispagestyle{empty}
\setcounter{page}{0}
We analyze panel data regressions in a high-dimensional setting where the number of time-varying covariates can be very large and potentially exceed the sample size. We leverage on the structured sparsity approach using sparse-group LASSO (sg-LASSO) regularization for time series data with dictionaries. The advantages of this approach for individual time series data, potentially sampled at mixed frequencies, have been recently reported in babii2020machine, who focus on nowcasting the US GDP growth in a data-rich environment. In this paper, we first show how to leverage on the sparse group regularization in a panel data setting. Second, we study the benefits of using the cross-sectional dimension for prediction with panel data paying particular attention to the issues of fat-tailed series which are relevant for the application involving financial time series. Third, we develop the debiased heteroskedasticity autocorrelation consistent (HAC) inference for regularized panel data regressions. Lastly, we provide an illustrative empirical example involving systematically predictable errors in analysts with individual firm earnings forecasts.
Our paper relates to the literature on high-dimensional panel data models and the (group) LASSO regularization; see harding2019panel, chiang2019post, chernozhukov2019demand, belloni2019high, belloni2016inference, lu2016shrinkage, kock2016oracle, su2016identifying, farrell2015robust, kock2013oracle, lamarche2010robust, koenker2004quantile, among others. However, to the best of our knowledge, the existing literature relates mostly to the microeconometric problems and does not address comprehensively (1) the advantages of long panels; (2) the performance of regularized panel data estimators with potentially heavy-tailed covariates and regression errors, (3) the debiased HAC inference for regularized panel data, and (4) the sg-LASSO regularization of simon2013sparse in a panel data setting.
We recognize that the economic and financial time series data are often persistent with fat tails. To that end, we introduce a new Fuk-Nagaev concentration inequality for long panels. Using this inequality, we obtain oracle inequalities for the sg-LASSO that shed new light on how the predictive performance of pooled and fixed effect estimators scales with $N$ (cross-section) and $T$ (time series), which is especially relevant for modern panel data applications, where both $N$ and $T$ can be large; see fernandez2016individual, hansen2007asymptotic, alvarez2003time, hahn2002asymptotically, and phillips1999linear, among others. Importantly, our theory covers the LASSO and the group-LASSO estimators as special cases of sg-LASSO.
In our empirical application we revisit a topic raised by ball2018automated and carabias2018real, but not resolved via formal inference in a high-dimensional setting. Namely, their empirical findings suggest that analysts tend to focus on their firm/industry when making earnings predictions while not fully taking into account the macroeconomic events affecting their firm/industry. More broadly, ball2018automated argue that analysts do not fully exploit information embedded in high-dimensional data and therefore {\it leave money on the table}. Thanks to the theoretical contributions in the current paper we can formally test that hypothesis in a data-rich environment. Note that, as ball2018automated point out, it is important to take into account the mixed frequency nature of the data flow, which is why the machine learning panel regression methods presented in the paper apply to mixed frequency data. We use 26 predictors, including traditional macro and financial series as well as non-standard series generated by textual analysis of financial news. Using such a rich set of covariates, we test whether analyst' consensus earnings prediction errors are systematically related to either one of the aforementioned variables.
The paper is organized as follows. Section (ref) introduces the models and estimators. Oracle inequalities for sg-LASSO panel data regressions appear in Section (ref). Section (ref) develops the debiased HAC inference for regularized panel data regressions. Monte Carlo simulations are reported in Section (ref). The results of our empirical application are reported in Section (ref). Section (ref) concludes. All technical details and detailed data descriptions appear in the Appendix and the Online Appendix.
\paragraph{Notation:} For a random variable $X\in\ensuremath{\mathbf{R}}$, let $\|X\|_q=(\ensuremath{\mathbb{E}}|X|^q)^{1/q}$ be its $L_q$ norm with $q\geq 1$. For $p\in\ensuremath{\mathbf{N}}$, put $[p] = \{1,2,\dots,p\}$. For a vector $\Delta\in\ensuremath{\mathbf{R}}^p$ and a subset $J\subset [p]$, let $\Delta_J$ be a vector in $\ensuremath{\mathbf{R}}^p$ with the same coordinates as $\Delta$ on $J$ and zero coordinates on $J^c$. Let $\mathcal{G}$ be a partition of $[p]$ defining the group structure, which is assumed to be known to the econometrician. For a vector $\beta\in\ensuremath{\mathbf{R}}^p$, the sparse-group structure is described by a pair $(S_0,\mathcal{G}_0)$, where $S_0=\{j\in[p]:\;\beta_j\ne 0 \}$ and $\mathcal{G}_0 = \left\{G\in\mathcal{G}:\; \beta_{G} \ne 0\right\}$ are the support and respectively the group support of $\beta$.
We also use $|S|$ to denote the cardinality of a set $S$. For $b\in\ensuremath{\mathbf{R}}^p$, its $\ell_q$ norm is denoted as $|b|_q$ = $(\sum_{j\in[p]}|b_j|^q)^{1/q}$ if $q\in[1,\infty)$ and $|b|_\infty = \max_{j\in[p]}|b_j|$ if $q=\infty$. For a group structure $\mathcal{G}$, the $\ell_{2,1}$ group norm of $b\in\ensuremath{\mathbf{R}}^p$ is defined as $\|b\|_{2,1}=\sum_{G\in\mathcal{G}}|b_G|_2$. For $\mathbf{u},\mathbf{v}\in\ensuremath{\mathbf{R}}^J$, the empirical inner product is defined as $\langle \mathbf{u},\mathbf{v}\rangle_J = J^{-1}\sum_{j=1}^J u_jv_j$ with the induced empirical norm $\|.\|_J^2=\langle.,.\rangle_J=|.|_2^2/J$. For a symmetric $p\times p$ matrix $A$, let $\mathrm{vech}(A)\in\ensuremath{\mathbf{R}}^{p(p+1)/2}$ be its vectorization consisting of the lower triangular and the diagonal elements. Let $A_G$ be a sub-matrix consisting of rows of $A$ corresponding to indices in $G\subset[p]$. If $G=\{j\}$ for some $j\in[p]$, then we simply write $A_G=A_j$. Let $\|A\|_\infty=\max_{j\in[p]}|A_j|$ be the matrix norm. For $a,b\in\ensuremath{\mathbf{R}}$, we put $a\vee b = \max\{a,b\}$ and $a\wedge b = \min\{a,b\}$. Lastly, we write $a_n\lesssim b_n$ if there exists a (sufficiently large) absolute constant $C$ such that $a_n\leq C b_n$ for all $n\geq 1$ and $a_n\sim b_n$ if $a_n\lesssim b_n$ and $b_n\lesssim a_n$.
Motivated by our empirical application, we allow the high-dimensional set of predictors to be sampled at a higher frequency than the target variable. Let $K$ be the total number of time-varying predictors $\{x_{i,t-(j-1)/m,k}:i\in[N],t\in[T],j\in[m],k\in[K]\}$ possibly measured at some higher frequency with $m$ observations for every low-frequency period $t\in[T]$ and every entity $i\in[N]$. Consider the following (mixed frequency) panel data regression
where $h\geq 0$ is the prediction horizon, $\alpha_i$ is the entity-specific intercept, and
is a high-frequency lag polynomial with $\beta_k=(\beta_{1,k},\dots,\beta_{m,k})^\top\in\ensuremath{\mathbf{R}}^m$. More generally, the frequency can also be specific to the predictor $k\in[K]$, in which case we would have $m_k$ instead of $m$. We can also absorb the (low-frequency) lags of $y_{i,t}$ in covariates. When $m$ = 1, we retain the standard panel data regression model
while $m>1$ signifies that the high-frequency lags of $x_{i,t,k}$ are also included. The large number of predictors $K$ with potentially large number of high-frequency measurements $m$ can be a rich source of predictive information, yet at the same time, estimating $N + m\times K$ parameters is costly and may reduce the predictive performance in small samples.
To reduce the proliferation of lag parameters, we follow the MIDAS literature; see ghysels2006predicting, ghysels:mrf, and babii2020inference,babii2020machine. Instead of estimating $m$ individual slopes of high-frequency covariate $k\in[K]$ in equation ((ref)), with some abuse of notation, we estimate a weight function $\omega$ parameterized by $\beta_k\in\ensuremath{\mathbf{R}}^L$ with $L< m$
where
and $(w_l)_{l\geq 0}$ is a collection of $L$ approximating functions, called the dictionary. An example of a dictionary is the set of orthogonal Legendre polynomials on $[0,1]$ that can be computed via the Rodrigues' formula $w_l(s)=\frac{1}{l!}\frac{\ensuremath{\mathrm{d}}^l}{\ensuremath{\mathrm{d}} s^l}(s^2-s)^l.$\footnote{The Legendre polynomials have the universal approximation property and can approximate any continuous function uniformly on $[0,1]$. At the same time they can generate a rich family of MIDAS weights with a relatively small number of parameters which is attractive in time series applications where the signal-to-noise ratio is often low.} For instance, the first five elements are
More generally, we can use Gegenbauer polynomials, trigonometric polynomials, or wavelets. The orthogonal polynomials usually have better numerical properties than their popular non-orthogonal counterpart, such as the almon1965distributed lag structure. The attractive feature of linear in parameters dictionaries is that we can map the MIDAS regression to the linear regression framework that can be solved via a convex optimization. To that end, define $\mathbf{x}_i = (X_{i,1}W,\dots,X_{i,K}W)$, where for each $k\in[K]$, $X_{i,k} = (x_{i,t-(j-1)/m,k})_{t\in[T],j\in[m]}$ is a $T\times m$ matrix of predictors and $W=(w_l((j-1)/m)/m)_{j\in[m],0\leq l\leq L-1}$ is an $m\times L$ matrix corresponding to the dictionary $(w_l)_{l\geq 0}$. In addition, let $\mathbf{y}_i$ = $(y_{i,1+h},\dots,y_{i,T+h})^\top$ and $\mathbf{u}_i$ = $(u_{i,1},\dots,u_{i,T})^\top.$ Then the regression equation after stacking time series observations for each $i\in[N]$ is
where $\iota\in\ensuremath{\mathbf{R}}^T$ is the all-ones vector and $\beta\in\ensuremath{\mathbf{R}}^{LK}$ is a vector of slopes. Lastly, put $\mathbf{y} = (\mathbf{y}_1^\top,\dots, \mathbf{y}_N^\top)^\top$, $\mathbf{X}=(\mathbf{x}_1^\top, \dots, \mathbf{x}_N^\top)^\top$, and $\mathbf{u} = (\mathbf{u}_1^\top,\dots,\mathbf{u}_N^\top)^\top$. Then the regression equation after stacking all cross-sectional observations is
where $B=I_N\otimes\iota$, $\alpha=(\alpha_1,\dots,\alpha_N)$, and $\otimes$ is the Kronecker product.
The MIDAS approach allows us to effectively reduce the dimensionality pertaining to the high-frequency lags. Alternatively, we may apply what is known as the UMIDAS scheme, see e.g., foroni2015unrestricted, and directly estimate the coefficients associated with each high-frequency covariate lags separately (see equation ((ref)) in Section (ref) for example). Such a strategy, which as foroni2015unrestricted argue works in single regressions when the ratio high to low-frequency sampling is small, may not be appealing in high-dimensional cases, as the estimation and prediction performance deteriorates due to the potentially large number of coefficients; see babii2020machine for further discussion. Also, while assuming that the individual lag coefficients in equation ((ref)) are approximately sparse is highly restrictive, the approximate sparsity of slopes of the dictionary elements $(w_l)_{l\geq 0}$ is plausible. For instance, if $w_0(s)=1$ with $\beta_{0,k}\ne 0$ and $\beta_{l,k}=0,\forall l\geq 1$, we recover the averaging of high-frequency lags of covariate $k$ as a special case. More generally, the weight $\omega$ may be a decreasing function over lags and we may want to learn its shape from the data maximizing the predictive performance.\footnote{See ball2013dissecting and ball2018mixed for further discussion on interpreting the shape of MIDAS polynomials in accounting data applications considered in our empirical application.}
Given that the number of potential predictors $K$ can be large, additional regularization can improve the predictive performance in small samples. To that end, we take advantage of the sg-LASSO regularization that was shown to be attractive for individual time series ML regressions in babii2020machine. The fixed effects panel data estimator with sparse-group regularization solves
where $\|.\|_{NT}^2 = |.|^2/(NT)$ is the empirical norm and
is a regularizing functional, which is a linear combination of LASSO and group LASSO penalties. The parameter $\gamma\in[0,1]$ determines the relative weights of the $\ell_1$ (sparsity) and the $\ell_{2,1}$ (group sparsity) norms, while the amount of regularization is controlled by the regularization parameter $\lambda\geq 0$. Recall also that for a group structure $\mathcal{G}$ described as a partition of $[p]=\{1,2,\dots,p\}$, the group LASSO norm is computed as $\|b\|_{2,1}=\sum_{G\in\mathcal{G}}|b_G|_2$. The group structure is assumed to be known to the econometrician, which in our setting corresponds to time series lags of covariates. More generally, we may also combine covariates of a similar nature in groups. Throughout the paper we assume that groups have fixed size, which is well-justified in our empirical applications.\footnote{See babii2020high for a continuous-time mixed-frequency regression where the group size is allowed to increase with the sample size under the in-fill asymptotics.} Therefore, the selection of covariates is performed by the group LASSO penalty, which encourages sparsity between groups. In addition, the $\ell_1$ LASSO norm promotes sparsity within groups and allows us to learn the shape of the MIDAS weights from the data.
It is worth mentioning that the linear in parameters approximation to the MIDAS weight function leads to the convex optimization parameter problem in equation ((ref)) that can be solved efficiently, e.g., via the proximal gradient descent algorithm, or its block-coordinate descent versions. In contrast, a popular beta weights leads to a nonlinear non-convex optimization problem that becomes challenging to solve in high-dimensions; cf.\ marsilli2014variable and khalaf2020dynamic.
In this section, we provide the theoretical analysis of predictive performance of regularized panel data regressions with the sg-LASSO regularization, including the standard LASSO and the group LASSO regularizations as special cases. It is worth stressing that the analysis of this section is not tied to the mixed-frequency data setting and applies to the generic high-dimensional panel data regularized with the sg-LASSO penalty function. Importantly, we focus on panels consisting of potentially persistent $\tau$-mixing time series with polynomial tails. Consider a generic panel data projection with a countable number of predictors
This model subsumes the mixed-frequency data regressions as a special case, in which case covariates are obtained, e.g., from the aggregation with Legendre polynomials. The covariates may also include the time-varying covariates common for all entities (macroeconomic factors), lags of $y_{i,t}$, the intercept, as well as additional lags of a baseline covariate.
We measure the persistence of the data with $\tau$-mixing coefficients. For a $\sigma$-algebra $\mathcal{M}$ and a random vector $\xi\in\ensuremath{\mathbf{R}}^l$, put
where $\mathrm{Lip}_1=\{f:\ensuremath{\mathbf{R}}^l\to\ensuremath{\mathbf{R}}:\; |f(x) - f(y)| \leq |x-y|_1 \}$ is a set of $1$-Lipschitz functions from $\ensuremath{\mathbf{R}}^l$ to $\ensuremath{\mathbf{R}}$.\footnote{See dedecker2004coupling and dedecker2005new for equivalent definitions.} For a stochastic process $(\xi_t)_{t\in\ensuremath{\mathbf{Z}}}$ with a natural filtration generated by its past $\mathcal{M}_t=\sigma(\xi_t,\xi_{t-1},\dots)$, the $\tau$-mixing coefficients are defined as
where the supremum is taken over all $t,t_1,\dots,t_j\in\ensuremath{\mathbf{Z}}$. If $\tau_k\downarrow0$, as $k\uparrow\infty$ then the process is called $\tau$-mixing. The class of $\tau$-mixing processes can be placed somewhere between the $\alpha$-mixing processes and mixingales --- the $\tau$-mixing condition is less restrictive than the $\alpha$-mixing condition,\footnote{The class of $\alpha$-mixing processes is too restrictive for the predictive linear projection model with covariates and autoregressive lags; see also babii2020machine, Proposition A.3.1.} yet at the same time, there exists a convenient for us coupling result for $\tau$-mixing processes, which is not the case for the mixingales or near-epoch dependent processes; see dedecker2003new and dedecker2004coupling,dedecker2005new for more details. This allows us to obtain concentration inequalities and performance guarantees for the sg-LASSO estimator; see Appendix (ref) for more details.
For pooled regressions, we assume that all entities share the same intercept parameter $\alpha_1=\dots=\alpha_N=\alpha$. The pooled sg-LASSO estimator $\hat\rho=(\hat\alpha,\hat\beta^\top)^\top$ solves
Define (a) $z_{i,t}=(1,x_{i,t}^\top)^\top$, where $x_{i,t}\in\ensuremath{\mathbf{R}}^p$ is a vector of predictors, (b) $u_i=(u_{i,1},\dots,u_{i,T})$ and (c) $x_i = (x_{i,1}^\top,\dots,x_{i,T}^\top)^\top$ for $i\in[N]$. The following assumption imposes mild restrictions on the data.
Note that we do not impose stationarity over $t\in\ensuremath{\mathbf{Z}}$ and require that only $2+\epsilon$ moments exist with $\epsilon>0$, which is a realistic assumption in our empirical application and more generally for datasets encountered in time series and financial econometrics applications. Note also that the time series dependence is assumed to fade away relatively slowly --- at a polynomial rate as measured by the $\tau$-mixing coefficients.
Next, we assume that the $(1+p)\times(1+p)$ matrix $\Sigma_{N,T}=\frac{1}{NT}\sum_{i=1}^N\sum_{t=1}^T\ensuremath{\mathbb{E}}[z_{i,t}z_{i,t}^\top]$ exists and is non-singular uniformly over $N,T,p$:
Assumption (ref) is satisfied for the spiked identity and Topelitz covariance structures. It can be interpreted as a completeness condition, see babii2017completeness, and can also be relaxed to the restricted eigenvalue condition imposed on the population covariance matrix $\Sigma_{N,T}$; see babii2020machine. We can also allow for $\gamma_{\min}\downarrow0$ as $N,T,p\uparrow\infty$, in which case $\gamma_{\min}^{-1}$ would slow down the convergence rates in oracle inequalities and could be interpreted as a measure of ill-posedness; see also carrasco2007linear.
Lastly, we assume that the regularization parameter $\lambda$ scales appropriately with the number of covariates $p$, the length of the panel $T$, the size of the cross-section $N$, and a certain exponent $\kappa$ that depends on the tail parameter $q$ and the persistence parameter $a$. The precise order of the regularization parameter is described by the Fuk-Nagaev inequality for long panels appearing in the Appendix; see Theorem (ref).
Our first result is the oracle inequality for the pooled sg-LASSO estimator described in equation ((ref)). The result allows for misspecified regressions with a non-trivial approximation error in the sense that we consider more generally
where $\mathbf{m}\in\ensuremath{\mathbf{R}}^{NT}$ is approximated with $\mathbf{Z}\rho$, $\mathbf{Z}=(\iota,\mathbf{X})$, $\iota\in\ensuremath{\mathbf{R}}^{NT}$ is all-ones vector, and $\rho = (\alpha,\beta^\top)^\top$. The approximation error $\mathbf{m}-\mathbf{Z}\rho$ might come from the fact that the MIDAS weight function may not have the exact expansion in terms of the specified dictionary or from the fact that some of the relevant predictors are not included in the regression equation. To state the result, let $S_0=\{j\in[p]:\; \beta_j\ne 0 \}$ be the support of $\beta$ and let $\mathcal{G}_0=\{G\in\mathcal{G}:\;\beta_G\ne 0 \}$ be the group support of $\beta$. Consider the effective sparsity of the sparse-group structure, defined as $s^{1/2} = \gamma\sqrt{{|S_0|}} + (1-\gamma)\sqrt{{|\mathcal{G}_0|}}$. Note that $s$ is proportional to the sparsity $|S_0|$, when $\gamma=1$ and to the group sparsity $|\mathcal{G}_0|$ when $\gamma = 0$. Define $r_{N,T}^{\rm pooled} = s^{\tilde\kappa}p^2/(NT)^{\tilde \kappa - 1} + p^2\exp(-cNT/s^2)$.
The proof of this result can be found in the Appendix. Theorem (ref) describes the non-asymptotic oracle inequalities for the prediction and the estimation accuracy in the environment where the number of regressors $p$ is allowed to scale with the effective sample size $NT$. Importantly, the result is stated under the weak tail and persistence conditions in Assumption (ref). Parameters $\kappa$ and $\tilde\kappa$ are the dependence-tails exponents for stochastic processes driving the regression score and the covariance matrix respectively. Theorem (ref) shows that the prediction and the estimation accuracy of pooled panel data regressions improves when the sparse-group structure is taken into account. Indeed, for the LASSO regression, the effective sparsity reduces to $s^{1/2}=\sqrt{|S_0|}$, which is larger than $\gamma\sqrt{|S_0|} + (1-\gamma)\sqrt{|\mathcal{G}_0|}$ in the case of sg-LASSO.
Next, we consider the convergence rates of the prediction and estimation errors. The following assumption considers a simplified setting, where the approximation error vanishes sufficiently fast, and the total number of regressors vanishes sufficiently fast with the effective sample size $NT$.
Note that Assumption (ref) allows for (1) $N\to\infty$ while $T$ is fixed; (2) $T\to\infty$ while $N$ is fixed; and (3) both $N\to\infty$ and $T\to\infty$ without restricting the relative growth of the two. The following result describes the prediction and the estimation convergence rates in the asymptotic environment outlined in Assumption (ref) and is an immediate consequence of Theorem (ref).
Corollary (ref) describes the prediction and the estimation accuracy of pooled sparse-group panel data regressions. It suggests that the predictive performance of the sg-LASSO (and consequently LASSO and group LASSO) regressions may deteriorate when regression errors and/or predictors are heavy-tailed or when the data are extremely persistent. However, for geometrically ergodic Markov processes, e.g., stationary AR(1) process, the $\tau$-mixing coefficients decline geometrically fast, so that $\kappa\approx q$ and $\tilde\kappa\approx \tilde q$. In this case, the prediction accuracy scales approximately at the rate $O_P\left(\frac{p^{2/q}}{(NT)^{2-2/q}}\vee \frac{\log p}{NT}\right)$ and the predictive performance may be affected only by the tails constant $q$.
If additionally the data are sub-Gaussian, then moments of all order $q\geq 2$ exist, and for any particular effective sample size $NT$, the first term can be made arbitrarily small relatively to the second term. In this case we recover the $O_P\left(\frac{\log p}{NT}\right)$ rate typically obtained for sub-Gaussian data. On the other hand, if the polynomial tail dominates, then we need $p = o((NT)^{q -1})$ for the prediction and the estimation consistency provided that $\tilde q\geq 2q-1$ and the sparsity constant $s$ is fixed. In this case, we have a significantly weaker requirement than the $p=o(T^{q - 1})$ condition needed for time series regressions in babii2020machine. Moreover, since $q>2$, $p=o((NT)^{q-1})$ can be significantly weaker than the $p=o(NT)$ condition typically needed for QMLE/GMM estimators without regularization.
Theorem (ref) and Corollary (ref) imply two practical consequences: (1) one may want to exclude (or suitably transform) the heavy-tailed series from the high-dimensional predictive regressions based on the preliminary estimates of the tail index, e.g., using the Hill estimator; (2) if the individual heterogeneity can be ignored, then pooling panel data can improve significantly the predictive performance. In the latter case, one can also preliminary cluster similar series in groups, e.g., based on the unsupervised clustering algorithms, which may strike a good balance between the pooling benefits and heterogeneity.
Pooled regressions are attractive since the effective sample size $NT$ can be huge, yet the heterogeneity of individual time series may be lost. If the underlying series have a substantial heterogeneity over $i\in[N]$, then taking this into account might reduce the projection error and improve the predictive accuracy. At a very extreme side, the cross-sectional structure can be completely ignored and individual time series regressions can be used for prediction. The fixed effects panel data regressions strike a good balance between the two extremes controlling for heterogeneity with entity-specific intercepts.
The fixed effects sg-LASSO estimator $\hat\rho=(\hat\alpha^\top,\hat\beta^\top)^\top$ solves
where $B=I_N\otimes\iota$, $I_N$ is $N\times N$ identity matrix, $\iota\in\ensuremath{\mathbf{R}}^T$ is an all-ones vector, and $\Omega$ is the sg-LASSO regularizing functional. It is worth stressing that the design matrix $\mathbf{X}$ does not include the intercept and that we do not penalize the fixed effects, that are typically not sparse. By Fermat's rule, the first-order conditions are
for some $z^*\in\partial\Omega(\hat\beta)$, where $b\mapsto \partial\Omega(b)$ is the subdifferential of $\Omega$ and $M_B = I - B(B^\top B)^{-1}B^\top$ is the orthogonal projection matrix. It is easy to see from the first-order conditions that the estimator of $\hat\beta$ is equivalent to 1) penalized GLS estimator for the first-differenced regression; 2) penalized OLS estimator for the regression written in the deviation from time means; and 3) penalized OLS estimator where the fixed effects are partialled-out. Therefore, the equivalence between the three approaches is not affected by the penalization; cf.\ arellano2003panel for low-dimensional panels.
With some abuse of notation, redefine
We will assume that the smallest eigenvalue of $\Sigma_{N,T}$ is uniformly bounded away from zero by some constant. Note that if $x_{i,t}\sim N(0,I_p)$, then $\Sigma_{N,T}=I_{N+p}$ and this assumption is trivially satisfied.
The order of the regularization parameter is governed by the Fuk-Nagaev inequality for long panels; see Appendix, Theorem (ref).
Similarly to the pooled regressions, we state the oracle inequality allowing for the approximation error. For fixed effects regressions, with some abuse of notation we redefine $\mathbf{Z} = (B,\mathbf{X})$ and $\rho = (\alpha^\top,\beta^\top)^\top$. Put also $r_{N,T}^{\rm fe}=p(s\vee N)^{\tilde \kappa} T^{1-\tilde \kappa}(N^{1-\tilde \kappa/2} + pN^{1-\tilde\kappa}) + p(p\vee N)e^{-cNT/(s\vee N)^2}$ with $\tilde\kappa = ((\tilde a + 1)\tilde q - 1) / (\tilde a + \tilde q - 1)$ and some $c>0$.
Theorem (ref) states a non-asymptotic oracle inequality for the prediction error in the fixed effects panel data regressions estimated with the sg-LASSO. To see clearly, how the prediction accuracy scales with the sample size, we make the following assumption.
The following corollary is an immediate consequence of Theorem (ref).
Corollary (ref) allows for $s,p,N,T\to\infty$ at appropriate rates. However, we pay an additional price for estimating $N$ fixed effects which plays a similar role to the effective dimension of covariates. An immediate practical implication is that to achieve accurate predictions with high-dimensional fixed effect regressions, the panel has to be sufficiently long to offset the estimation error of the individual fixed effects. Likewise, the tails and the persistence of the data may also reduce the prediction accuracy in small samples through $\kappa$, which is approximately equal to $q$ for geometrically decaying $\tau$-mixing coefficients.
In this section, we develop the debiased inferential methods for pooled panel data regressions. For a vector $\rho\in\ensuremath{\mathbf{R}}^{p+1}$, we use $\rho_G\in\ensuremath{\mathbf{R}}^{|G|}$ to denote the subvector of elements of $\rho\in\ensuremath{\mathbf{R}}^{p+1}$ indexed by $G\subset[p+1]$. Let $B = \hat\Theta\mathbf{Z}^\top(\mathbf{y} - \mathbf{Z}\hat\rho) / NT$ denote the bias-correction for the sg-LASSO estimator, where $\hat\Theta$ is the nodewise LASSO estimator of the precision matrix $\Theta=\Sigma^{-1}$, where $\Sigma=\ensuremath{\mathbb{E}}[z_{i,t}z_{i,t}^\top]$. For pooled panel data, this estimator can be obtained as follows:
Let $v_{i,t,j}=z_{i,t,j}-\sum_{k\ne j}\mu_{j,k}z_{i,t,k}$ be the regression error for $j^{\rm th}$ nodewise LASSO regression. Let $s_j$ be the number of non-zero elements in $j^{\rm th}$ row of precision matrix $\Theta_j$, and put $S=\max_{j\in G}s_j$, and $s^*=s\vee S$.
The following assumption describes an additional set of conditions for the debiased central limit theorem.
Assumption (ref) (i) requires that the conditional variance of the regression error is bounded. Condition (ii) requires that the rows of the precision matrix have bounded $\ell_1$ norm and is a plausible assumption in the high-dimensional setting, where the inverse covariance matrix is often sparse. Condition (iii) is a mild restriction needed for the consistency of the sample variance of regression errors. The rate conditions in (iv) are similar to the condition used in babii2020inference. Lastly, condition (v) is trivially satisfied when the projection coefficients are sparse and, more generally, it requires that the misspecification error vanishes asymptotically sufficiently fast.
The following result describes a large-sample approximation to the distribution of the debiased sg-LASSO estimator with serially correlated heavy-tailed errors.
Theorem (ref) applies to panel data consisting of non-Gaussian, heavy-tailed, and persistent time series under the large $N$ and $T$ large sample approximation. In contrast to the fixed $T$ approximations, Theorem (ref) leads to more precise inference, e.g., the standard errors and the length of confidence intervals would scale at $O(1/\sqrt{NT})$ rate instead of $O(1/\sqrt{N})$ that we typically encounter for fixed $T$ approximations.
To estimate $\Xi_G$, we can use the following pooled HAC estimator
where $\hat{\Gamma}_{k,i} = \hat\Theta_G\left(\frac{1}{T}\sum_{t=1}^{T-k}\hat u_{i,t}\hat u_{i,t+k} x_{i,t}x_{i,t+k}^\top\right)\hat\Theta_G^\top$, $\hat u_{i,t}$ is the sg-LASSO residual, and $\hat\Gamma_{-k,i}=\hat{\Gamma}_{k,i}^\top$. The kernel function $K:\ensuremath{\mathbf{R}}\to[-1,1]$ with $K(0)=1$ is puts less weight on more distant noisy covariances, while $M_T\uparrow\infty$ is a bandwidth (or lag truncation) parameter; see babii2020inference for more details as well as formal results on the validity of HAC-based inference using sg-LASSO residuals.
In this section, we assess the finite sample performance of the Granger causality tests for high-dimensional pooled panel data MIDAS regressions. A first subsection describes the design, followed by a second reporting the findings.
We simulate the data from the following DGP:
where $i\in[N]$, $t\in[T]$, $\alpha$ is the the common intercept, $\frac{1}{m}\sum_{j=1}^m \omega((j-1)/m;\beta_k)$ is the weight function for $k$-th high-frequency covariate and the error term is $u_{i,t} \sim_{i.i.d.}N(0,4)$. The DGP corresponds to the target variable of interest $y_{i,t}$ driven by one autoregressive lag augmented with high-frequency series. The DGP is therefore a pooled MIDAS panel data model.
We set $\rho=0.15$ and take the first high-frequency regressor, $k$ = 1, as relevant, i.e.\ the first regressor Granger causes the response variable. We are interested in quarterly/monthly data, and use four quarters of data for the high-frequency regressors so that $m$ = 12. The high-frequency regressors are generated as $K$ i.i.d.\ realizations of univariate autoregressive (AR) processes $x_h = \rho x_{h-1}+ \varepsilon_h,$ where $\rho=0.7$ and $\varepsilon_h\sim_{i.i.d.}N(0,1)$, where $h$ denotes the high-frequency sampling. For the DGP we rely on a commonly used weighting scheme in the MIDAS literature, namely the weights $\omega(s;\beta_k)$ for the only relevant high-frequency regressor $k=1$ determined by the beta density, $\mathrm{Beta}(3,3)$; see ghysels2007midas or ghysels2019estimating, for further details. The empirical estimation involves MIDAS regressions with Legendre polynomials of degree $L=3$. Lastly, we draw the intercepts $\alpha \sim \text{Uniform}(-4,4).$ Throughout the experiment, we fix the sample sizes to $T$ = 50 and $N$ = 30.
We compare the empirical size and power of the Granger causality test under different structures placed on the regression models.
First, we compare sg-LASSO-MIDAS with LASSO-UMIDAS pooled panel data models. The former exploits the group structure of covariates by applying the sg-LASSO penalty function and a flexible way to model lags for each covariate using the MIDAS weight functions parametrized by low-dimensional coefficients. The latter pertains to the unstructured LASSO estimator together with the UMIDAS scheme. Introduced by foroni2015unrestricted, UMIDAS consists of estimating a regression coefficient for each high-frequency lag separately, and therefore the weight function for each covariate is
where $b_{j,k}$ is a regression coefficient associated with each high-frequency lag. We estimate regression coefficients by applying the standard unstructured LASSO estimator; hence we call the model LASSO-UMIDAS.
Second, we compare the pooled panel with individual time series regressions, for sg-LASSO-MIDAS and LASSO-UMIDAS, where the former exploits the benefits of the panel structure and the latter does not. In this case, we take the first sample $i=1$ to compute empirical size and power of the Granger test for the individual regression models. babii2020inference propose tests of Granger causality in univariate regularized regressions and high-dimensional data.
In Table (ref), we report the empirical rejection frequency (ERF) for the Granger causality test based on the HAC estimator with two different kernel functions, Parzen and Quadratic spectral, and two different estimation strategies, sg-LASSO-MIDAS and LASSO-UMIDAS. We test whether the first high-frequency covariate Granger causes the low-frequency series, which corresponds to the DGP potential causal pattern. We report results for a set of bandwidth parameters, denoted $M_T$ = 10, 20 and 30. The reported results are based on 2000 Monte Carlo replications.
To assess the performance we scale the Beta density function by multiplying it with a constant $a\in\{0,1/5,1/4,1/3\}$, i.e.\ the weight function for the relevant covariate is:
For $a=0$, the ERF shows the empirical size of the test for the nominal level of 5%, while $a\in\{1/5,1/4,1/3\}$ the ERF shows the empirical power of the Granger causality test. For the larger scaling constant $a,$ the alternatives are separated further away from the null hypothesis and the Granger causality test is expected to perform better.
The results reported in Table (ref) show that the Granger causality test based on the sg-LASSO-MIDAS has empirical size close to the nominal level of $5\%$. In contrast, the LASSO-UMIDAS leads to undersized Granger causality tests with size distortions around 0.01. The Granger causality test based on the sg-LASSO-MIDAS has also better empirical power against each of the alternative hypotheses $a\in\{1/5,1/4,1/3\}$. Additionally, it approaches 1 much faster as opposed to the LASSO-UMIDAS.
The results for individual regressions reveal worse performance compared to pooled panel data regressions, hence showing the usefulness of pooling the data. The empirical size shows considerable size distortions of around 0.05. Tests for individual regressions have worse power compared to the pooled panel data cases. Nonetheless, similar to the pooled panel data cases, the sg-LASSO-MIDAS estimation method seems to have better empirical power when comparing to LASSO-UMIDAS.
Overall, the results of the Monte Carlo experiments indicate that the structured regularization leads to better Granger causality tests in small samples and that pooling individual series improves the results even further.
In this section we revisit a topic raised by ball2018automated and carabias2018real. Their empirical findings suggest that analysts tend to focus on their firm/industry when making earnings predictions while not fully taking into account the impact of macroeconomic events. While their findings were suggestive, there was no formal testing in a data-rich environment. The theory established in the previous sections allows us to do so.
More specifically, we consider the earnings of 210 US firms using a set of predictors sampled at mixed frequencies --- quarterly, monthly and daily series. We use 26 predictors (and their lags), including traditional macro and financial series as well as non-standard series generated by textual analysis of financial news.
The full sample consists of observations between the \(1^{st}\) of January, 2000 and the \(30^{th}\) of June, 2017. Due to the lagged dependent variables in the models, our effective sample starts at the third fiscal quarter of 2000. We collected data from CRSP and I/B/E/S to compute quarterly earnings and firm-specific financial covariates; RavenPack was used to compute daily firm-level textual-analysis-based data; real-time monthly macroeconomic series are from the ALFRED; FRED is used to compute daily financial markets data and, lastly, monthly news attention series extracted from the {\it Wall Street Journal} articles were retrieved from bybee2019structure.\footnote{The dataset is publicly available at \href{http://www.structureofnews.com/}{http://www.structureofnews.com/}.} Table (ref) provides a list of the variables used in our analysis, whereas Online Appendix Section (ref) covers a detailed description of the RavenPack data. Finally, the list of all firms we consider in our analysis appears in Online Appendix Table (ref). Table (ref) has six panels, namely three panels of firm-level series: A1 -- describes earnings data, B1 -- describes daily firm-level stock market data, and C1 -- describes daily firm-level sentiment data series. The remaining three panels are: A2 -- describes real-time monthly macro series, B2 -- describes daily financial markets data, and C2 -- describes monthly news attention series. In the models we include 365 daily lags, 12 monthly lags and 4 quarterly lags respectively.
Whether analysts leave money on the table amounts to testing whether forecast errors in earnings can be predicted by current information variables. Hence, this amounts to performing something akin to the Granger causality test. In our empirical application we are dealing with a panel, and it is important to exploit the multivariate data structure to perform such tests.
We analyze the difference between realized earnings and analysts' predictions, i.e., the response variable $y_{i,t+1}$ is computed by taking the difference between realized earnings, denoted $e_{i,t+1}$, and the median of analysts' predictions for the quarter $t+1$, denoted $f_{i,t+1|t}$,
We then fit the following pooled panel data MIDAS model using sg-LASSO estimator:
We test which factors Granger cause future errors of earnings forecasts made by the analysts. In the sg-LASSO, groups are defined as all lags of a single covariate $k;$ Legendre polynomials up to degree three are applied to all weight functions $\psi(L^{1/m};\beta_k).$ We use 10-fold cross-validation to tune both $\lambda$ and $\gamma$, where we define folds as adjacent blocks over the time series dimension to take into account the time series dependence. Similarly, we estimate the precision matrix using nodewise LASSO regressions selecting the tuning parameter in a similar vein. The results are reported in Table (ref).
In Panel (A) of Table (ref) we find that the AR(1) lag is significant, leading us to conclude that the prediction errors made by the analysts are persistent. The autoregressive coefficient is significant throughout all specifications of the models, including in a simple pooled AR(1) model. In the latter case, the AR(1) coefficient is estimated to be 0.147.
Panel (B) of Table (ref) reports that beyond the AR(1) we find that the highly significant covariates are TED rate, CPI inflation and real GDP growth. These results support previous findings that analysts tend to miss information associated with macroeconomic conditions --- including real GDP growth and the TED spread, which is an indicator of measure credit risk. The latter is rather surprising, as it indicates that analysts tend to miss out on credit risk information at the macro level in their earnings forecasts. Lastly, the term spread (10-year less 3-month treasury yield), often viewed as a business cycle indicator, is also significant at the 10% level.
Finally, in Panel (C) of Table (ref) we report results based on the unstructured LASSO applying UMIDAS for the lag polynomials of each covariate. The findings reveal similar results for the TED rate, but notably miss real GDP and CPI inflation as significant covariates.
In Table (ref) we show results based on a different way of pooling analysts' prediction errors $y_{i,t+1}$. We split the data into two parts based on how large the average disagreement among analysts is. For each firm, we compute the forecast disagreement as the difference between 95% and 5% percentile of the empirical forecast distribution and take the average over the sample. We sort from high to low disagreement and split the sample of firms into two subsamples of equal size. The results show that macro variables which are significant for the full sample are also significant for the large disagreement subsample. On the other hand, little significance is reported for the low disagreement subsample. In this case, only the AR(1) lag and stock returns are significant at the 5% significance level.
Lastly, in Figure (ref) we plot the ratio of firms for which we find Granger causality based on individual regressions versus panel models. In Panel (a) we plot the ratios for sg-LASSO estimator using MIDAS weighting scheme while in Panel (b) we plot the ratios for the LASSO estimator with UMIDAS scheme. The plot shows ratios for each covariate representing the fraction with respect to sg-LASSO (Panel (a)) or LASSO with UMIDAS (Panel (b)) each covariate is significant by running individual regressions. For example, the AR(1) lag is significant for around 30% (0.3) of firms when running individual sg-LASSO-MIDAS regressions. Some covariates that are not significant in pooled panels are significant for some firms; therefore, we show results for all covariates, including those that are not significant in pooled panel cases. We also show how the ratios differ for low (dark-gray color) versus high disagreement (light-gray color) firms. They represent whether a specific firm we run an individual regression for is in the high-disagreement versus low-disagreement subsample. Interestingly, the largest ratios are for AR(1), TED rate, Real GDP, CPI inflation and term spread in the case of sg-LASSO-MIDAS. Moreover, the portion of firms in the high disagreement subsample seem to have the largest ratios. In the case of LASSO-UMIDAS, the ratios show a less clear pattern, with only the AR(1) and TED rate covariates significant for a larger number of firms.
This paper introduced a new class of high-dimensional panel data regression models with dictionaries and sg-LASSO regularization. This type of regularization is an especially attractive choice for predictive panel data regressions, where the low- and/or the high-frequency lags define a clear group structure. The estimator nests the LASSO and the group LASSO estimators as special cases. Our theoretical treatment allows for heavy-tailed data frequently encountered in financial time series. To that end, we obtain a new panel data concentration inequality of the Fuk-Nagaev type for $\tau$-mixing processes, which allows us to establish oracle inequalities that are used subsequently to develop the debiased HAC inference for the panel data sg-LASSO estimator.
Using the theory of HAC-based inference for pooled panel data regressions developed in our paper, our empirical analysis revisits a topic raised by earlier literature that analysts tend to focus on firm and/or industry information when forming earnings forecasts, while not fully taking into account the macroeconomic data. Our results suggest that indeed analysts tend to miss on macro information, i.e., macro variables turn out to be significant in pooled panel regression models.