High-Dimensional Granger Causality Tests with an Application to VIX and News

High-Dimensional Granger Causality Tests with an Application to VIX and News

\def\spacingset#1{ {#1}} \spacingset{1}

{\it Keywords:} Granger causality, high-dimensional time series, heavy tails, inference, tau-mixing, HAC estimator, sparse-group LASSO, Fuk-Nagaev inequality. \thispagestyle{empty} \setcounter{page}{0} \spacingset{1.5}

Introduction

Modern time series analysis is increasingly using high-dimensional datasets, typically available at different frequencies. Conventional time series are often supplemented with non-traditional data, such as the high-dimensional data coming from the natural language processing. For instance, bybee2019structure extract 180 topic attention series from the over 800,000 daily Wall Street Journal news articles during 1984-2017 that have shown by babii2020midasml to be a useful supplement to more traditional macroeconomic and financial datasets for nowcasting US GDP growth.

In his seminal paper, Clive Granger defined causality in terms of high-dimensional time series data. His formal definition, see granger1969investigating, considered all the information accumulated in the universe up to time $t - 1$ (a process he called $U_t$) and examined predictability using $U_t$ with and without a specific series of interest $Y_t.$ It is still an open question how to implement Granger's test in a high-dimensional time series setting. It is the purpose of this paper to do this via regularized regressions using HAC-based inference. In a sense, we are trying to implement Granger's original idea of causality.\footnote{There exists an extensive literature on causal inference with machine learning methods within the static Neyman-Rubin's potential outcomes framework; see athey2019machine for the excellent review and further references.}

Following babii2020midasml, we focus on the structured sparsity approach based on the sparse-group LASSO (sg-LASSO) regularization for the high-dimensional time series analysis. The sg-LASSO allows capturing the group structures present in high-dimensional time series regressions where a single covariate with its lags constitutes a group. Alternatively, we can also combine covariates of similar nature in groups. An attractive feature of this estimator is that it encompasses the LASSO and the group LASSO as special cases, hence, it allows improving upon the unstructured LASSO in the high-dimensional time-series setting. At the same time, the sg-LASSO can learn the distribution of time series lags in a data-driven way solving elegantly the model selection problem that dates back to fisher1937note.\footnote{The distributed lag literature can be traced back to fisher1925our; see also almon1965distributed, sims1971discrete, and shiller1973distributed, as well as more recent mixed frequency data sampling (MIDAS) approach in ghysels2006predicting, ghysels2007midas, and andreou2013should.} In particular, the group structure can also accommodate data sampled at different frequencies as discussed in detail by babii2020midasml.

The proper inference for time-series data relies on the heteroskedasticity and autocorrelation consistent (HAC) estimation of the long-run variance; see eicker1963asymptotic, huber1967behavior, white1980heteroskedasticity, gallant1987nonlinear, newey1987simple, andrews1991heteroskedasticity, among others.\footnote{For stationary time series, the HAC estimation of the long-run variance is the same problem as the estimation of the value of the spectral density at zero which itself has even longer history dating back to the smoothed periodogram estimators; see daniell1946discussion, bartlett1948smoothing, and parzen1957consistent.} Despite the increasing popularity of the LASSO in finance and more generally in time series empirical research, to the best of our knowledge, the validity of HAC-based inference for LASSO has not been established in the relevant literature.\footnote{See chernozhukov2019lasso for LASSO inference with dependent data and feng2019taming for an asset pricing application; see also belloni2014inference and van2014asymptotically for i.i.d.\ data; and chiang2019lasso for exchangeable arrays.} The HAC-based inference is robust to the model misspecification and leads to the valid Granger causality tests even when the fitted regression function has only projection interpretation which is the case for the projection-based definition of the Granger causality. Developing the asymptotic theory for the linear projection model with autoregressive lags and covariates, however, is challenging because the underlying processes are typically not $\beta$-mixing.\footnote{More generally, it is known that the linear transformations based on infinitely many lags do not preserve the $\alpha$- or $\beta$-mixing property.}

In this paper, we obtain the debiased central limit theorem with explicit bias correction for the sg-LASSO estimator and time series data, which to the best of our knowledge is new. Next, we establish the formal statistical properties of the HAC estimator based on the sg-LASSO residuals in the high-dimensional environment when the number of covariates can increase faster than the sample size. We show that the convergence rate of the HAC estimator can be affected by the tails and the persistence of the data, which is a new phenomenon compared to low-dimensional regressions. For the practical implementation, this implies that the optimal choice of the bandwidth parameter for the HAC estimator should scale appropriately with the number of covariates, the tails, and the persistence of the data. These results allow us to perform inference for groups of coefficients, including the (mixed-frequency) Granger causality tests.

Our asymptotic theory applies to the heavy-tailed time series data, which is often observed in financial and economic applications. To that end, we establish a new Fuk-Nagaev inequality, see fuk1971probability, for $\tau$-mixing processes with polynomial tails. The class of $\tau$-mixing processes is flexible enough for developing the asymptotic theory for the linear projection model and, at the same time, it contains the class of $\alpha$-mixing processes as a special case.

The paper is organized as follows. We start with the large sample approximation to the distribution of the sg-LASSO estimator (and as a consequence of the LASSO and the group LASSO) with $\tau$-mixing data in section (ref). Next, we consider the HAC estimator of the asymptotic long-run variance based on the sg-LASSO residuals and study the inference for groups of regression coefficients. In section (ref), we establish a suitable version of the Fuk-Nagaev inequality for $\tau$-mixing processes. We report on a Monte Carlo study in section (ref) which provides further insights about the validity of our theoretical analysis in finite sample settings typically encountered in empirical applications. Section (ref) covers an empirical application examining the Granger causal relations between the VIX and financial news. Conclusions appear in section (ref). Proofs and supplementary results appear in the appendix and the supplementary material.

\paragraph{Notation:} For a random variable $X\in\ensuremath{\mathbf{R}}$ and $q\geq 1$, let $\|X\|_q=(\ensuremath{\mathbb{E}}|X|^q)^{1/q}$ be its $L_q$ norm. 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}=\{G_g:\; g\geq 1 \}$ be a partition of $[p]$ defining groups. For a vector of regression coefficients $\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 \}$ is the support of $\beta$ and $\mathcal{G}_0 = \left\{G\in\mathcal{G}:\; \beta_{G} \ne 0\right\}$ is its group support. For $b\in\ensuremath{\mathbf{R}}^p$ and $q\geq 1$, its $\ell_q$ norm is denoted $|b|_q = \left(\sum_{j\in [p]}|b_j|^q\right)^{1/q}$ if $q<\infty$ and $|b|_\infty = \max_{j\in[p]}|b_j|$ if $q=\infty$. For $\mathbf{u},\mathbf{v}\in\ensuremath{\mathbf{R}}^T$, the empirical inner product is defined as $\langle \mathbf{u},\mathbf{v}\rangle_T = \frac{1}{T}\sum_{t=1}^T u_tv_t$ with the induced empirical norm $\|.\|_T^2=\langle.,.\rangle_T=|.|_2^2/T$. 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 part. 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 put $A_G=A_j$. For a $p\times p$ matrix $A$, let $\|A\|_\infty = \max_{j\in[p]}|A_j|_1$ be its 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$.

