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.
92,211 characters · 14 sections · 113 citation commands
Granger Causality Testing in High-Dimensional VARs: a Post-Double-Selection Procedure
\doublespacing
Economics, statistics and finance have seen a rapid increase of applications involving time series in high-dimensional systems. Central to many of these applications is the vector autoregressive (VAR) model that allows for a flexible modelling of dynamic interactions between multiple time series. In this paper we develop a simple method to test for Granger causality in high-dimensional VARs (HD-VARs) with potentially many variables.
Many financial applications consider Granger causality analysis, especially for constructing high-dimensional networks. Networks of financial firms' intedependencies are investigated in basu2015network, gao2017efficient, demirer2018estimating and barigozzi2019nets. Similarly, spillovers and contagion among stock returns are investigated in networks using Granger causality analysis in lin2017regularized, vyrost2015granger and corsi2018measuring.
Most of the econometric literature has traditionally been focused on allowing for high dimensionality in VARs through the use of factor models bernanke2005measuring,chudik2016theory or Bayesian methods banbura2010large. For instance billio2012econometric develops measures of connectedness to assess systemic risk propagation among institutions in the financial system using principal component analysis and Granger causality networks. Recent years have seen an increase in regularized, or penalized, estimation of sparse VARs based on popular methods from statistics such as the lasso tibshirani1996regression and elastic net zou2005regularization, which impose sparsity by setting a (data-driven) selection of the coefficients to zero.
Compared to factor models, such sparsity-seeking methods have often an advantage of interpretability, as in many economic applications, it appears natural to believe that the most important dynamic interactions among a large set of variables can be adequately captured by a relatively small -- but unknown -- number of `key' variables. As such, the use of these methods for estimating HD-VAR models has also increased significantly in recent years, see e.g. nicholson2017varx, basu2019low, billio2019bayesian, wilms2018algorithm, korobilis2019adaptive).
Regularized estimation theory for high-dimensional time series and VAR models is now well established, see among others song2011large, basu2015regularized, kock2015oracle, davis2016sparse, medeiros2016l1, audrino2018oracle and masini2019regularized and wong2020lasso; kock2020penalized provide a recent review. However, performing inference on HD-VARs, such as testing for Granger causality, still remains a non-trivial matter. As is well known, performing inference after model selection (post-selection inference) is complicated as the selection step invalidates `standard' inference where the uncertainty regarding the selection is ignored leeb2005model. Complexities introduced by the temporal and cross-sectional dependencies in the VAR mean that most recently developed post-selection inference methods are not automatically applicable.
Most existing literature on Granger causality testing in HD-VARs therefore has so far not considered post-selection inferential procedures. wilms2016predictive propose a bootstrap Granger causality test in HD-VARs, but do not account for post-selection issues. Similarly, skripnikov2018joint investigate the problem of jointly estimating multiple network Granger causal models in VARs with sparse transition matrices using lasso-type methods, but focus mostly on estimation rather than testing. song2019better focus on statistical procedures for testing indirect/spurious causality in high-dimensional scenarios, but consider factor models rather than regularized regression techniques. lin2017regularized consider high-dimensional multi-block VARs derived from a two-blocks recursive linear dynamical system and use a maximum likelihood (ML) estimator for Gaussian data. In order to obtain the ML estimates for the system transition matrices and the precision matrix, respectively the lasso and graphical lasso on the residuals are iterated until convergence. krampe2018bootstrap develops bootstrap techniques for sparse VAR models combining a model-based bootstrap procedure and the de-sparsified lasso (see van2014asymptotically) to perform inference on the autoregressive parameters. chaudhry2017uncertainty look at de-biased estimators as in javanmard2014confidence, for Gaussian and sub-Gaussian VAR processes with a focus on Granger-causality and control of the false discovery rate.
In this paper we build on the post-double-selection approach proposed by belloni2014inference, to develop a valid post-selection test of Granger causality in HD-VARs. The finite-sample performance depends heavily on the exact implementation of the method. In particular, the tuning parameter selection in the penalized estimation is crucial. We therefore perform an extensive simulation study to investigate the finite-sample performance of the different ways to set up the test in order to be able to give some practical recommendations. In addition, we investigate the construction of networks of realized volatilities using a sample of 30 financial stocks modeled as a vector heterogeneous VAR corsi2009simple. We are able to demonstrate how our approach allows for obtaining much sharper conclusions than standard low-dimensional VAR techniques.
The remainder of the paper is as follows: Section (ref) introduces the high-dimensional VAR model and Granger causality tests. In Section (ref) we propose our estimation and inferential framework. Section (ref) establishes the asymptotic properties of our method and discusses the assumptions required for the theory to hold. Section (ref) reports the results of the Monte Carlo simulations. We apply our method in Section (ref) to construct volatility spillover networks. Section (ref) concludes. Proofs and supplemental results can be found in the appendix.
A few words on notation. For any $n$-dimensional vector $\bm x$, we let ${\left\lVert\bm x\right\rVert}_p = \left(\sum_{i=1}^n {\left\lvertx_i\right\rvert}^p \right)^{1/p}$ denote the $\ell_p$-norm. For any index set $S \subseteq \{1, \ldots, n\}$, let $\bm x_{S}$ denote the sub-vector of $\bm x_t$ containing only those elements $x_i$ such that $i \in S$. ${\left\lvertS\right\rvert}$ denotes the cardinality of the set $S$. We use $\xrightarrow{p}$ and $\xrightarrow{d}$ to denote convergence in probability and distribution, respectively.
Loosely speaking, the notion of Granger causality captures predictability given a particular information set granger1969investigating,granger1980testing. If the addition of variable $X$ to the given information set $\Omega$ alters the conditional distribution of another variable $Y$, and both $X$ and $\Omega$ are observed prior to $Y$, then $X$ improves predictability of $Y$, and is said to Granger cause $Y$ with respect to $\Omega$. granger1969investigating originally envisioned the information set $\Omega$ “be all the information in the universe” (p. 428), which is of course not a workable concept. Yet clearly the choice of information set has a major effect on the interpretation of the finding of (non-)Granger causality, as discussed in granger1980testing. In particular, spurious Granger causality from $X$ to $Y$ may be found when both $X$ and $Y$ are Granger caused by $Z$, but $Z$ is omitted from $\Omega$. As such, one might want to include as many potentially relevant variables in the information set as possible in order to avoid finding spurious causality due to omitted variables, thereby moving as much as possible towards the universal information set envisioned by Granger. However, conditioning on so many variables leads to obvious problems of high-dimensionality rendering many standard statistical techniques invalid.
In this paper we focus on testing Granger causality in mean using linear models, in which setup the VAR model is the natural tool to investigate this problem. However, to enlarge the information set means estimating a VAR with an increasing number of variables. The number of parameters in a VAR increases quadratically with the number of time series included; an unrestricted VAR($p$) has $K^2p$ coefficients to be estimated, where $K$ is the number of series and $p$ is the lag-length. As the time series dimension $T$ is typically fairly small for many economic applications, the data do not contain sufficient information to estimate the parameters and consequently standard least squares and maximum likelihood methods suffer from the curse of dimensionality, resulting in estimators with high variance that overfit the data.
Let $\bm y_1,\ldots,\bm y_T$ be a $K$-dimensional multiple time series process, where $\bm y_t=(y_{1,t},\ldots,y_{K,t})^{'}$ is generated by a VAR($p$) process
where for notational simplicity we assume the variables have zero mean; if not they can be demeaned prior to the analysis, or equivalently a vector of intercepts is added. $\bm A_1,\ldots,\bm A_{p}$ are $K\times K$ parameter matrices and $\bm u_t$ is a martingale difference sequence (mds) of error terms. We consider weakly stationary VAR models, as formalized in Assumption (ref) below.
In the VAR model (ref) we are interested in testing whether variables in the set $J$ Granger cause variables in the set $I$ in mean, conditional on all the other variables, where $J, I \subset \{1,\ldots, K\}$ and $J \cap I = \emptyset$. Let $N_I = {\left\lvertI\right\rvert}$ and $N_J = {\left\lvertN_J\right\rvert}$ denote the number of variables in $I$ and $J$ respectively. We describe our proecdure here in general form for testing blocks of variables. For any sets $S_1, S_2 \subseteq \{1, \ldots, K\}$ of variables define the best linear predictor in $L_2$-norm of $\bm y_{S_1,t}$ given $\bm x_{S_2,t-1}^{(p)} = (\bm y_{S_2,t-1}^\prime, \ldots, \bm y_{S_2, t-p}^\prime)^\prime$ as $\mathcal{P} (\bm y_{S_1,t}|\bm x_{S_2,t-1}^{(p)}) = \bm \varGamma^{*} \bm x_{S_2,t-1}^{(p)}$, where $\bm \varGamma^* = \min_{\bm \varGamma} \mathbb{E} \left[\left.{\left\lVert\bm y_{S_1,t} - \bm \varGamma \bm x_{S_2, t-1}\right\rVert}_2^2 \right. \right]$. Then we say that $\bm y_{J,t}$ does not Granger cause $\bm y_{I,t}$ conditionally on $\bm x_{J^c,t}$ if
for any value of $\bm x_{J^c,t}$. In other words, conditional on $\bm x_{J^c,t}$, addition of the lags of $\bm y_{J,t}$ to the information set does not improve predictability of $\bm y_{I,t}$. Note that Granger (non-)causality as defined in (ref) is a property of the population. In the VAR (ref) this means that testing for Granger causality can be done via testing the joint significance of the blocks of coefficients in the matrices $\bm A_1, \ldots, \bm A_p$ corresponding to the impact of variables $J$ on $I$.
To illustrate, consider (ref) with $p=1$ lag, and assume without loss of generality that the variables in $\bm y_t$ are ordered such that $\bm y_t = \left(\bm y_{I,t}^\prime, \bm y_{J, t}^\prime, \bm y_{-(I \cup J), t}^\prime\right)^\prime$, where $-(I \cup J))$ refers to all variables not in $J$ or $I$. Then we can write
where $\bm A$ is partitioned conformably with the blocks in $\bm y_t$. In this case, the best linear predictors in (ref) are given by
For any arbitrary value of $\bm y_{t-1}$, these can only coincide if $\bm A_{I, J} = \bm 0$. Hence, the null hypothesis of no Granger causality from $J$ to $I$ in the VAR($1$) model can be formulated in terms of $\bm A_{I, J} = \bm 0$. This is easily extended to $p>1$ by simply testing if the $(I, J)$-block of all $p$ lag matrices is equal to zero.
In the remainder of the paper, we will be working with a stacked representation of (ref) for the variables in $I$. Specifically, let $\bm Y = \left(\bm y_{p+1}, \ldots, \bm y_T\right)^\prime$ and let $\bm y_{I} = \operatorname{vec}\left(\bm Y_{I}\right)$ denote the $N_I \times 1$ stacked vector containing all observations corresponding to the variables in $I$. Similarly, let $\bm u_{I} = \operatorname{vec}(\bm U_{I})$, where $\bm U = \left(\bm u_{p+1}, \ldots, \bm u_T\right)^\prime$. Let $\bm X = \left(\bm x_{p}^{(p)}, \ldots, \bm x_{T-1}^{(p)}\right)^\prime$ and $\bm X^{\otimes} = \bm I_{N_I} \otimes \bm X$, while defining the stacked parameter vector $\bm \beta = \operatorname{vec}((\bm A_1, \ldots, \bm A_p)^\prime)$. Then we can write
where $\bm X_{GC}^{\otimes} = \bm I_{N_I} \otimes \bm X_{GC}$, and $\bm X_{GC} = \left(\bm x_{J, p}^{(p)}, \ldots, \bm x_{J, T-1}^{(p)}\right)^\prime$ contains those columns of $\bm X$ corresponding to the potentially Granger causing variables in $J$; $\bm X_{-GC}$ and $\bm X_{-GC}^{\otimes}$ are then defined similarly but containing the remaining variables.\footnote{Note that if $I = \{i\}$ for one particular value of interest, then (ref) simply corresponds to a single equation from the VAR in (ref).} Testing for no Granger causality is then equivalent to testing $H_0: \bm \beta_{GC}=\bm 0$ against $H_1: \bm \beta_{GC}\neq \bm 0$.
Define $N_J = {\left\lvertJ\right\rvert}$ and $N_I = {\left\lvertI\right\rvert}$. Note that $\bm \beta_{-GC}$ has $\left(K - N_J \right) \times N_I \times p$ elements, which we assume large through having a large number of variables $K$. On the other hand, throughout the paper we assume that $N_J$, $N_I$ and $p$ are small, or more precisely, fixed when sample size increases to infinity. As $\bm \beta_{GC}$ has $N_{GC} = N_J \times N_I \times p$ elements, these are also implied to be fixed. While theoretically it is possible to consider an increasing number of elements in $\bm \beta_{GC}$ (see Remark (ref) for details), it would not be required for typical applications. $J$ and $I$ are under the researcher's control and in most applications it is natural to consider a small number of variables of interest; often both $J$ and $I$ will only consist of a single variable, as in our application.
For $p$ it may appear more restrictive to assume it small. However, large $p$ in univariate regressions or small systems often arise from neglected dynamics with omitted variables hecq2016univariate. As our HD-VAR attempts to include many more variables than typical small systems, we hope to alleviate the omitted variable issue, and thereby also directly making smaller $p$ much more realistic. Of course, $p$ is generally unknown in practice. However, in many applications it is possible to give a reasonable (and small) upper bound on $p$, which is sufficient for our algorithm. If not, $p$ has to be estimated. We discuss two ways in the next section.
In this section we introduce our inferential procedure to the Granger causality tests in high-dimensional VARs. We first discuss the lasso, which we use in the initial stage to select relevant variables. Next we discuss how naive use of the lasso introduces post-selection problems for inference, and we propose our algorithm to remedy this.
As $\bm \beta$ is high-dimensional when $Kp$ is large relative to $T$, least squares estimation is not appropriate, and a structure must be imposed on $\bm \beta$ to be able to estimate it consistently. We assume sparsity of $\bm \beta$; that is, we assume that $\bm \beta$ can accurately be approximated by a coefficient vector with a (significant) portion of the coefficients equal to zero.
The sparsity assumption validates the use of variable selection methods, thereby reducing the dimensionality of the system without having to sacrifice predictability. For a general $n$-dimensional vector of responses $\bm y$ and $n \times M$-dimensional matrix of covariates $\bm X$, the (weighted) lasso simultaneously performs variable selection and estimation of the parameters by solving
where $\lambda$ is a non-negative tuning parameter determining the strength of the penalty, and $\{w_m\}_{m=1}^{M}$ are non-negative weights corresponding to the parameters in $\bm \beta$. For the standard lasso the weights are either equal to one, or equal to zero (if this parameter should not be penalized). The notation $\hat{\bm \beta}(\lambda)$ highlights that the solution to the minimization problem depends on $\lambda$, which has te be selected as well (see Section (ref)). When no confusion can arise, we simply write $\hat{\bm \beta}$.
One may also consider the adaptive lasso zou2006adaptive with parameter-specific weights $w_j$ in (ref) based on an initial estimation of $\bm \beta$, which is able to delete more irrelevant variables. However, for our purpose such oracle properties are not very relevant; we wish to eliminate the effects of the other “nuisance” variables on the relation between the variables tested for Granger causality, but we do not need to identify which of these nuisance variables matter.
Theoretical properties of lasso estimation in stable VAR models have now been studied extensively. We here non-exhaustively mention some of the key results for our setting; see kock2020penalized for a thorough review. kock2015oracle derive oracle properties of the adaptive lasso for VAR models. basu2015regularized establish restricted eigenvalue conditions for VAR models and show their sufficiency for estimation consistency. medeiros2016l1 relax the Gaussianity assumptions of these papers by considering conditionally heteroskedastic errors, and demonstrate that the adaptive lasso retains oracle properties in time series settings. Finally, masini2019regularized derive bounds on estimation errors in approximately sparse VAR models under very general conditions, allowing for heavy tails and dependence in the error terms. In particular, they show that several commonly used volatility processes in financial research satisfy these assmuptions, thereby formally establishing the suitability of the lasso for many financial applications of VAR models.
One might be tempted to simply perform the (adaptive) lasso as in (ref) on (ref), setting $w_{GC} = 0$, and then testing whether $\bm \beta_{GC}=0$, potentially after re-estimating the model by OLS on only the selected variables. However, this ignores the fact that the final, selected, model is random and a function of the data. The randomness contained in the selection step means the post-selection estimators do not converge uniformly to a normal distribution, as the potential omitted variable bias from omitting (weakly) relevant variables in the selection step is too large to maintain uniformly valid inference.
In a sequence of papers leeb2005model, Leeb and P\"otscher address these issues, showing that distributions of post-selection estimators only converge point-wise but not uniformly in the parameter space to normal distributions. Therefore, “standard” asymptotics fail to deliver a proper approximation of finite-sample behavior due to the presence of small, hard to detect parameters, whose omitted variable bias is too large to ignore asymptotically. As such, post-selection based on oracle properties is only appropriate if one a priori rules out small parameters conditions (via beta-min conditions, see e.g. van2011adaptive) thus obtaining a sharp separation of non-zero from zero coefficients. This is typically far too strong to be reasonable in applications, and methods explicitly accounting for selection are required.
Several approaches to valid post-selection inference, also referred to as honest inference, have been developed in recent years based on various philosophies, such as simultaneous inference across models berk2013valid, inference conditional on selected models lee2016exact, or debiasing (desparsifying) the lasso estimates van2014asymptotically,zhang2014confidence. We focus on the double selection approach developed by Belloni, Chernozhukov and co-authors; see e.g. belloni2014high for an overview. This approach is tailored for the lasso, easy to implement, and can be extended to dependent data.
belloni2014uniform develop a post-double-selection approach to construct uniform inference for treatment effects in partially linear models with high-dimensional controls using the lasso. Two initial lasso estimations of both the outcome and the treatment variable on all the controls are performed, and a final post-selection least squares estimation is conducted of the outcome variable on the treatment variable and all the controls selected in at least one of the two steps. The double variable selection step substantially diminishes the omitted variable bias and ensures the errors of the final model are (close enough to) orthogonal with respect to the treatment. The authors proved uniform validity of the procedure under a wide range of DGPs, including heteroskedastic and non-Gaussian errors.
chernozhukov2019lasso extend the analysis of estimation and inference for highly-dimensional systems in regressions, allowing for (weak) temporal and cross-sectional dependency. Regularization techniques for dimensionality reduction are applied iteratively in the system and the overall penalty is jointly chosen by a block multiplier bootstrap procedure. Oracle properties and bootstrap consistency of the test procedure are derived. Furthermore, simultaneous valid inference is obtained via algorithms employing least square or least absolute deviation after (double) lasso selection step(s). Although our approach is closely related to that of chernozhukov2019lasso, it differs in a number of ways. Our method is simpler and faster to implement as it does not rely on bootstrap methods. Also, chernozhukov2019lasso focus on general systems of equations and general ways of performing inference, which is different from our specific focus on Granger causality and VAR models. Third, we consider a different set of assumptions to establish the validity of our method, where we specifically focus on the relevance of these assumptions for applications in financial econometrics.
We here describe how to implement the post-double-selection procedure in a VAR context. Let $\bm x_{GC,j}$, $j=1,\ldots, N_X$, where $N_X = p N_J$, denote the $j$-th column of $\bm X_{GC}$ and consider the partial regressions:
where $\bm \gamma_j, \; j = 0,\ldots, N_X$, are the best linear prediction coefficients\footnote{Note that Assumption (ref)((ref)) implies that $\left(\mathbb{E} \bm x_{-GC,t-1} \bm x_{-GC,t-1} \right)^{-1}$ and hence $\left(\mathbb{E} \bm X_{-GC,t-1}^{\otimes} \bm X_{-GC,t-1}^{\otimes \prime} \right)^{-1}$ exist.}
where $\bm X_{-GC,t}^{\otimes} = \bm I_{N_I} \otimes \bm x_{-GC,t-1}$. As the errors $\bm e_0, \ldots, \bm e_{N_X}$ are orthogonal to $\bm X_{-GC}$, partialling out the effects of these variables would allow for a valid test of Granger causality. Of course, (ref) and (ref) are still high-dimensional and cannot be estimated by least squares. However, we can select the relevant variables from lasso estimation of (ref) and (ref) and collect all these for the final estimation of $\bm y_{I}$ on $\bm X_{GC}^{\otimes}$ plus only those relevant variables.
Intuitively, this works because to cause omitted variable bias on the coefficients of $\bm X_{GC}$, a particular variable in $\bm X_{-GC}$ must have a nonzero coefficient in both (ref) and one of the regressions in (ref). If its coefficient is zero in (ref), it has no effect on $\bm y_{I}$ and is therefore not wrongfully omitted. If it has a zero coefficient in all regressions in (ref), it is not correlated with any variables of interest, and omitting it will not result in a bias. By including all variables that are selected in at least a single of these regressions, we essentially allow for “one free mistake” by the lasso in failing to select a relevant variable. That is, omitted variable bias will only occur if the lasso fails to select a relevant variable in both regressions simultaneously. As the probability of this occurring decreases quadratically, this is sufficient to be negligible asymptotically and allow for uniformly valid inference. We provide a formal justification in Section (ref).
We now state the details of our algorithm which executes the post-double-section along the lines described above, and conclude this section with some remarks.
Appropriate selection of the lasso tuning parameter $\lambda$ in (ref) is crucial to achieve good performance. Many different data-driven methods exist giving wildly varying results. We provide a systematic comparison of several popular methods discussed in the literature in a simulation study. To the best of our knowledge, this is the first such comparison in the context of post-selection inference. We now introduce the methods considered in our study.
One option is to minimize an information criterion (IC) to determine an appropriate data-driven $\lambda$. Let $\hat{S}(\lambda) = \left\{m \in \{1,\ldots,Kp\}: {\left\lvert\hat{\beta}_m (\lambda)\right\rvert} > 0 \right\}$ denote the set of active variables in the lasso solution for a given $\lambda$. For a generic response vector $\bm y$ and predictor matrix $\bm X$, the value $\lambda^{IC}$ is found as
where $C_T$ is the penalty specific to each criterion. We consider the Akaike information criterion (AIC) by akaike1974new with $C_T=2$, the Bayesian information criterion (BIC) by schwarz1978estimating with $C_T=\ln(T)$, and the Extended Bayesian information criterion (EBIC) by chen2008extended with $C_T = \ln(T) + 2 \gamma \ln(Kp)$ with $\gamma=0.5$ proposed by chen2012extended who argue that BIC fails to select the correct variables when the number of parameters is larger than the sample size.
An alternative approach is to plug in estimates of theoretically optimal values bickel2009simultaneous,belloni2013least,belloni2011square. The lasso requires that $\lambda\geq c{\left\lVert\bm X'\bm u\right\rVert}_{\infty}/T$ for some constant $c>0$ with “high probability”. The central limit theorem motivates a Gaussian approximation where one chooses $\lambda^{th}=\frac{2c\hat{\sigma}}{\sqrt{T}}\Phi^{-1}\bigg(1-\frac{\alpha}{2N}\bigg)$ for a small $\alpha=o(1)$, where $\Phi^{-1}(\cdot)$ is the inverse of standard Gaussian cumulative distribution function and $\hat{\sigma}$ is an estimate the variance of $\bm u$. In this paper we set $\alpha =0.05/\ln(T) $ and $c = 0.5$, while we follow belloni2012sparse in the estimation of $\sigma$. Specifically, we obtain an initial (conservative) estimate by least squares estimation of $\bm y$ on the five most correlated regressors. This estimate is then updated iteratively, for details see belloni2012sparse.
Perhaps the most popular way to choose the tuning parameter is cross-validation (CV), although CV is not always appropriate in the time series setup without modifications bergmeir2018note. To estimate the tuning parameter with CV in a time series setup (TSCV) we use an expanding window out-of-sample forecasting scheme and minimize its squared forecasting error. The rolling window is set up with $80\%$ of the sample for training and $20\%$ for testing. Cross-validation is appealing since it does not require any plug-in estimates, however, as observed in chetverikov2020cross it typically yields small values of $\lambda$ thus still gaining fast convergence rate but at the price of less variable selection.
In this section we derive the asymptotic properties of our method. We first present and discuss our general high-level assumption under which the properties are derived, and then state our main results.
Assumption (ref) is a high-level assumption that allows for much flexibility on the underlying DGP and the used estimators in the first step. We now discuss each part in turn. Part (a) assumes that the minimum eigenvalue of $\bm \varSigma$ is bounded. This is required for application of lasso methods, as well as for the inverse covariance matrix $\bm \varSigma^{-1}$ and the projection coefficients in (ref) and (ref) to exist. Part ((ref)) assumes that a central limit theorem and weak law of large numbers hold. Essentially this require that the process is sufficiently well-behaved in terms of moments and dependence allowed. Although for convenience we assume martingale difference errors in Assumption (ref), ((ref)) holds under much weaker conditions such as mixing errors; see e.g. davidson1994stochastic.
Part ((ref)) is closely related to ((ref)), but additionally controls the tail behavior of the empirical process. Results of this kind are standard in the lasso literature and can be derived using a variety of tail bounds depending on the properties of the random variables of interest, see e.g. kock2015oracle and medeiros2016l1 for results relevant to VAR and time series models. Of particular interest for financial applications, masini2019regularized show that this condition is satisfied for VAR models with general weakly dependent erorrs that include many popular multivariate volatility models. The boundedness assumption in (ref) is not very restrictive, and with $N_{GC}$ fixed follows directly if the parameter space of $\bm \beta$ is a compact set.
Part ((ref)) imposes an appropriate consistency rate on the predictions coming from the first-stage estimator. Such prediction consistency is a standard result for lasso estimators; in particular, wong2020lasso obtain it for a very general class of VAR models allowing for conditional heteroskedasticity and dependence in the error terms. adamek2020lasso derive consistency of the lasso under misspecified time series models, and show that their setting covers (among others) the first-step regressions of the relevant predictors in $\bm X_{GC}$ on the other regressions, which are inherently misspecified in a VAR setup due to the missing lags; see their Remark 3 for further details.
Next to consistency, we also require sparsity of the DGP and the estimator, as controlled by part ((ref)). The assumption of exact sparsity in the DGP for the initial regressions can be relaxed to approximate sparsity as in belloni2014inference. For the sake of expositional clarity we do not work under that assumption here but stick to the simpler exact sparsity. 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. Sparsity of lasso estimators is analysed in belloni2013least, while kock2015oracle and medeiros2016l1 provide results for adaptive lasso for time series. Importantly, we do not require consistent model selection; the selection method used is allowed to make “persistent” mistakes, allowing for both variables to be incorrectly included and relevant variables to be missed, as long as the estimator remains sufficiently sparse and consistency is guaranteed. Unlike belloni2014inference, we allow for the order of sparsity of the estimator to differ from the true sparsity thereby opening the way for conservative selection procedures.
Given the assumptions above, the eigenvalue assumption in ((ref)) becomes almost superfluous, as it is generally needed to establish ((ref)) and ((ref)) for lasso-type estimators; see e.g. belloni2013least and medeiros2016l1 for details. It requires that for sufficiently sparse vectors, the eigenvalues of the subset of the Gram matrix corresponding to their non-zero support do not decrease to zero too fast. Such assumptions are standard in the lasso literature in various guises as restricted eigenvalue conditions, and can typically be derived by making similar conditions on the population covariance matrix $\bm \varSigma_{-GC, -GC}$ coupled with a convergence result of the Gram matrix $\bm X_{-GC}^\prime \bm X_{GC}$ to $\bm \varSigma_{-GC, -GC}$. basu2015network, masini2019regularized and wong2020lasso establish the plausibility of such restricted eigenvalue conditions for various VAR models. We state the condition here explicitly as it is needed directly in the proofs.
Finally, note that the restrictions on tail behavior (via $\gamma_T$), sparsity (via $\bar{s}_T$) and minimum eigenvalues (via $\phi_{T,\min}$) are meaningless if no rates on these sequences are imposed. Part ((ref)) therefore is the key part which connects all assumptions with explicit rates needed for the validity of the PDS method. The restrictions here represents a trade-off between sparsity, thickness of tails and minimum eigenvalues. For example, if, as often assumed $\phi_{T,\min}$ is fixed and $\bm u_t$ is Gaussian, tails are sufficiently thin that $\gamma_T$ can be chosen as roughly the order of $\sqrt{\ln(K^2 p)}$ kock2015oracle, leaving room for either almost exponentially large $K$ relative to $T$, or a fairly non-sparse model. On the other hand, if only $m$ moments of $\bm u_t$ exist, $\gamma_t$ should be taken roughly of the order $(K^2 p)^{2/m}$ masini2019regularized, requiring polynomial growth of $K$ compared to $T$ and sparser models.
The most restrictive and crucial assumption needed on the underlying DGP for satisfying Assumption (ref) is the sparsity of the underlying DGP formulated in part ((ref)). The plausibility of this assumption highly depends on the specific application. In many financial applications sparsity (or its approximate version) is natural, for example in portfolio selection when the number of assets is large and the estimation of high-dimensional volatility matrices in financial risk assessment (see fan2011sparse for an overview), as well as in our investigation of Granger causality in networks of realized volatilities in Section (ref). The volatility of one particular stock is likely to have specific channels of contagion rather than affecting the whole stock market at the same time. Shocks to one asset therefore likely propagate through the system via specific channels, which corresponds to sparse lag polynomials. One might worry about systemic shocks affecting many assets; however, the dense covariance matrix $\bm \varSigma_u$ can accommodate simultaneous common shocks. Moreover, the dynamic of such shocks can generally well be captured through a sparse combination of the most important and most affected assets. Similarly, in macroeconomic applications it has been found that a few important variables can capture the effects of unobserved common factors, leading sparse models to perform as well as common factors demol2008forecasting,smeekes2018macroeconomic.
We are now ready to state our main asymptotic result of this section in Theorem (ref) which establishes the asymptotic normality of the post-lasso (generalized) least squares estimator. Here we slightly deviate from the LM test in Algorithm (ref); after the double selection procedure carried out in Step [1], we regress the transformed outcome variables $\widetilde \bm y_t = (\bm G_T \otimes \bm I_T) \bm y_{I}$ on both the Granger causing $\widetilde\bm X_{GC}^\otimes = (\bm G_T \otimes \bm I_T) \bm X_{GC}^\otimes$ and selected variables $\widetilde \bm X_{\hat{S}^\otimes}^\otimes = (\bm G_T \otimes \bm I_T) \bm X_{\hat{S}^{\otimes}}^{\otimes} $
The transformation by the matrix $\bm G_T$ allows for the GLS esitmation needed in the LM procedure by taking $\bm G_T = \hat\bm \varSigma_{u,I}^{-1/2}$, while OLS is performed with $\bm G_T = \bm I_{N_I}$. In the latter case the theorem provides the foundation for the Wald test discussed in Remark (ref) (minus the required variance estimation for that test). We state this result separately as it is interesting in its own right, and can be used to establish validity of other tests such as the Wald test.
Theorem (ref) establishes the asymptotic normality of the post-double-selection OLS estimators. The statement `uniformly in all DGPs that satisfy Assumption (ref)' should be interpreted as the theorem holding uniformly over a parameter space that is defined such that Assumption (ref) holds for all parameters in that parameter space. Importantly, no beta-min conditions on the smallest magnitude of parameters are required, thus alleviating the post-selection inference problem. We refer to Comments 3.4 and 3.5 in belloni2014inference for further details regarding the uniformity. The limit distribution of the LM test now follows straightforwardly from Theorem (ref), and is stated in the corollary below.
Theorem (ref) establishes the limiting distribution of the PDS-LM test under an additional condition on the (co)variances of the partial regression errors, which is satisfied if the errors are iid. To allow for heteroskedaticity the LM test has to be modified, which would only lead to more cumbersome proofs without adding any novelty specific to the high-dimensional case. Therefore we focus on the homoskedastic case here, although we do consider a heteroskedasticity-robust version of the test in the volatility application in Section (ref).\footnote{Note that this is no different for the Wald test, for which the variance estimation has to be adjusted as well.}
We now evaluate the finite-sample performance of our proposed Granger causality test. We consider three Data Generating Processes (DGPs) inspired by kock2015oracle:
The diagonal VAR in DGP1 respects the sparsity assumption while in DGP2 the entries are set to decrease exponentially fast in the distance from the main diagonal and hence the sparsity assumption is not met. DGP2 could be empirically motivated by looking e.g. at financial interconnectedness. Financial institution, such as banks, lend to and borrow from one another becoming interconnected through interbank credit exposures. The financial distress experienced by one bank is likely to be most heavily transmitted the closer the connections are as well as less transmitted, the weaker the connections. DGP3 is a block-diagonal system. Such a structure is motivated by e.g. typical quarterly macroeconomic models capturing business cycle dynamic and monetary and fiscal policy effects. One such example is DSGE models, where the dynamic of the economy through time is monitored on quarterly frequency. Note that as written above, DGP1 satisfies the null of no Granger causality from unit 2 to 1, while DGP2 and DGP3 do not. Therefore, we adapt DGP 1 for the power analysis by setting the coefficient in position $(2,1)$ equal to 0.2. Conversely, we set the same coefficient equal to zero for DGP2 and DGP3 for the size analysis.
We choose our series of interest as $I=\{2\}$ and $J = \{1\}$, therby focusing on the case where we have single variables of interest for both elements of the test. Here we consider for simplicity $p=1$ lag, namely the same lag-length as in the DGPs, so $j=1$. The equation of interest can then be written as
Hence, for each DGP we test $H_0:\beta_{GC}=0$ against $H_1:\beta_{GC}\neq 0$ using our proposed PDS-LM test.
Table (ref) reports the size and power of the test for 1000 replications by using different combinations of time series length $T=(50,100,200,500)$ and number of variables in the system $K=(10,20,50,100)$ and a fixed lag-length $p=1$. All the rejection frequencies are reported using a burn-in period of fifty observations. For each scenario, AIC, BIC and EBIC are compared with the theoretical choice of the tuning parameter $\lambda^{th}$ and time series cross validation $\lambda^{TSCV}$ as described in Subsection (ref).
Simulations are also 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|}$ by using two scenarios of correlation: $\rho=(0,0.7)$. The first case corresponds to no correlation, and is equivalent to set $\bm \varSigma = \bm I_{K}$.
In the Appendix we provide some additional simulation results. First, Table (ref) reports the simulation results for all three DGPs using $\Sigma_{i,j}=0.7^{|i-j|}$. Second, we investigate the Wald version of our test in Table (ref). Third, in Table (ref) we investigate the effects of miss-specification of the lag length by estimating the over-specified VAR$(p+1)$ instead of the true-order VAR$(p)$.\footnote{For both the Wald test and the over-specified VAR$(p+1)$ we report the simulations for $\Sigma_{i,j}=0.7^{|i-j|}$ and DGP1 only. Results for the other DGPs are available upon request.} Fourth, in Table (ref) we report the results for the size of a bivariate Granger causality test for a non-sparse DGP when using a standard Wald ($F$) test. This test is obviously sensitive to omitted variable bias, and our goal is to demonstrate its effect. Finally, although all results reported here use the finite sample correction in Step 3b of the algorithm, we also investigated the differences with Step 3a. We comment on these results in the next subsection. All results not reported in this paper are available from the authors upon request.
Our proposed approach shows a good performance in terms of size and (unadjusted) power for all DGPs considered. Both for the setting of no correlation and high correlation of errors, sizes are in the vicinity of 5% and power is increasing with the sample size $T$.
Only moderate size distortion is visible in large systems for small samples (e.g. $K\ge 50,\; T=50$). As expected, the test procedure works remarkably well for the sparse DGP1 in high dimensions. However, size properties under the non-sparse DGP2 do not deviate much from its sparse counterpart, although for both DGP2 and DGP3 we do observe a slight deterioration of size when the dimension of the system increases.
Interestingly, the three different information criteria show substantially different behavior. EBIC, due to its very stringent nature, tends to perform well only in very large systems, while it is essentially equal to a bivariate Granger causality test in small systems. We have to add though that the good performance of AIC in particular is somewhat inflated by the imposed lower bound on the penalty; unreported simulations show that without the lower bound AIC performs significantly worse, often selecting too many variables rendering the post-OLS estimation infeasible. The one advantage of using EBIC as information criterion to tune $\lambda$ in the $K>>T$ settings when $T$ is small (e.g. $T= 50,100$) is the possibility to avoid the lower bound on the penalty. However, since this comes at a price of more size distortion, we recommend the use of BIC instead, along with the lower bound on the penalty. When comparing the different choices of the tuning parameter we can narrow down the best performing ones (in terms of size and power) to BIC and $\lambda^{th}$. However, in terms of computational time, estimating the tuning parameter using information criteria is considerably faster.
Comparing our test to the bivariate VAR in Table (ref), it is clear that our proposed PDS-LM is very robust to omitted variable bias, unlike the bivariate test, whose size distortions increase with both the sample size and the number of variables, with sizes of 45% observed for the sample sizes we consider in our application in Section (ref). There we will also further elaborate on this difference between our method and the bivariate test. Table (ref) shows that for sample sizes smaller than $T=500$, rarely the power exceeds 90%. However, one must keep in mind that the powers are not size-adjusted, and thus the high reported power of the low-dimensional test is an artefact of the huge size distortions rather than genuine power. It also seems unreasonable to expect that PDS-LM test has vey high power if $T$ is small; we are still considering large systems with many parameters to estimate, and there seems to be no way around this if one desires to test Granger causality in large systems with many (control) variables. In that sense we may fully expect the bivariate test to also have higher size-corrected power; yet with all its disadvantages and sensitivity to omitted variables this is not a good comparison. All in all, we believe our test still has sufficiently adequate power properties to be useful in practice.
The results of robustness to misspecification of the lag length order with $p=2$ instead of $p=1$ are reported in Table (ref) in Appendix (ref). As the size distortions across the range of considered DGPs are only marginally higher for large $K$ and $T$ comparatively small, the test appears to be quite robust to this misspecification. Again, BIC seems to be the best choice for tuning the penalty for all DGPs. Unreported simulations (available upon request) further show that the finite sample adjustment for the test performed in Step 3b of the algorithm is able to substantially reduce size distortions in small samples compared to the asymptotic version of Step 3a.
We also compare it to a four-variate VAR (4-VAR) model
where we also consider the variation of nominal interest rates (3-months treasury bill) $\Delta(TB3MS)$ and variation of inflation $\Delta^2\log(CPI)$. Finally, we compare our method with the factor-augmented VAR (FAVAR) as introduced by bernanke2005measuring. Following mccracken2016fred, we estimate the static factors with PCA from the original (standardized) dataset, excluding the series of GDP and M1. The significant factors are selected by means of the $PC_p=\frac{\log(\min(K,T))}{\min(K,T)}$ criterion of bai2002determining. We obtain a total of eight factors ($F_1,\ldots,F_8$); the top three series most correlated with each factor are reported in the Appendix B. We then test for Granger non causality using the FAVAR model where we stacked the eight extracted factors along with M1 and GDP.
For 2-VAR, 4-VAR and FAVAR models we select the lag-length via AIC or BIC. We exclude EBIC since BIC is already parsimonious enough for selecting lags in low-dimensional VARs.
Table (ref) reports the results from the Granger causality tests. We find that our PDS-LM method provides strong evidence of Granger causality from real money to real output when the lasso is appropriately tuned with (E)BIC or the theoretical plug-in method. When AIC or time series cross-validation is used, too many variables are selected, resulting in a counter-intuitive high $p$-value.\footnote{In order to obtain the $p$-value for the PDS-LM test with AIC we had to adjust the penalty lower bound to $0.5(T)-4$.} Granger causality in the opposite direction is instead always clearly rejected.
The two small VAR models provide conflicting evidence on Granger causality from money to GDP, as well as the other direction. The 2-VAR shows a weak significance in the $M1\rightarrow GDP$ causal relation, yet the 4-VAR is very far from rejecting the null hypothesis, although the difference between the two VAR models might here be attributed to differences in lag length selection. In the opposite direction, the differences between the two small VAR models are even much larger, with the 2-VAR having $p$-values below 0.1 and the 4-VAR above 0.8, regardless of lag length selection.
The factor-augmented VAR appears to be very sensitive to the selected lag length. Adding a single lag from $p=1$ to $p=2$ decreases the $p$-value for Granger causality from M1 to GDP by nearly 70%, while in the other direction an increase is observed from highly significant to decisively not significant. It appears that, even though BIC selects one lag, this is insufficient to capture enough dynamics. We also observed that when manually increasing the lag length to four the $p$-values appear to remain stable, and are qualitatively similar as those of the appropriately tuned PDS-LM test. However, we have not investigated the estimation of the number of factors, which is notoriously difficult, instead going with the established choice of mccracken2016fred. It is likely that uncertainty about the number of factors will further increase the variability of the outcome of the test using the FAVAR.
To wrap up, our results show empirically that by increasing the information set by considering a high-dimensional VAR model, and thus allowing for the potential interaction of many other indicators, one is able to reduce the effect of the omitted variables and thereby retrieve a clearer picture of the causal relations between money and outcome. \fi
We first investigate the volatility transmission in stock return prices using the daily realized variances of 30 US assets. \footnote{We would like to thank Marcelo C. Meideiros for providing us with the high frequency data on stock prices that we have used to construct the realized variances. See Table (ref) for the stocks considered. The R package HDGCvar is available on the GitHub page of the corresponding author (\url{https://github.com/Marga8}).} Both the computational simplicity and the theoretical foundations make realized volatility measures (realized variance, bi-power variation, median realized variance, etc.) very attractive among practitioners and academics for modelling time varying volatilities and monitoring financial risk. We have considered 10-minute realized variances
using $j=1,\ldots, M$ intraday 10 minutes stock prices $P_{j,t}$. We consider 10 minute returns as this is the frequency that minimizes for our sample the microstructure noise (mcaleer2008realized).\footnote{To determine the optimal frequency, we computed realized variances using different frequencies of 1, 5, 10, 15, 30, 65 and 130 minutes, in addition to the estimation using daily returns. The latter estimation has the advantage of being unbiased but the drawback of being very noisy (pooter2008predicting). To find an optimal trade-off between bias and variance martens2004estimating), mean, variances and mean squared errors (MSE) were computed for each estimation frequency in a similar way as pooter2008predicting, and it was found that the frequency of 10 minutes minimizes the MSE.} We investigate the period from March 2008 until February 2017 (2236 trading days).
Given the time series of realized volatilities as defined in ((ref)), we employ a multivariate version of the heterogeneous autoregressive model (VHAR) of corsi2009simple to model their joint behavior (see also cubadda201911). To formally define the VHAR model, we log-transform the series and we stack the logarithmic RV into a vector $y_t$. The VHAR specification is given by the following model:
where $\bm y_{t}^{(week)} = \frac{1}{5} \sum_{j=0}^{4} \bm y_{t-j}$ and $\bm y_{t}^{(month)} = \frac{1}{22} \sum_{j=0}^{21} \bm y_{t-j}$ are the vectors containing the average volatility over the last 5 (week) and 22 (month) days. Granger causality in this context represents contagion, or spillover, of volatility from one asset to another. To test for the null hypothesis of no Granger causality / no volatility spillovers from $y_{k,t}$ to $y_{i,t}$ against the alternative of spillovers, we test
where $\beta_{i,k}^{(1)}$ is the $(i,k)$-th element of $B^{(1)}$. We perform this test for every $(i,k)$-pair to obtain the full $29\times 29$ network of spillover effects. As heteroskedasticity is likely present in these data, we robustify the PDS-LM procedure by implementing the heteroskedasticity-robust LM test such as for example described in wooldridge2015introductory. The full algorithm for the heteroskedasticity-robust PDS-LM test is given in Appendix (ref).\footnote{In the presence of heteroskedasticity, one might prefer the Wald version of the test, as this can be corrected in the standard way by using heteroskedasticty-robust standard errors. Empirically we found hardly any differences between the LM and Wald versions.}
We now report the results of our spillover tests for the volatility network. We use BIC to select the tuning parameter of the lasso, and perform the Granger causality tests with a 1% significance level.\footnote{We do not perform a correction for multiple testing, as this would only qualitatively affect our results. Moreover, our goal is not to identify exactly the set of spillovers, but to get a feeling of the relations between two variables at a time. As such, we believe a multiple testing correction is not needed, though it can be easily implemented.} Figure (ref) reports the transmission networks of volatilities estimated with the high-dimensional HVAR (PDS-LM HVAR), bivariate Granger causality tests (BiHVAR) for each pair of stocks, Granger causality tests from a full-system VAR (FullHVAR). The latter is feasible because of our large time series dimension with $T=2236$. For all methods we consider heteroskedasticity-robust variants.
The drawback of that log transform is that the interpretation of the original series, in our case the volatilities, is lost as the combinations involve nonlinear transforms of both realized variances and covariances. This is not compatible with the aim of this paper.
The second approach uses the log realized volatilities and the correlations separately, as done by for instance oh2016high. The underlying idea, following the DCC model of engle2002dynamic, is to decompose $ RC_{t}^{(d)}=D_{t}^{(d)}R_{t}^{(d)}D_{t}^{(d)}$ with $D_{t}^{(d)}$ a diagonal matrix with the square root of the individual realized variance and $R_{t}^{(d)}$ the realized correlation matrix. oh2016high use the HVAR model structure for each realized volatilities, they consequently assume no Granger causality across volatilities.
We propose something which is, to some extent, in between these two approaches. We look at two separate objects as in the DCC model, but stack the log of the realized variances $\bm y_{1t}^{(d)^{\prime }}$ and $z$-transforms $\bm y_{2t}^{(d)}=\mathrm{arc}\tanh \left( \mathrm{vech}(\bm R_{t}^{(d)})\right) $ of the realized correlations in a larger vector $\bm y_{t}^{(d)}=(\underset{1\times 30}{\underbrace{\bm y_{1t}^{(d)^{\prime }}}}, \underset{1\times 435}{\underbrace{\bm y_{2t}^{(d)^{\prime }}}})^{\prime }$ on which we estimate a VHAR of dimension 465.
In this HVAR each of the 465 equations depends on 1395 dynamic parameters plus the constant. We focus on the 30 equations corresponding to $\bm y_{1t}^{(d)}$ volatilities and consequently the bivariate causalities between these realized volatilities as in the previous section. Figure (ref) reports a total of 113 connections, which is about twice the connections in Figure (ref). In red we highlighted the 31 common connections with Figure (ref). Interestingly, adding more variables therefore allows us to uncover more relations. It seems that this allows us to uncover partial effects that were previously obscured by counteracting effects of the correlations. Importantly, the number of connections is still far less than compared to the BiHVAR in Figure (ref), and the PDS-LM HVAR is still able to deliver a clear picture of the causal connections when the system considered is high-dimensional. While the different connections found here obviously also lead to a different clustering, Figure (ref) shows that the clustering is quite similar, certainly regarding qualitative conclusions.
We propose an LM test in order to test for Granger causality in high-dimensional VAR models. We employ a post-double selection procedure using the lasso to select the set of relevant covariates in the system. The double selection step allows to substantially reduce the omitted variable bias and thereby allowing for valid post-selection inference on the parameters.
We provide an extensive simulation study to evaluate the performance of our method in finite samples, paying particular attention to the tuning of the penalty parameter. We compare different information criteria, time series cross-validation and a plug-in method based on theoretical arguments, and find that generally BIC and the theoretically tuned penalty perform best. However, to use information criteria in systems with a significantly larger number of variables than observations, a lower bound on the penalty parameter is needed to prevent too many variables being selected. The simulations also show that, when properly tuned, our proposed PDS-LM test attains good results both for size and power under different DGPs. Especially, it is shown to be robust both to non-sparse settings as well as to lag-length overspecification.
We also empirically investigate the usefulness of our method in a study where we apply our PDS-LM method to a high-dimensional VHAR process in order to construct a contagion network of volatility spillovers for 30 large capital stocks, also accounting for effects from changing correlations. We find that by increasing the information set through considering a high-dimensional VAR model instead of bivariate models, we are able to obtain more realistic effects than in low-dimensional models. Furthermore, even when the sample size is not large enough to use standard full-system VAR techniques, our method remains reliable and delivers accurate results.
Note that unlike belloni2014high, we do not give a “truly” causal interpretation to the established Granger causalities. In how far Granger causality is a useful concept to study true causality is (and has long been) open to debate, see for example eichler2013causal and the references therein. Moreover, though it appears desirable and in line with granger1969investigating's (granger1969investigating) original intentions to make the information set as large as possible, it is well known in the literature on graphical models eichler2013causal for causality that considering only the full model is not sufficient for establishing true causal relations from Granger causal ones. For instance, one-period Granger causality in systems with more than two variables cannot capture indirect causal chains spanning over multiple periods. However, the analysis of the full model is a necessary ingredient for any study of causality in a graphical framework. It would therefore be an interesting avenue for further research to study how the method proposed here could fit into such a graphical framework.