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.
88,498 characters · 9 sections · 68 citation commands
Inference in Non-stationary High-Dimensional VARs
\doublespacing
In this paper we construct an inferential procedure for Granger causality in high-dimensional non-stationary vector autoregressive (VAR) models. Investigating causes and effects in time series models has a long and rich history, dating back to the seminal work of granger1969investigating. The statistical assessment of the directional predictability among two (or blocks of) time series can have important consequences for decision making processes. Applications of Granger causality range over from macroeconomics, finance, network theory, climatology and even the neuroscience. For instance, the evidence of causality between money and gross domestic product is a long-debated issue in the macroeconomic literature. This was first discussed in sims1990inference and stock1989interpreting and it is still argued to this day, see e.g., miao2020high. Financial applications of Granger causality include, among others, billio2012econometric who study the connectedness among monthly returns of hedge funds, banks, broker/dealers, and insurance companies. hecq2021granger build Granger causality networks from high-dimensional vector autoregressive models, describing the dynamic volatility spillovers among a large set of stock returns. Many applications are also found in climate science, for instance in trying to understand and disentangle the causes of climate change. Among others, stern2014anthropogenic investigate causality between greenhouse gas emissions and temperature. In neuroscience, Granger causality is employed to understand mechanisms underlying complex brain function and cognition, with examples in the field of functional neuroimaging seth2015granger, friston2013analysing.
With the increased availability of larger and richer datasets, these causality concepts have recently been extended to a high-dimensional setting where they can benefit from the inclusion of many more series within the information set. In the original concept of Granger causality granger1969investigating, conditioning on a given information set plays a central role.\footnote{Throughout the paper we focus on Granger causality in mean. While there are clearly many other forms of causal or predictive relations possible, Granger causality in mean is the most prominent and therefore the focus of our attention. It should therefore be understood that whenever we mention Granger causality in the paper, we refer to Granger causality in mean.} After all, Granger causality is a study of predictability. Only by considering predictability given a specific conditioning set, is possible to attach some sort of causal meaning to the outcome (where the nature of that causality is still up for debate, of course). To identify direct `truly causal' links between variables would require one to condition upon all possible variables that may be related to the variables of interest. Otherwise, any discovered relation may simply be an artifact of omitted variable bias which would invalidate a causal interpretation. Indeed, Granger himself envisioned the information set as “all the knowledge in the universe available at that time" granger1980testing. While this concept cannot be operationalized with any finite amount of data, the availability of increasingly high-dimensional datasets, along with the econometric techniques to analyze them, provide a great opportunity for turning Granger causality from `just predictability' into a concept to which at least some form of causal interpretation can be attached.
Recently, hecq2021granger (HMS henceforth) proposed a test for Granger causality in high-dimensional vector autoregressive models. This combines dimensionality reduction techniques based on penalized regression such as the lasso of tibshirani1996regression, with the post-double selection procedure of belloni2014inference designed to guarantee uniform asymptotic validity of the post-selection least squares estimator. However, HMS assumed stationarity of all the time series considered. This is a typical assumption in much of the literature on regularization methods, in particular when inference is considered. While the literature on estimation and forecasting with high-dimensional non-stationary processes is growing SmeekesWijler20, this is not the case for inference due to the complexities arising with unit roots and cointegration, which already have severe effects in low-dimensional settings.
Working with (data possibly transformed to yield) stationary time series avoids these complications in the asymptotic analysis and allows to invoke standard (Gaussian) limit theory, thereby enabling the use of standard inferential procedures such as $t$ and $F$-tests. On the other hand, it assumes prior knowledge of the order of integration of all the time series entering the model. This prior knowledge is usually acquired via unit root tests such as the augmented Dickey-Fuller (ADF) test dickey1979distribution. However, these tests are sensitive to their specifications, such as the inclusion of a deterministic time trend in the regression equation and the choice of the lag length. Even when performing what seems to be the same test in, for example, different R packages, it may actually lead to different outcomes SmeekesWilms20. That does not even take into account that there are many different unit root tests, none of which is clearly preferred to the others, as well as how specific data features (and their treatment) such as seasonality (adjustments), structural breaks and outliers may affect these tests. As such, the outcome of a unit root test is ambiguous at best; even more so if one also takes into account their well-documented low power cochrane1991critique on the one hand and the accumulation of the probability of obtaining false positives through multiple testing -- especially in high dimensions -- on the other hand.
But even if one did know the true orders of integration, transformations such as differencing to achieve stationarity are not innocuous. In particular, information about long-run relations such as cointegration is deleted from the series. Yet the error correction mechanisms at play for the movement towards long-run equilibria may well induce (Granger) causal relations that are not apparent anymore in the transformed data. While this is not necessarily a problem if one is interested in the transformed data as such (for example the growth rate of economic series may be an object of interest in themselves), in many applications we are interested in the level of the series and the transformations are only done to avoid issues with non-stationarity. It would therefore be very beneficial to practitioners to have high-dimensional inference methods available that are robust to (unknown) orders of (co)integration.
In this paper we develop a method which allows for testing Granger causality in high-dimensional VAR models, irrespective of the (co)integration properties of the time series in the VAR. We avoid any bias coming from unit root and cointegration pre-testing and instead use the VAR in levels directly to perform inference on the (Granger causality) parameters of interest. The procedure we design builds on the work of toda1995statistical who first considered a simple lag augmentation of the system, which they showed provides asymptotically normal estimators irrespective of potential unit roots and cointegration dolado1996making. The same approach has also recently been used to develop inference on impulse response functions inoue2020uniform,montiel2021local that is robust to unit roots. We modify the approach outlined above by confining lag augmentation to only the variable(s) tested as Granger causing, instead of adding lags of all variables. We show that this modification, which makes it suitable for application to high-dimensional VARs, does not affect the asymptotic properties. In addition, we modify the post-double selection procedure of HMS to prevent the possibility of spurious regression, thereby extending its uniform validity to data with potential (co)integrated time series.
The remainder of the paper is organized as follows: Section (ref) introduces the model, the Granger causality testing and our lag-augmented post-double-selection approach. The theoretical properties of our method are studied in Section (ref). In Section (ref) we investigate finite-sample performance through a simulation study, while Section (ref) proposes a data-driven method to find a sensible upper bound for the lag length in the VAR. Section (ref) uses the proposed testing framework to investigate the causes and effects of economic uncertainty in the context of the FRED-MD dataset. Section (ref) concludes. In Appendix (ref) the theory underlying the lag augmentation and the asymptotic properties is developed. Introductory Lemmas and complementary results to the following two appendices are also given. Appendix (ref) is devoted to the validation of high-level assumptions, while additional empirical results are presented in Appendix (ref).
A few words on notation. For any $n$-dimensional vector $\boldsymbol x$, we let ${\left\lVert\boldsymbol x\right\rVert}_p = \left(\sum_{i=1}^n |x_i|^p \right)^{1/p}$ denote the $\ell_p$-norm. For any index set $S \subseteq \{1, \ldots, n\}$, let $\boldsymbol x_{S}$ denote the sub-vector of $\boldsymbol x_t$ containing only those elements $x_i$ such that $i \in S$. $|S|$ denotes the cardinality of the set $S$. We use $\xrightarrow{p}$ and $\xrightarrow{d}$ to denote convergence in probability and distribution, respectively.
In this section we propose our model, our strategy for lag augmentation of the Granger causality tests, and the post-double-selection procedure needed to achieve uniformly valid inference. Section (ref) first presents the model and lag-augmented Granger causality test. Next, (ref) sets out the post-double-selection procedure.
Let $\boldsymbol z_1,\ldots,\boldsymbol z_T$ be a $K$-dimensional multiple time series process, where $\boldsymbol z_t=(y_t, x_t, \boldsymbol w_t)^{'}$. Here $y_t$ is the series we would like to test being Granger caused, $x_t$ the potentially Granger causing series and $\boldsymbol w_t$ is a $K-2$ dimensional vector of controls constituting the information set. We allow for $K$ to be large, potentially larger than (and growing with) the sample size $T$. We assume $\boldsymbol z_t$ is generated by a VAR($p$) process as
where $\boldsymbol A_1,\ldots,\boldsymbol A_{p}$ are $K\times K$ parameter matrices and $\boldsymbol u_t$ is a $K\times 1$ vector of error terms.
Note that Assumption (ref)((ref)) allows for the time series to have unit roots and be cointegrated. We specifically allow elements of $\boldsymbol z_t$ to be integrated of order $d$: $I(d)$ for $d=0,1,2$ and possibly cointegrated of order $d,b$: $CI(d,b)$ with $0<b\leq d$. We formulate more specific assumptions on (co)integration properties in Section (ref) and Appendix (ref).
We are interested in testing the null hypothesis of Granger non-causality in mean between the Granger causing series $x_t$ and the Granger caused $y_t$, conditional on all the series in $\boldsymbol w_t$. For the moment we assume the lag-length $p$ in (ref) to be known; we shall further elaborate on data-driven ways to select $p$ in Section (ref). Also, in order to ease the notation, we omit both the intercept and any polynomial time trend from the model; the results we derive easily extend to those cases as well. It is convenient to introduce the following stacked notation. Let $\boldsymbol X_{-p} = (\boldsymbol x_{-1}, \ldots, \boldsymbol x_{-p})$ denote the $T\times p$ matrix containing the $p$ lags of the Granger causing variable $x_t$, where $\boldsymbol x_{-j}$ is the vector containing the observations corresponding to its $j$-th lag.\footnote{To keep $T$ observations for the lags, one can replace the missing values for $t\leq 0$ with zeros. Alternatively, the vectors are shortened to contain $T-p$ observations only. Both approaches are equivalent asymptotically, and for notational simplicity in the following we do not explicitly distinguish between them.} Similarly, we define $\boldsymbol Y_{-p}$ containing the $p$ lags of the Granger caused variable $y_t$ and $\boldsymbol W_{-p}$ for the conditioning set such that $\boldsymbol W_{-p}$ is a $T\times (K-2)p$ matrix. Then, our equation of interest for the Granger causality testing can be stated as
where for notational simplicity we define $\boldsymbol V = (\boldsymbol Y_{-p}, \boldsymbol W_{-p})$ as the matrix containing all control variables and $\boldsymbol \delta = (\boldsymbol \delta_1', \boldsymbol \delta_2')'$ its coefficients.
Testing for no Granger causality is then equivalent to testing the following null hypothesis:
To account for the potential unit roots in the system, we follow the approach pioneered by toda1995statistical and dolado1996making of augmenting the regression of interest with redundant lags of the variables. However, in contrast to the existing approaches, we only augment the lags of the Granger causing series $x_t$. That is, we consider the lag-augmented regression
where $\boldsymbol X_{-(p+d)} = (\boldsymbol X_{-p}, \boldsymbol x_{-(p+1)}, \ldots, \boldsymbol x_{-(p+d)})$ contains the $p+d$ lags of $x_t$. Here $d$ represents the maximum order of integration one suspects the series are having. Note that in fact $\beta_{p+j} = 0$ for all $j\geq 1$ such that $\boldsymbol \beta_{+}=(\boldsymbol \beta', \boldsymbol 0_d')'$, as we are adding redundant variables. It should be clear that by testing whether the first $p$ elements of $\boldsymbol \beta_{+}$ -- those corresponding to $\boldsymbol \beta$ in (ref) -- we can perform the same test of Granger non-causality as in the original setup.
By adding “free” lags of $x_t$, we in essence allow this variable to “difference itself” into the correct order to remove the unit roots. As a simple illustration, consider the following DGP where $x_t$ may be $I(1)$:
Testing whether $\beta=0$ by estimating this regression directly, is complicated as the limit distribution of the least squares estimator and the corresponding test changes depending on whether $x_t$ is $I(1)$ or $I(0)$. By adding a redundant lag we can write (ref) as
suggesting that regressing on $x_{t-1}$ and $x_{t-2}$ is equivalent to regressing on $u_{2,t-1}$ and $x_{t-2}$. Moreover, if $u_{2,t}$ were observed, we could test $\beta=0$ directly using its regression coefficient. In Appendix (ref) we show formally that not only the equivalence in (ref) continues to hold in the presence of higher order lags as well as control variables, but also that testing for Granger causality in the lag-augmented regression is in fact equivalent to testing on the appropriately transformed variables. This, in turn, allows us to retrieve the asymptotic normality of the least squares estimator regardless of the order of integration, provided that the lag order $p$ of the VAR and maximum order of integration $d$ are correctly specified.
The main difference compared to the original results of toda1995statistical is that we only augment the Granger causing series $x_t$. While this difference is not of importance when $K$ is small, it opens the door for high-dimensional applications where $K$ is large, and even larger than $T$. As we cannot estimate such regressions with least squares anymore, we combine lag-augmentation with the post-double-selection framework of HMS to construct Granger causality tests in high-dimensional models that are robust to unknown orders of integration. We describe the method in the next section.
Assuming that $\boldsymbol \delta$ is sparse, in the sense that many elements of the vector are equal to zero, one might consider estimating (ref) with a technique that does variable selection such as the lasso, and then re-estimating the equation including only the selected control variables. However, doing so we run into the problems with inference after model selection leeb2005model, where even if selection is done consistently, the (diminishing) probability of omitting a relevant variable causes a sufficiently large bias to prevent uniform convergence of the post-selection least squares estimator to a Gaussian limit distribution. To circumvent this problem belloni2014inference introduced the post-double-selection (PDS) method that not only selects relevant variables on the outcome variable, but also on the treatment variable. This double selection reduces the probability of admitting relevant variables to such an extent that it does not affect the asymptotic distribution anymore.
The PDS framework was used by HMS to develop high-dimensional tests for Granger causality in sparse VAR models. In this section we show how to adapt the post-double selection framework of HMS to the unit root non-stationary framework with lag augmentation. In this framework we first perform a set of initial regressions -- to be estimated with variable selection techniques such as the lasso -- of the dependent variable plus the explanatory variables of interest (here the $p$ lags of $x_t$) on all other variables. To properly account for potential unit roots and avoid spurious regression in these initial regressions, we need to slightly adapt the setup of HMS.
Let $\boldsymbol X_{-p,\backslash\{j\}}$ denote the matrix $\boldsymbol X_{-p}$ from which the $j$-th column, corresponding to the $j$-th lag of $x_t$, has been removed. Let $\boldsymbol Z_{-p} =(\boldsymbol X_{-p}, \boldsymbol Y_{-p}, \boldsymbol W_{-p})$ denote the matrix of all (non-augmented) variables and $\boldsymbol Z_{-p,\backslash\{j\}}=(\boldsymbol X_{-p,\backslash\{j\}}, \boldsymbol Y_{-p}, \boldsymbol W_{-p})$. Then we consider the following regressions in the first step:
In contrast to HMS, who follow belloni2014inference by only regressing $\boldsymbol y$ and $\boldsymbol x_{-j}$ on the controls $\boldsymbol Y_{-p}$ and $\boldsymbol W_{-p}$, we also add all other lags of $x_t$. This is done to avoid that the error terms in (ref) contain unit roots and the regressions become spurious.
Although the concept is tricky with integrated variables, the coefficients $\boldsymbol \eta_j$ in (ref) can be thought of as best linear predictors. Indeed, as the error terms $\boldsymbol e_j$ are $I(0)$, these coefficients are ensured to exist and can be defined as the probability limits of the respective least squares estimators while keeping the dimension $K$ fixed.
As will be formalized in Assumption (ref), we assume that the coefficients $\boldsymbol \eta_j$ are sparse. Let $\mathcal{S}_j = \{m>p: \eta_{m,j} \neq 0\}$, $j=0,1, \ldots, p$, be the sets of active variables in (ref).\footnote{We exclude the variables corresponding to the lags of $x_t$ from $\mathcal{S}_j$ as these will always be included in the second stage.} Then, we need that the cardinality of these sets is small; that is, smaller than $T$ and at most growing at a slow rate of $T$. In fact, the union over all these sets, $\mathcal{S} := \bigcup_{j=0}^{p} \mathcal{S}_j$ is our set of interest. This set represents all variables needed to control for when regressing $\boldsymbol y$ on $\boldsymbol X_{-p}$, as they either have non-zero coefficients in (ref) or are correlated with $x_t$. In fact, variables would need to have both properties to cause omitted variable bias if they were not included in the final regression. By aiming to recover $\mathcal{S}$, we therefore have two opportunities to select relevant variables. This is sufficient to reduce the probability of missing them to be asymptotically negligible.
Let $\boldsymbol V_{\mathcal{S}}$ denote the matrix consisting of only those columns of $\boldsymbol V$ that corresponds to the selected variables in $\mathcal{S}$, and $\boldsymbol \delta_{\mathcal{S}}$ the corresponding parameter vector. Our goal of the PDS regression is then to recover the regression
in the second stage, and base inference on this.
To obtain an estimate for the set $\mathcal{S}$, one can use the lasso or any of its relatives that similarly ensure variable selection. The lasso tibshirani1996regression simultaneously performs variable selection and estimation of the parameters in (ref) by solving the following minimization problems
where $\lambda$ is a non-negative tuning parameter determining the strength of the penalty. Here, we follow the framework of HMS and use the Bayesian information criterion (BIC) in selecting the tuning parameter, coupled with a penalty lower bound ensuring a maximum of selected variables per estimated equation (see Remark (ref) for details). Minimizing an information criterion (IC) in order to determine an appropriate data-driven $\lambda$ is one way to deal with dependent data (see HMS for an overview of other methods and their finite sample behaviors).
The first step of our PDS method is then to estimate each of the regressions in (ref) using a penalized regression technique such as in (ref). Let $\hat{\mathcal{S}}_j = \{m>p: \hat\eta_{m,j} \neq 0\}$ represent the corresponding sets of active variables retained in each of these regressions, and let $\hat{\mathcal{S}} = \bigcup_{j=0}^{p} \hat{\mathcal{S}}_j$ denote the set of all active variables. Then, we base our inference on the second-stage regression
which may now be treated as if no selection took place, and can simply be estimated by least squares. As we will show in the next section, the OLS estimator converges uniformly to a normal distribution, which in turn makes standard tests such as the Wald test or the LM test applicable with their regular limiting distributions.
Algorithm (ref) describes the main steps of our post-double selection lag-augmented (PDS-LA) Granger causality test implemented through the Lagrange Multiplier (LM) or Wald test.
In this section we present the main theoretical result. We show the post-selection, lag-augmented, least squares estimator $\hat \boldsymbol \beta_{p}$ is asymptotically Gaussian uniformly over the parameter space. Therefore, tests for Granger causality are $\chi^2$ distributed. For the PDS-LA method to deliver uniformly valid inference, a set of assumptions is necessary. Among others, sparsity needs to be assumed on the high-dimensional vector $\boldsymbol \delta$. Also, a restricted (sparse) eigenvalue condition needs to be assumed on the (scaled) Gram matrix, bounding away its smallest eigenvalue over a subset (cone) of $\mathbb{R}^K$. This condition essentially guarantees that over a sufficiently large subset of the parameter space the Gram matrix is well behaved and thus invertible. An empirical process bound is also required. This states that the (scaled) process $\boldsymbol Z_{-p}'\boldsymbol u$ uniformly concentrates around zero.
These conditions are standard in the literature to prove the consistency of the lasso. However, when nonstationary processes are considered, more refined results are required. In particular, scaling becomes more complicated as time series of different orders need to be scaled with different rates. In addition, when time series are cointegrated, we first need to separate the stochastic trends from the stationary components before the appropriate scaling can be applied.
For example, suppose that $\boldsymbol z_{t}$ can be written as the cointegrated system
where $\boldsymbol u_t$ is $I(0)$ and $\boldsymbol A$ and $\boldsymbol B$ are $n\times r$ matrices with $r < n$. Then define $\boldsymbol \zeta_t := \boldsymbol z_t \boldsymbol Q$ where $\boldsymbol Q := (\boldsymbol B, \boldsymbol A_{\perp})'$ and $\boldsymbol A_\perp$ is the $n \times n - r$ orthognal complement of $\boldsymbol A$, such that $\boldsymbol A_\perp' \boldsymbol A = \boldsymbol 0$. Then we can write
Now the first $r$ variables in $\boldsymbol \zeta_t$ are $I(0)$, while the remaining ones are $I(1)$. Further details on constructing such a matrix $\boldsymbol Q$ can be found in, e.g., lutkepohl2005new.
Given the existence of such a matrix $\boldsymbol Q$, we can without loss of generality assume that we may partition $\boldsymbol \zeta_t$ as $\boldsymbol \zeta_t = (\boldsymbol \zeta_{0,t}, \boldsymbol \zeta_{1,t}, \boldsymbol \zeta_{2,t})'$, where $\boldsymbol \zeta_{i,t}$ contains the $I(i)$ variables. We then consider a scaling matrix
where $n_i$ corresponds to the number of $I(i)$ variables in $\boldsymbol \zeta_t$. Finally, we let $\boldsymbol G_T := \boldsymbol Q \boldsymbol D_T$ as the final `scaling and rotation' matrix.
Condition ((ref)) bounds the empirical process with high probability. Such deviation bounds for non-stationary time series can be found in, among others, smeekes2021automated, mei2022lasso and wijler2022restricted. Condition ((ref)) is standard and assumes compactness of the parameter space of the vector $\boldsymbol{\beta}$ which in turn implies the boundedness.
Condition ((ref)) supposes consistency of the first-stage estimator. Consistency is thereby established using an error bound or “oracle inequality". Such inequality gives the rate of convergence of the estimator as a function of: the tuning parameter $\lambda$, the deviation bound, the cardinality $s$ of the active set and the restricted (sparse) eigenvalue. As such, this condition is intimately related to ((ref)), ((ref)) and ((ref)). For stationary time series many results exist under a variety of settings; see e.g. masini2022regularized for VAR models or adamek2022lasso for general time series models. Nonstationary time series are treated in smeekes2021automated, mei2022lasso and wijler2022restricted.
Condition ((ref)) requires sparsity of the population parameters and the estimator. Sparsity of the first-stage estimator is needed in our framework as we perform OLS on the selected variables from the first-stage regressions. If the selected variables are not sparse enough, too many variables will be selected for OLS to be feasible. For simplicity we work under the assumption of exact sparsity; however, at the expense of more complicated notation and proofs this can be relaxed to approximate sparsity following belloni2014inference, where it is assumed that the exact sparse model is a (good) approximation to the true DGP, or weak sparsity following adamek2022lasso, where many non-zero but small coefficients are allowed.
Condition ((ref)) requires that for sufficiently sparse vectors, the eigenvalues of the subset of the Gram matrix corresponding to their nonzero support do not decrease to zero too fast. While this is a standard (and easily verifiable) assumption for stationary time series (see e.g., adamek2022lasso and masini2022regularized), time series with unit roots require more care. Indeed, smeekes2021automated, mei2022lasso and wijler2022restricted show for unit root regressors that the standard rate of $T^{-2}$ coming from the $\boldsymbol D_T$ scaling still results in at least a factor of $\bar{s}_T^{-1}$ in $\kappa_{T,\min}$. As we allow $\kappa_{T,\min}$ to depend on the sample size, this can be accommodated, as long as the interplay of the rates in condition ((ref)) is satisfied. This condition links the sparsity $\bar{s}_T$, the tails thickness of the empirical process $\bar{\gamma}_T$ and the minimum eigenvalue $\kappa_{T,\min}$. The exact rates allowed for are a compromise between the number of moments existing, the strength of the dependence, the growth rate of the dimension $K$ and the sparsity of the parameters.\footnote{For an illustration of the complexities of these relations, we refer to Figure C.1 in adamek2022lasso which visualises the feasible combinations with regards to lasso consistency.}
We can now state our first, and main, theoretical result. This establishes that the first stage regressions allow for sufficiently accurate variable selection and estimation such that the second stage is not affected by the variable selection performed.
Theorem (ref) establishes the asymptotic equivalence of the estimator in the feasible second-stage regression (ref) based on the estimated active set $\hat{\mathcal{S}}$, and the estimator in the infeasible second-stage regression (ref) based on the true, unobserved, active set $\mathcal{S}$. Note that the equivalence is only established for the estimators of the coefficients of the lags of the Granger causing variable -- minus the augmented lags -- this is however all that is needed to establish the validity of the Granger causality tests. Before stating the result about the asymptotic distribution, we need another assumption on the asymptotic behaviour of the relevant variables. For a finite number of relevant variables -- measured by $s$ in Assumption (ref)((ref)) -- this assumption follows directly from well-known results in the unit root and cointegration literature; see Appendix (ref). We here state the assumption more generally to also accommodate sparsity that increases with the sample size.
Assumption (ref) appears rather abstract with the sequences $\phi_{T,\cdot}$, $\kappa_{T,\cdot}$ and $\gamma_{T,\cdot}$ bounded in a complicated nonlinear way. However, note that the actual assumption that is required is that if we combine the blocks in (a)-(d) appropriately -- first through (i)--(iii), then via (ref) -- the terms become negligible. If we assume that the number of variables $n$ is finite, which in our setting follows from the sparsity $s$ not growing with the sample size, we can express $\phi_{T,\cdot}$, $\kappa_{T,\cdot}$ and $\gamma_{T,\cdot}$ in well-known rates needed to obtain limiting distributions of sums of products of integrated variables, see e.g. hamilton1994time for an overview. We show in Appendix (ref) how Assumption (ref) can be verified in this case.
One way to extend the result to a growing number of variables would be through a Gaussian approximation theorem; such an approach is considered in, e.g., zhang2019identifying and smeekes2021automated. This is a relatively crude approach that puts significant limitations on the allowed growth rate of $s$. However, as the growth rate of $s$ is anyway only allowed to be limited compared to the sample size $T$, such an approach would be feasible here without imposing strong additional restrictions (unlike in the general case such as for establishing a result like the deviation bound (ref)((ref)) where it would significantly restrict the allowed growth rate of $K$). With this assumption in place we can now state our second theoretical result, which establishes the limit distribution of the Granger causality tests.
We now evaluate the finite-sample performance of our proposed PDS-LA-LM Granger causality test. We consider the following Data Generating Processes (DGPs) in first differences inspired by kock2015oracle:
with $a=0.3$. The diagonal VAR(1) for DGP1 allows the sparsity assumption to be met. Instead, for DGP2 the coefficients decrease with exponential pace departing from the main diagonal and hence although the farthest coefficients are small, the exact sparsity assumption is not met. We report simulations for Granger causality tests from the first variable to the second variable. Therefore, also when integrated out to nonstationary, DGP1 automatically satisfies the null of no Granger causality from unit 2 to 1, however DGP2 does not. Therefore, for the power analysis of both DGP1 and DGP2 we set the coefficient in position $(2,1)$ equal to 0.2. We set the same coefficient equal to zero for DGP2 for the size analysis.
Simulations are reported for different types of covariance matrices of the error terms. We employ a Toepliz-version for calculating the covariance matrix as $\Sigma_{i,j}=\rho^{|i-j|}$, where $(i,j)$ refer to row $i$, column $j$ of the matrix $\Sigma_{u}$. We cover two scenarios of correlation: $\rho=(0,0.7)$.
The lag length is fixed to $p=2$, while we employ a double ($d=2$) augmentation of the Granger causing variable. Having $p\geq 2$ guarantees that no spurious results occur in the selection steps. Following the recommendation in HMS, we employ the BIC in selecting the tuning parameter $\lambda$ for the lasso.
Table (ref) reports the size and power of the PDS-LA-LM test out of 1000 replications. We use different combinations of time series length $T=(50,100,200,500,1000)$ and number of variables in the system $K=(10,20,50,100)$. All the rejection frequencies are reported using a burn-in period of fifty observations.
Our PDS-LA-LM test shows good performance in terms of size and (unadjusted) power for all DGPs considered. The setting of no correlation is handled remarkably well by all DGPs and only moderate size distortion is visible in large systems for small samples. Whenever high correlation of errors is present, sizes are still in the vicinity of 5% for DGP1 where the sparsity assumption is met. However, we notice how for DGP2, for which the sparsity assumption is not met, some residual size distortion remains visible even in large systems. However, the power of the test is always increasing with the sample size $T$ for all the considered cases.
Up until this point, we considered the lag length $p$ as given. In reality, this is typically not the case. In this section we propose a simple, data-driven method to estimate $p$. Standard techniques for tuning the lag length such as information criteria or sequential testing fail when applied directly to the high-dimensional VAR. The standard approach with penalized regression methods would be to set $p$ as a generous upper bound, and let the method decide which lags are needed. This however provides complications for our approach as a large $p$ means many first-stage regressions need to be done, with the potential of selecting too many variables. On top of that, we need to augment the second stage with $d$ additional lags.
It is therefore important for our approach to have a reasonable data-driven selection of the lag length. We are essentially looking for an `informative upper bound' on the lag length; while mild over-specification is not a problem, under-specification breaks the lag augmentation, and must be avoided. We therefore base selection on univariate autoregressions for all time series in our dataset, and we apply an information criterion to those. This approach is motivated through the final equations representation of a VAR model, which implies that a VAR($p$) with $K$ variables generates individual ARMA models with maximal orders $(Kp,(K-1)p)$. As such, lag length selection based on individual autoregressions is likely to yield lag lengths larger than $p$ zellner1974time,cubadda2009studying. On the other hand, cubadda2009studying observed that individual autoregressions often need much smaller lag lengths than the large orders implied by the final equations representation, thus making it plausible that the orders found are not overly conservative.
We implement this method as follows. First, we estimate the autoregressions
by OLS. Then, letting $\hat\omega_i = \frac{1}{T} \sum_{t=1}^T \hat{\varepsilon}_{i,p,t}^2$, we set $\hat{\boldsymbol \varOmega}_p = \text{diag}(\hat\omega_1, \ldots, \hat\omega_K)$ and choose the $p$ that minimizes
where $C_T$ takes the standard values for well-known information criteria, e.g. $C_T = \ln T$ yields BIC$^*$, while $C_T = 2$ yields AIC$^*$. Note that next to the dimensionality reduction offered by running individual autoregressions instead of a VAR, further reduction is achieved by letting $\hat{\boldsymbol \varOmega}$ be diagonal.
We investigate the performance of this method in a small simulation study. We simulate VAR models using DGP2 with $\rho=0$ as in Section (ref) for the combinations of $K=(10,20,50,100)$ and $T=(50,100,200,500,1000)$.\footnote{We investigated DGP1 and $\rho=0.7$ as well. The method showed similar or better performance there. Results are available on request.} Knowing the true value of $p=2$, we apply the model selection procedure on a grid of 10 values for $p$.
We show the results in Figure (ref) for AIC$^*$ and BIC$^*$. Both criteria succeed in barely ever underestimating the lag length. Moreover, BIC$^*$ is quite accurate, only overestimating the lag length by more than one when $T=1000$, where larger values of $p$ are not problematic. AIC$^*$ is more liberal but still performs well in smaller samples, where it matters most. This simple way of selecting lag lengths therefore works remarkably well in finite examples, and is therefore recommended in combination with the PDS-LA method.
The role of economic uncertainty is a hot topic in macroeconomics. In particular, it is a matter of much debate whether economic uncertainty should be seen as an exogenous shock, causing economic conditions such as business cycle fluctuations, or an endogenous response to economic fundamentals ludvigson2021uncertainty. A complicating issue is that there is no universal definition or measurement of uncertainty. In this application we address both issues: we investigate whether uncertainty (Granger) causes or is caused by other variables in the economy, and whether different measurements contain the same information.
Since the seminal paper of bloom2009impact economic uncertainty has been at the forefront of macroeconomic debate. bloom2009impact built a measure of uncertainty from the Chicago Board Options Exchange (CBOE) S&P 100 Volatility Index (VXO), arguing that stock market volatility expectations are a good proxy for overall economic uncertainty. The FRED-MD dataset of mccracken2016fred has introduced the series VXOCLSx in its September 2015 vintage. This series is constructed by splicing a synthetic historical VXO series obtained from Nicholas Bloom’s website and VXOCLS from FRED.\footnote{In its December 2021 vintage, FRED-MD removed VXOCLSx and replaced it with VIXCLSx as the former has been discontinued from the source. VIX is a similar measure of implied volatility based on the S&P 500. While I made analysis focus on a pre-Covid vintage incorporating VXOCLSx, in Appendix (ref) We extend our analysis to a recent post-Covid vintage with VIXCLSx included.}
An alternative, very popular measure of uncertainty was proposed by baker2016measuring, who constructed the Economic Policy Uncertainty (EPU) index based on newspaper coverage frequency. This accounts for the frequency with which certain strings of keywords related to economic uncertainty appear in the 10 leading U.S. newspapers. The Federal Reserve Economic Division (FRED) itself issues a Global Economic Policy Uncertainty Index (GEPUCURRENT) which is a GDP-weighted average of national EPU indices for 20 countries. However, the FRED-MD dataset does not contain any EPU uncertainty index. This naturally raises the question whether EPU can serve as a better measure of economic uncertainty within FRED-MD or VXOCLSx is “good enough". If EPU and VXOCLSx measure the same content, one would expect that adding EPU to FRED-MD would not lead to much predictive power for EPU if VXOCLSx is accounted for in the information set.
Our PDS-LA method can thus be used in this context with a twofold purpose. First, by estimating Granger causal relationships between FRED macroeconomic variables and the uncertainty index and counting the significant number of outgoing (index$\to$FRED) and incoming (FRED$\to$index) links, we contribute in the understanding of whether such indexes are respectively exogenous sources of business cycle (index$\to$FRED) or rather endogenous responses to economic fundamentals (FRED$\to$index). Second, by investigating the links between a second index (EPU) and the FRED variables conditional on VXOCLSx in the dataset, we can investigate whether EPU measures information not already present in VXOCLSx.
In our analysis we are particularly going to pay attention to the effect of how to treat potential nonstationarity in the data. The FRED-MD dataset comes with a detailed appendix mccracken2016fred and Matlab/R routines which allow not only to clean the data from missing values and outliers, but also to take the appropriate transformations to render all the time series stationary. This presents the practitioner with an easy and popular solution to deal with non-stationarity, however there are two potential issues with this approach. First, not every time series can be clearly categorized regarding the order of integration, and there are several variables for which a case can be made for two different orders. Indeed, as illustrated by SmeekesWijler20 and SmeekesWilms20, unit root tests may give ambiguous results and depending on the type of test employed, a different classification may arise. Second, even with a correct classification, differenced time series loose long-run information on any cointegrating relations, which may change the Granger causal relations. Our method provides an alternative way of performing the tests directly on the levels of the data, thereby avoiding the issue of transformations altogether.\footnote{To run these analyses we used the authors R package HDGCvar available at \url{https://github.com/Marga8/HDGCvar}}.
For the Economic Policy Uncertainty index we gather the data directly from baker2016measuring.\footnote{The data is freely available on the website of Scott R. Baker, Nick Bloom and Steven J. Davis at \url{https://www.policyuncertainty.com/index.html}.} Specifically, we use the three component index which combines (i) news coverage about policy-related economic uncertainty, (ii) tax code expiration data and (iii) economic forecaster disagreement.\footnote{We also repeated the analysis using the News-based Policy Uncertainty index as described in baker2016measuring but the results are nearly identical to those reported here.} For details on how the index is computed we refer to the given reference. As the EPU index is only available starting January 1985, we accordingly use the FRED-MD data from January 1985. We split the analysis considering two different endpoints of the sample. First we consider the series until November 2019, thus intentionally excluding from the sample both the Covid pandemic and the war in Ukraine, but including the 2008 financial crisis. This makes for a total of 117 variables (when EPU is included) and 419 observations. In Appendix (ref) we extend the analysis to a recent vintage which includes the aforementioned crises and spans until September 2022.
First, we do not add yet US-EPU to the information set and in Figure (ref) and (ref) we instead just investigate the predictive relations to and from VXOCLSx with all the macroeconomic series of FRED-MD using PDS-LA-LM.
At significance level $\alpha=0.01$, eight macroeconomic series are found to Granger cause VXOCLSx while the other way around VXOCLSx Granger-causes $41$ macroeconomic series. Before adding US-EPU to the dataset and investigate the changes, we repeat the analysis, this time applying the recommended stationary transformations (ST) from FRED-MD and we use the PDS-LM test of HMS instead of PDS-LA-LM to investigate the same relations.
The difference between Figure (ref), (ref) and Figure (ref), (ref) is quite striking, especially for the case VXOCLSx$\rightarrow$ FRED. In fact, the number of significant connections drops dramatically from $39$ to just two and some differences can be found in the connections for the FRED$\rightarrow$VXOCLSx too. To zoom in on the connections and the difference between the lag-augmented (LA) and stationary transformed (ST) cases, in Figure (ref) we loosen the p-value threshold up to $10\%$ and group the variables by their FRED-MD sector classification. The bars now indicate the $p$-values of the test, with a full bar equaling a $p$-value of 0, and an empty bar indicating a $p$-value above 10%.
The results again suggest how stationary transforming variables can have a profound impact on the inference performed, leading to very different pictures of Granger causality. It is particularly noticeable that sectors thought to be highly affected by uncertainty, such as the Output sector, barely has any significant connections left when testing for Granger causality from VXOCLSx. Similarly, the absence of causality from the Stocks category to VXOCLSx is surprising, given the latter's construction.
Next, we add US-EPU to the dataset and again repeat the analysis where now the focus is connections from and to US-EPU with all the macroeconomic series of FRED-MD using PDS-LA-LM. Importantly, now the results will be conditioned on VXOCLSx which remains in the information set.
In Figure (ref) we only find five connections when it come to Granger causality relations from FRED-MD series to US-EPU. Vice-versa, in Figure (ref) there are 13 causal paths from US-EPU to the other macroeconomic series. The results are quantitatively similar to FRED$\leftrightarrow$VXOCLSx in Figure (ref), (ref), though less pronounced. Both US-EPU and VXOCLSx are Granger caused by the S&P500 and S&P:indust but the remaining connections are not equal. Importantly, VXOCLSx is found to be Granger causal for US-EPU. Given the network US-EPU$\rightarrow$FRED is conditional on VXOCLSx, the non-overlapping connections found in Figure (ref) that are not occurring in Figure (ref) are new connections that could make US-EPU an interesting uncertainty index to add to FRED-MD. A cross-check reveals that only UMCSENTx is new, all others connections were already uncovered using VXOCLSx. This would seem to give credit to FRED-MD in using VXOCLSx as main index of economic uncertainty without the need to further include US-EPU. On the other hand, finding the connections from US-EPU despite the presence of VXOCLSx in the information set, indicates that US-EPU may still add additional predictive information about the variables predictable by VXOCLSx, suggesting that could still be worth including it in an empirical analysis.
As before, in Figure (ref), (ref) we repeat the analysis, this time around however we apply the recommended stationary transformations (ST) from FRED-MD and we use hecq2021granger PDS-LM test to investigate the same relations. Again, a decrease in the number of connections in both directions is evident from the network results. In Figure (ref) we again zoom in on the $p$-values and group the results by sector. The red bar refers to the VXOCLSx series which is found to be Granger causal for EPU in the analysis in levels. Interestingly, we find the opposite result with stationary transformed setting. Again, the differences in the relations found indicate that transforming variables to stationarity may delete useful information about predictability.
Our analysis suggests that uncertainty accompanying a wide variety of global events, if measured in terms of expected volatility on the financial market, is primarily a cause rather than an effect of variations in the economic activity. This is in line with the recent work of ludvigson2021uncertainty who find uncertainty about financial markets to be a source of output fluctuations. In addition, our analysis shows that even if different measures of uncertainty may be closely related, they are not identical and the choice of which one to use can affect further analysis. It may therefore be worthwhile to consider multiple measures of uncertainty in empirical applications.
We propose an inferential procedure for Granger causality testing in high-dimensional non-stationary VAR models which avoids any knowledge or pre-tests about integration or cointegration in the data. To do so we adapt the toda1995statistical approach of augmenting the lag length of the system and we show that by reducing this augmentation to only the variables of interest for the test, we are able to minimize parameter proliferation in high dimensions. We develop a post-double selection LM test which is based on penalized least squares estimators to partial-out those variables having no influence on the variables tested on while safeguarding from omitted variable bias using a double-selection mechanism.
We prove that the augmentation of the Granger causing variables has no effect on the null hypothesis tested, yet provides an automatic differencing mechanism letting the OLS estimator have standard asymptotic results. Also, we extend the relevant assumptions needed for the post-double selection estimator to work in the context of potential unit roots. We derive the asymptotics of the post-selection, augmented estimator, showing it attains standard asymptotic normality hence allowing for a valid test with standard $\chi^2$ limiting distribution.
Our proposed test shows good finite sample properties over different DGPs. We also give practical recommendations on both the optimal augmentation $d$ and on how to estimate the lag-length $p$. We argue that $d=2$ lags is the optimal augmentation in order to take into account possible I(2) as well as near I(2) variables that could compromise the size of the test. In order to estimate the lag length $p$ we propose to reduce the original VAR to a diagonal VAR. This reduces the high-dimensional system to a sequence of low-dimensional autoregressions, to which an information criteria is applied in order to select the correct lag length.
Finally, we investigate how our test performs in practice by analysing the causes and effects of economic uncertainty. Using the FRED-MD dataset directly, without needing to apply their recommended stationary transformations, we compare two different ways of measuring economic uncertainty (VXO and EPU) and their relationship with all macroeconomic variables within the dataset. We also compare the analysis with the stationary transformation case, highlighting how such transformation can profoundly impact the results and give a very different picture of the causal structure. Our results suggest that uncertainty is primarily a cause rather than an effect of variations in the economic activity.