HAC-based inference for sg-LASSO

In this section, we cover the large sample approximation to the distribution of the sg-LASSO (LASSO and group LASSO) estimator for $\tau$-mixing processes. Next, we consider the HAC estimator of the asymptotic long-run variance based on the sg-LASSO residuals and consider the Granger causality tests. In the first subsection, we cover the debiased central limit theorem. The next subsection covers the HAC estimator and the final subsection pertains to the Granger causality tests.

Debiased central limit theorem

Consider a generic linear projection model

where $(y_t)_{t\in\ensuremath{\mathbf{Z}}}$ is a real-valued stochastic process and predictors may include the intercept, some covariates, (mixed-frequency) lags of covariates up to a certain order, as well as lags of the dependent variable. For a sample of size $T$, in the vector notation, we write

where $\mathbf{y}=(y_1,\dots,y_T)^\top$, $\mathbf{m} = (m_1,\dots,m_T)^\top$ with $m_t = \sum_{j=1}^\infty\beta_j x_{t,j}$, and $\mathbf{u}=(u_1,\dots,u_T)^\top$. We approximate $m_t$ with $x_t^\top\beta = \sum_{j=1}^p\beta_jx_{t,j}$ and put $\mathbf{X}\beta$, where $\mathbf{X}$ is $T\times p$ design matrix and $\beta\in\ensuremath{\mathbf{R}}^p$ is the unknown projection parameter. This approximation can be constructed from lagged values of $y_t$, some covariates, as well as lagged values of covariates measured at a higher frequency, in which case, we obtain the autoregressive distributed lag mixed frequency data sampling model (ARDL-MIDAS) described as

where $\phi(L) = I - \rho_1L - \rho_2L^2 - \dots - \rho_JL^J$ is a low frequency lag polynomial and the MIDAS part $\psi(L^{1/m};\beta_k)x_{t,k} = \frac{1}{m}\sum_{j=1}^m\beta_{k,j}x_{t-(j-1)/m,k}$ is a high-frequency lag polynomial; see andreou2013should and babii2020midasml. Note that when $m$ = 1 we have all data sampled at the same frequency and recover the standard autoregressive distributed lag (ARDL) model. The ARDL-MIDAS regression has a group structure where a single group is defined as all lags of $x_{t,k}$ or all lags of $y_t$ and following babii2020midasml, we focus on the sparse-group LASSO (sg-LASSO) regularized estimator.\footnote{The sg-LASSO estimator allows selecting groups and important group members at the same time.} The leading example here is the MIDAS regression involving the projection of future low frequency series onto its own lags and lags of high frequency data aggregated via some dictionary, e.g., the set of Legendre polynomials. The setup also covers what is sometimes called the reverse MIDAS, see foroni2018using and mixed frequency VAR, see ghysels2016macroeconomics, involving the projection of high frequency data onto its own (high frequency) lags and low frequency data. Such regressions, which appear in the empirical application of the paper, simply amount to a different group structure.

The sg-LASSO, denoted $\hat\beta$, solves the regularized least-squares problem

with the regularization functional

where $|b|_1 = \sum_{j=1}^p|b_j|$ is the $\ell_1$ norm corresponding to the LASSO penalty, $\|b\|_{2,1}=\sum_{G\in\mathcal{G}}|b_{G}|_2$ is the group LASSO penalty, and the group structure $\mathcal{G}$ is a partition of $[p]=\{1,2,\dots,p\}$ specified by the econometrician.

We measure the persistence of the series 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. Let $(\xi_t)_{t\in\ensuremath{\mathbf{Z}}}$ be a stochastic process and let $\mathcal{M}_t=\sigma(\xi_t,\xi_{t-1},\dots)$ be its natural filtration. The $\tau$-mixing coefficient is defined as

where the supremum is taken over $t$ and $(t_1,\dots,t_j)$. The process is called $\tau$-mixing if $\tau_k\downarrow0$ as $k\uparrow\infty$; see Lemma (ref) for the comparison of this coefficient to the mixingale and the $\alpha$-mixing coefficients. The following assumptions impose tail and moment conditions on the series of interest.

Assumption (ref) can be relaxed to non-stationary data with stable variances of partial sums at the cost of heavier notation. It allows for heavy-tailed and persistent data. For instance, it requires that either both covariates and the error process have at least $4+\epsilon$ finite moments, or that the error process has at least $2+\epsilon$ finite moments, whenever covariates are sufficiently integrable. It is also known that the $\tau$-mixing coefficients decline exponentially fast for geometrically ergodic Markov chains, including the stationary AR(1) process, so condition (ii) allows for relatively persistent data. Next, we require that the covariance matrix of covariates is invertible.

Assumption (ref) ensures that the precision matrix $\Theta = \Sigma^{-1}$ exists and rules out perfect multicollinearity. It also requires that the smallest eigenvalue of $\Sigma$ is bounded away from zero by $\gamma$ independently of the dimension $p$ which is the case, e.g., for the spiked identity and the Toeplitz covariance structures. Strictly, speaking this condition can be relaxed to $\gamma\downarrow 0$ as $p\uparrow\infty$ at the cost of slower convergence rates and more involved conditions on rates, in which case $\gamma$ can be interpreted as a measure of ill-posedness; see carrasco2007linear. The next assumption describes the rate of the regularization parameter, which is governed by the Fuk-Nagaev inequality; see Theorem (ref) and Eq. (ref).

Lastly, we impose the following condition on the misspecification error, the number of covariates $p$, the sparsity constant $s_\alpha$, and the sample size $T$.

The effective sparsity constant $\sqrt{s_\alpha} = \alpha\sqrt{|S_0|} + (1-\alpha)\sqrt{|\mathcal{G}_0|}$ is a linear combination of the sparsity $|S_0|$ (number of non-zero coefficients) and the group sparsity $|\mathcal{G}_0|$ (number of active groups). It reflects the finite sample advantages of imposing the sparse-group structure as $|\mathcal{G}_0|$ can be significantly smaller than $|S_0|$ that appears in the theory of the standard LASSO estimator. Throughout the paper we assume that the groups have fixed size, which is well-justified in time-series applications of interest.

The four assumptions listed above are needed for the prediction and estimation consistency of the sg-LASSO estimator; see Theorem (ref) in the supplementary material. Next, let $v_{t,j}$ be the regression error in $j^{th}$ nodewise LASSO regression; see the following subsection for more details. Put also $s=s_\alpha\vee S$, $S=\max_{j\in G}S_j$, where $S_j$ is the number of non-zero coefficients in the $j^{th}$ row of $\Theta$. The following assumption describes an additional set of sufficient conditions for the debiased central limit theorem.

Assumption (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, e.g., in the Gaussian graphical model. Condition (iii) is a mild restriction needed for the consistency of the sample variance of regression errors. The sparsity condition $s^2\log^2p/T\to 0$ is also used in van2014asymptotically who assume that the data are sub-Gaussian, see their Corollary 2.1. On the other hand, the rate condition on the dimension $p/\sqrt{T^{\kappa - 2}\log^\kappa p}\to 0$, is additional condition needed in our setting when regression errors are not Gaussian and may only have a certain number of finite moments. 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. Conditions of this type are standard in nonparametric literature.

Let $B = \hat\Theta\mathbf{X}^\top(\mathbf{y} - \mathbf{X}\hat\beta) / T$ denote the bias-correction for the sg-LASSO estimator, where $\hat\Theta$ is the nodewise LASSO estimator of the precision matrix $\Theta$; see the following subsection for more details. The following result describes a large-sample approximation to the distribution of the debiased sg-LASSO estimator with serially correlated non-Gaussian regression errors.

It is worth mentioning that since the group $G$ has fixed size and the rows of $\Theta$ have finite $\ell_1$ norm, the long-run variance $\Xi_G$ exists under the maintained assumptions; see Proposition (ref) in the Appendix for a precise statement of this result.

Theorem (ref) extends van2014asymptotically to non-Gaussian, heavy-tailed and persistent time series data and describes the long run asymptotic variance for the low-dimensional group of regression coefficients estimated with the sg-LASSO. One could also consider Gaussian approximations for groups of increasing size, which requires an appropriate high-dimensional Gaussian approximation result for $\tau$-mixing processes and is left for future research; see chernozhukov2013gaussian for a comprehensive review of related coupling results in the i.i.d.\ case.

Nodewise LASSO

The bias-correction term $B$ and the expression of the long-run variance in Theorem (ref) depend on the appropriate estimator of the precision matrix $\Theta = \Sigma^{-1}$. We focus on the nodewise LASSO estimator of $\Theta$, following the work of meinshausen2006high and van2014asymptotically. The estimator is based on the observation that the covariance matrix of the partitioned vector $X = (X_j,X_{-j}^\top)^\top\in\ensuremath{\mathbf{R}}\times\ensuremath{\mathbf{R}}^{p-1}$ can be written as

where $\Sigma_{j,j}=\ensuremath{\mathbb{E}}[X_j^2]$ and all other elements similarly defined. By the partitioned inverse formula, the $1^{st}$ row of the precision matrix $\Theta=\Sigma^{-1}$ is

where $\gamma_j = \Sigma_{-j,-j}^{-1}\Sigma_{-j,j}$ is the projection coefficient in the regression of $X_j$ on $X_{-j}$

and $\sigma_j^2 = \Sigma_{j,j} - \Sigma_{j,-j}\gamma_j = \ensuremath{\mathbb{E}}[v_j^2]$ is the variance of the projection error.\footnote{To ensure that the projection coefficient is well defined and does not change with the dimension of the model $p$, we can consider the limiting linear projection model and take into account the approximation error.} This suggests estimating the $1^{st}$ row of the precision matrix as $\hat\Theta_j = \hat\sigma_j^{-2}

$ with $\hat\gamma_j$ solving

and

where $\mathbf{X}_j\in\ensuremath{\mathbf{R}}^T$ is the column vector of observations of $x_j\in\ensuremath{\mathbf{R}}$ and $\mathbf{X}_{-j}$ is the $T\times (p-1)$ matrix of observations of $x_{-j}\in\ensuremath{\mathbf{R}}^{p-1}$. In the matrix notation, the nodewise LASSO estimator of $\Theta$ can be written then as $\hat\Theta = \hat B^{-1}\hat C$ with

HAC estimator

Next, we focus on the HAC estimator based on sg-LASSO residuals, covering the LASSO and the group LASSO as special cases. For a group $G\subset [p]$ of a fixed size, the HAC estimator of the long-run variance is

where $\hat{\Gamma}_k = \hat\Theta_G\left(\frac{1}{T}\sum_{t=1}^{T-k}\hat u_t\hat u_{t+k} x_tx_{t+k}^\top\right)\hat\Theta_G^\top$, $\hat u_t$ is the sg-LASSO residual, and $\hat\Gamma_{-k}=\hat{\Gamma}_k^\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 parzen1957consistent and andrews1991heteroskedasticity. Several choices of the kernel function are possible, for example, the Parzen kernel is

It is worth recalling that the Parzen and the Quadratic spectral kernels are high-order kernels that superior to the Bartlett kernel used in newey1987simple; see appendix for more details on the choice of the kernel.

Note that under stationarity, the long-run variance in Theorem (ref) simplifies to

where $\Gamma_k = \Theta_G\ensuremath{\mathbb{E}}[u_tx_tu_{t+k}x_{t+k}^\top]\Theta_G^\top$ and $\Gamma_{-k}=\Gamma^\top_k$. The following result characterizes the convergence rate of the HAC estimator pertaining to a group of regression coefficients $G\subset[p]$ based on the sg-LASSO residuals.

The first term in the inner parentheses is of the same order as the estimation error of the maximum between the estimation errors of the sg-LASSO and the nodewise LASSO. Theorem (ref) suggests that the optimal choice of the bandwidth parameter should scale appropriately with the number of covariates $p$, the sparsity constant $s$, and the dependence-tails exponent $\kappa$.\footnote{A comprehensive study of the optimal bandwidth choice based on higher-order asymptotic expansions is beyond the scope of this paper and is left for future research, see, e.g., lazarus2018har for the recent literature review and practical recommendations in the low-dimensional case.} This contrasts sharply with the HAC theory for regressions without regularization developed in andrews1991heteroskedasticity, see also li2018uniform, and allows for faster convergence rates of the HAC estimator.

High-dimensional Granger causality tests

Consider a linear projection model

where $h\geq 0$ is the horizon, $G\subset[p]$ is a group of regression coefficients of interest, $x_t=\{x_{t,j}:\; j\in G \}$ represents the series for which we wish to test the Granger causality, and $\{x_{t,j}:j\in G^c \}$ represents all the remaining information available at time $t$. For instance, $x_t$ may contain $L$ low-frequency lags of some series $(z_t)_{t\in\ensuremath{\mathbf{Z}}}$, in which case $x_t = (z_t,z_{t-1},z_{t-2},\dots, z_{t-L})^\top$. Alternatively, it may contain low and/or high-frequency lags of $(z_t)_{t\in\ensuremath{\mathbf{Z}}}$ aggregated with dictionaries, e.g., Legendre polynomials as in babii2020midasml. In both cases the dimensionality of $x_t$ small. On the other hand, the set of controls representing all the information available at time $t$ is high-dimensional. The Granger causality test corresponds to the following hypotheses

where $\beta_G = \{\beta_j:\; j\in G \}$.

It is worth mentioning our framework is based on the weakest notion of the Granger causality corresponding to the marginal improvement in time series projections due to the information contained in $x_t$. A stronger notion of Granger non-causality corresponds to the independence on the conditional mean of $y_t$ of $x_t$. Yet, even stronger version of Granger non-causality requires that $y_t|x_t,\sigma(x_{t,j}:j\in G^c) =_d y_t|\sigma(x_{t,j}:j\in G^c)$, so that the entire conditional distribution of $y_t$ does not depend on $x_t$; see also florens1982note.

Let $R$ be $r\times |G|$ matrix of linear restrictions imposed on $\beta_G$. For the Granger causality test, we set $R=I_{|G|}$, but more generally, we might be interested in testing other linear restrictions implied by the economic theory. Assuming that $R$ is a full row rank matrix, consider the debiased Wald statistics

where $A^+$ is the generalized inverse of $A$. It follows from Theorems (ref) and (ref) that under $H_0$, $W_T\xrightarrow{d}\chi^2_r$. The Wald test rejects when $W_T>q_{1-\alpha}$, where $q_{1-\alpha}$ is the quantile of order $1-\alpha$ of $\chi^2_{r}$. More generally, the linear restrictions can be extended to the nonlinear restrictions by the usual Delta method argument.

For testing hypotheses on the increasing set of regression coefficients, it might be preferable to use the non-pivotal sup-norm based statistics, see ghysels2018testing, due to the remarkable deoendence on the dimension in the high-dimensional setting; see chernozhukov2013gaussian for high-dimensional Gaussian approximations with i.i.d.\ data.

Fuk-Nagaev inequality

In this section, we describe a suitable for us version of the Fuk-Nagaev concentration inequality for the maximum of high-dimensional sums. The inequality allows for the data with polynomial tails and $\tau$-mixing coefficients decreasing at a polynomial rate. The following result does not require that the series is stationary.

The inequality describes the mixture of the polynomial and Gaussian tails for the maximum of high-dimensional sums. In the limiting case of the i.i.d.\ data, as $a\to \infty$, the dependence-tails exponent $\kappa$ approaches $q$ and we recover the inequality for the independent data stated in fuk1971probability, Corollary 4 for $p=1$. In this sense, the inequality in Theorem (ref) is sharp. It is well-known that the Fuk-Nagaev inequality delivers sharper estimates of tail probabilities in contrast to Markov's bound in conjunction with Rosenthal's moment inequality, cf. nagaev1998some. The proof relies on the blocking technique, see bosq1993bernstein, and the coupling inequality for $\tau$-mixing sequences, see dedecker2004coupling, Lemma 5. In contrast to previous results, e.g., dedecker2004coupling, Theorem 2, the inequality reflects the mixture of the polynomial and the exponential tails.

For stationary processes, by Lemma (ref) in the appendix, $B_T^2 = O(T)$ as long as $a>(q-1)/(q-2)$, whence we obtain from Theorem (ref) that for every $\delta\in(0,1)$

where $C>0$ is some finite universal constant.

Monte Carlo experiments

In this section, we aim to assess the debiased HAC-based inferences for the low-dimensional parameter in a high dimensional data setting. To that end, we draw covariates $\{x_{t,j},j\in[p]\}$ independently from the AR(1) process

The regression error follows the AR(1) process

where errors are $\epsilon_{t,j},\nu_t \sim_{i.i.d.} N(0,1)$. The vector of population regression coefficients $\beta$ has the first five non-zero entries which are drawn from Uniform$(0,4)$ and all remaining entries are zero. The sample size is $T\in\{100,1000\}$ and the number of covariates is $p\in\{10,200\}$. We set the persistence parameter $\rho = 0.6$ and focus on the LASSO estimator to estimate coefficients $\hat \beta.$ Throughout the experiment, we choose the LASSO tuning parameters using the 10-fold cross-validation, defining folds as adjacent over time blocks.

We report the average coverage (av.\ cov) and the average length of confidence intervals for the nominal coverage of 0.95 and on a grid of values of the bandwidth parameter $M_T \in \{5, 10, \dots, 60\}$, using the Parzen kernel. We estimate the long run covariance matrix $\hat\Xi$ using the LASSO residuals, denoted $\hat u_t$. We also use the nodewise LASSO regressions to estimate the precision matrix $\Theta$. The first step is to compute scores $\hat V_t = \hat u_t x_t$, where $\hat u_t = y_t - x_t^\top\hat \beta$, and $\hat\beta$ is the LASSO estimator. Then we compute the high-dimensional HAC estimator using the formuala in equation ((ref)). We compute the pivotal statistics for each MC experiment $i\in[N]$ and each coefficient $j\in[p]$ as $\text{pivot}_{j}^{(i)}$ $\triangleq$ $(\hat \beta_{j}^{(i)} + B_{j}^{(i)} - \beta )/\sqrt{\hat\Xi_{j,j}^{(i)}/T}$, where $B_{j}^{(i)} = \hat\Theta^{(i)}_j \mathbf{X}^{\top(i)} \hat{\mathbf{u}}^{(i)}/T$, and $\mathbf{\hat u}^{(i)}=\mathbf{y}^{(i)}-\mathbf{X}^{(i)}\hat\beta^{(i)}$. Then we compute the empirical coverage as

and the average confidence interval length as $\text{length}_{j}$ = $\frac{1}{N} \sum_{i=1}^{N} 2 \times 1.96 \times \sqrt{\hat\Xi_{j,j}^{(i)}/T}.$ The number of Monte Carlo experiments is set to $N$= $5000.$

We report average results over the active and inactive sets of the vector of coefficients. Table (ref) shows results for the small sample size ($T=100$) and the large sample size ($T=1000$). We find that the value of the bandwidth parameter $M_T$ should be smaller when the number of regressors $p$ is larger. For $T=100$ (small sample size), for the active set of coefficients, the best coverage is achieved when the bandwidth parameter is set at $10$ when $p=10$ and at $5$ when $p=200$. Results for the inactive set and $T=1000$ are similar. We also see that the increase in $p$ relative to $T$ leads to worse performance. Furthermore, the coverage improves when the bandwidth increases with the sample size. Lastly, as the sample size increases, the average coverage approaches the nominal level of 0.95 and the confidence interval shrinks in size. Overall, the simulation results confirm our theoretical findings.

Testing Granger causality for VIX and financial news

The CBOE Volatility Index, known as the VIX, is a popular measure of market-based expectation of future volatility and is often referred to as the “fear index”. The VIX index quotes the expected annualized change in the S&P 500 index over the following 30 days, as computed from options-based theory and current options-market data.

There is a large literature studying the theoretical and empirical properties of the VIX and it is impossible to cite only a few papers to do justice to all the outstanding research output on the topic. Focusing on Granger causal patterns, there are several studies pertaining to causality between the VIX and VIX futures. For example, bollen2017tail suggest that the VIX futures lagged the VIX in the first few years after its introduction, and show an increasing dominance of VIX futures over time. Along similar lines, shu2012causality study price‐discovery between VIX futures the spot VIX index and find evidence of a bi-directional causal pattern.

We study the causal relationship between financial news and the VIX. There is also a substantial literature on the impact of news releases on financial markets (see, e.g., andersen2003micro). Traditionally, such analysis looks at news releases and studies the behavior of asset prices pre- and post-release. News is usually quantified numerically via the surprise component measured as the difference between an expectation prior to the release and the announcement. In the age of machine learning, the characterization of news has been expanded into the textual analysis of news coverage. To paraphrase the title of gentzkow2019text, the text is treated as data. It is in this spirit that we conduct our high-dimensional Granger causality analysis between the VIX and news.

We use a data set from bybee2019structure which contains 180 news attention monthly series, all of which potentially Granger cause future US equity market volatility.\footnote{We downloaded daily VIX data from St. Louis Fed FRED database and took the end-of-month values. The FRED mnemonic for the VIX is VIXCLS. Table with the full list of series appears in Appendix (ref).} We estimate the following time series regression model

where $y_{t+1}$ is the value of the VIX at the end of month $t+1$, $\psi(L^{1/m}; \beta)y_{t}$ is a MIDAS polynomial of 22 daily VIX lags where the first lag is the last day of the month $t$, and $x_{t,k}$ is the $k$-th news attention series. Note that we only take one lag for the news attention series to simplify the model (and also an empirically justified simplification). The MIDAS polynomial of daily lags of the VIX involves Legendre polynomials of degree 3. Note that the specification is what is sometimes called a reverse MIDAS regression as mentioned earlier in the paper. Prior to estimating the regression model, we time demean the response and covariates such that the intercept is zero. We further standardize all covariates to have a unit standard deviation. The daily VIX lags are standardized before the aggregation.

We apply the sg-LASSO estimator to estimate the slope coefficients and nodewise LASSO regressions to estimate the precision matrix. To fully exploit the group sparsity of sg-LASSO, we group all high-frequency lags of daily VIX, see babii2020midasml for further details on such grouping. The news attention series are monthly and we are interested in whether the most recent news Granger causes the VIX, hence we don't apply the group structure along the time dimension. Instead, we group news attention series based on a broader theme that is available for each series, see bybee2019structure for further details. Namely, the data set contains 24 broader topics which group each of the 180 news attention series.

Main results

We report the p-values for a range of $M_T$ values for series that appear to be significant at the 1% or 5% significance level for all $M_T\in\{20,40,60\}$ values and for two kernel functions, namely Parzen and Quadratic Spectral. The sample starts January 1990 January and ends June 2017, determined by the availability of the textual analysis data. Both sg-LASSO and LASSO tuning parameters are selected via 10-fold cross-validation, defining folds as adjacent blocks over the time dimension to take into account the time series nature of the data. Similarly, we tune nodewise LASSO regressions for the precision matrix estimation.

Granger causality of news topics

The results appear in Table (ref) which contains two main row blocks reporting results for the structured sg-LASSO and unstructured LASSO estimators, and two-column blocks, reporting results for two kernel functions. Irrespective of the initial estimator and kernel function, the lagged daily VIX and the Financial crisis news series are highly significant at 1% significance level. Comparing results for the initial estimator, in the case of LASSO we see more significant predictors than for the sg-LASSO case, while the subset of significant covariates using sg-LASSO is a subset of the LASSO significant predictors. Many more series are selected by the initial LASSO estimator compared to the sg-LASSO, see Table (ref). This suggests that relevant group structures are important, and may help in recovering salient relationships in the data.

Bi-directional Granger causality

We also test whether the daily VIX Granger causes Financial crisis news series. For this we run the following MIDAS regression model

where $x_{t+1,j}$ is the Financial crisis news series. We test whether daily VIX Granger causes future values of Financial crisis news series. Note that we only need to estimate the initial initial coefficient vector, since the precision matrix remains the same. The results appear in Table (ref). They show a rather weak predictability of future news series by daily VIX suggesting a unidirectional Granger causality pattern.

Granger causal clusters of news topics

The news attention series are classified into 24 broader meta topics that group the individual news series according to a common theme. We test which group of individual news series Granger causes future VIX values. The results are reported in Table (ref). They show that the group {\it Banks} is significant at 5% significance level. This group consists of news series pertaining to news about Mortgages, Bank loans, Credit ratings, Nonperforming loans, Savings & loans, and the Financial crisis.

Conclusion

This paper develops valid inferential methods for high-dimensional time series regressions estimated with the sparse-group LASSO (sg-LASSO) estimator that encompasses the LASSO and the group LASSO as special cases. We derive the debiased central limit theorem with the explicit bias correction for the sg-LASSO with serially correlated regression errors. Furthermore, we also study HAC estimators of the long-run variance for low dimensional groups of regression coefficients and characterize how the optimal bandwidth parameter should scale with the sample size, the temporal dependence, as well as tails of the data. These results lead to the valid t- and Wald tests for the low-dimensional subset of parameters, such as Granger causality tests. Our treatment relies on a new suitable variation of the Fuk-Nagaev inequality for $\tau$-mixing processes which allows us to handle the time series data with polynomial tails. An interesting avenue for future research is to study more carefully the problem of the optimal data-driven bandwidth choice based on higher-order asymptotic expansions, see, e.g., sun2008optimal for steps in this direction in low dimensional settings.

In an empirical application we use a high-dimensional news attention series to study causal patterns between the VIX, sometimes called the fear index, and financial news. We find that almost exclusively the topic of financial crisis exhibits unidirectional Granger causality for the VIX.

\setcounter{page}{1} \setcounter{section}{0} \setcounter{equation}{0} \setcounter{table}{0} \setcounter{figure}{0}

Proofs

Next, we focus on the HAC estimator based on LASSO residuals. Note that by construction of the precision matrix $\hat\Theta$, its $j^{th}$ row is $\hat\Theta_jx_t = \hat v_{t,j}/\hat\sigma^2_j$, where $\hat v_{t,j}$ is the regression residual from the $j^{th}$ nodewise LASSO regression and $\hat\sigma^2_j$ is the corresponding estimator of the variance of the regression error. Therefore, the HAC estimator based on the LASSO residuals in Eq. (ref) can be written as

where $\hat\Gamma_{k}$ has generic $(j,h)$-entry $\frac{1}{T}\sum_{t=1}^{T-k}\hat u_t\hat u_{t+k}\hat v_{t,j}\hat v_{t+k,h}\hat \sigma_j^{-2}\hat \sigma_h^{-2}.$

Similarly, we define

where $\tilde\Gamma_{k}$ has generic $(j,h)$-entry $\frac{1}{T}\sum_{t=1}^{T-k}u_tu_{t+k}v_{t,j}v_{t+k,h}\sigma_j^{-2}\sigma_h^{-2}$ and note that the long-run variance $\Xi_G$ has generic $(j,h)$-entry $\ensuremath{\mathbb{E}}[u_tu_{t+k}v_{t,j}v_{t+k,h}]\sigma_j^{-2}\sigma_h^{-2}$.

For a stationary process $(\xi_t)_{t\in\ensuremath{\mathbf{Z}}}$, let

be its $L_1$ mixingale coefficient with respect to the canonical filtration $\mathcal{M}_0=\sigma(\xi_0,\xi_{-1},\xi_{-2},\dots)$. Let $\alpha_k$ be the $\alpha$-mixing coefficient and let $Q$ be the quantile function of $|\xi_0|$. The following covariance inequality allows us controlling the autocovariances in terms of the $\tau$-mixing coefficient as well as comparing the latter to the mixingale and the $\alpha$-mixing coefficients.

The following result shows that the variance of partial sums can be controlled provided that the $\tau$-mixing coefficients decline sufficiently fast.

Lastly, we show that in the linear regression setting, the long-run variance for a group of projection coefficients $G\subset[p]$ of fixed size exists under mild conditions. Let $\xi_t=u_t\Theta_Gx_t$, where $\Theta_G$ are rows of the precision matrix $\Theta=\Sigma^{-1}$ corresponding to indices in $G$.

We recall first the convergence rates for the sg-LASSO with weakly dependent data that will be needed throughout the paper from babii2020midasml, Corollary 3.1.

Next, we consider the regularized estimator of the variance of the regression error

where $\hat\beta$ is the sg-LASSO estimator. While the regularization is not needed to have a consistent variance estimator, the LASSO version of the regularized estimator ($\alpha=1$) is needed to establish the CLT for the debiased sg-LASSO estimator. The following result describes the converges of this variance estimator to its population counterpart $\sigma^2 = \ensuremath{\mathbb{E}}\|\mathbf{u}\|^2_T$.

Next, we look at the estimator of the precision matrix. Consider nodewise LASSO regressions in equation ((ref)) for each $j\in[p]$. Put $S=\max_{j\in G}S_j$, where $S_j$ is the support of $\gamma_j$.

Next, we first derive the non-asymptotic Frobenius norm bound with explicit constants for a generic HAC estimator of the sample mean that holds uniformly over a class of distributions. We focus on the $p$-dimensional centered stochastic process $(V_t)_{t\in\ensuremath{\mathbf{Z}}}$ and put

where $\Gamma_{k}=\ensuremath{\mathbb{E}}[V_tV_{t+k}^\top]$ and $\tilde\Gamma_{k}=\frac{1}{T}\sum_{t=1}^{T-k}V_tV_{t+k}^\top.$ Put also $\Gamma = (\Gamma_{k})_{k\in\ensuremath{\mathbf{Z}}}$ and let $\langle.,.\rangle$ be the Frobenius inner product with corresponding Frobenius norm $\|.\|$. The following assumption describes the relevant class of distributions and kernel functions.

Condition (ii) describes the smoothness (or order) of the kernel in the neighborhood of zero. $\varsigma=1$ for the Bartlett kernel and $\varsigma=2$ for the Parzen, Tukey-Hanning, and Quadratic spectral kernels, see andrews1991heteroskedasticity. Since the bias of the HAC estimator is limited by the order of the kernel, it is typically not recommended to use the Bartlett kernel in practice. Higher-order kernels with $\varsigma>2$ do not ensure the positive definiteness of the HAC estimator and require additional spectral regularization, see politis2011higher. Condition (iv) describes the class of autocovariances that vanish rapidly enough. Note that if (iv) holds for some $\bar\varsigma$, then it also holds for every $\varsigma<\bar\varsigma$ and that if (ii) holds for some $\tilde\varsigma>\varsigma$, then it also holds for $\tilde\varsigma=\varsigma$. The covariance condition in (iv) can be justified under more primitive moment and summability conditions imposed on $L_1$-mixingale/$\tau$-mixing coefficients, see Proposition (ref) and andrews1991heteroskedasticity, Lemma 1. The following result gives a nonasymptotic risk bound uniformly over the class $\mathcal{G}$ and corresponds to the asymptotic convergence rates for the spectral density evaluated at zero derived in parzen1957consistent.

{\scriptsize

}