# Inference in Non-stationary High-Dimensional VARs

Hecq, A.†, Margaritella, L.‡111The second and third author thank NWO for
financial support under grant number 452-17-010. Previous versions of this
paper have been presented at the (EC)2 2020, ICEEE 2021, IAAE 2021 and NESG
2022 conferences as well as the 2022 Maastricht Dimensionality Reduction and
Inference in High-dimensional Time Series and the 2022 Maastricht-York
Economics workshops. We gratefully acknowledge the comments by participants at
these conferences. In addition, we thank Anders B. Kock, Joakim Westerlund,
Gianluca Cubadda, Etienne Wijler for their valuable feedback. All remaining
errors are our own. Address correspondence to: Luca Margaritella, Department
of Economics, Lund University, P.O. Box 7082, S-220 07 Lund, Sweden. E-mail:
<EMAIL_ADDRESS>Smeekes, S.† †Maastricht University, ‡Lund
University

###### Abstract

In this paper we construct an inferential procedure for Granger causality in
high-dimensional non-stationary vector autoregressive (VAR) models. Our method
does not require knowledge of the order of integration of the time series
under consideration. We augment the VAR with at least as many lags as the
suspected maximum order of integration, an approach which has been proven to
be robust against the presence of unit roots in low dimensions. We prove that
we can restrict the augmentation to only the variables of interest for the
testing, thereby making the approach suitable for high dimensions. We combine
this lag augmentation with a post-double-selection procedure in which a set of
initial penalized regressions is performed to select the relevant variables
for both the Granger causing and caused variables. We then establish uniform
asymptotic normality of a second-stage regression involving only the selected
variables. Finite sample simulations show good performance, an application to
investigate the (predictive) causes and effects of economic uncertainty
illustrates the need to allow for unknown orders of integration.

Keywords: Granger causality, Non-stationarity, Post-double-selection, Vector
autoregressive models, High-dimensional inference.
JEL codes: C55, C12, C32.

##  1 Introduction

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 Granger (1969). 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
Sims et al. (1990) and Stock and Watson (1989) and it is still argued to this
day, see e.g., Miao et al. (2023). Financial applications of Granger causality
include, among others, Billio et al. (2012) who study the connectedness among
monthly returns of hedge funds, banks, broker/dealers, and insurance
companies. Hecq et al. (2021) 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, Stern and Kaufmann (2014)
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 (see e.g., Seth et al., 2015, Friston et al., 2013,
for some reviews).

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 (Granger, 1969), conditioning on
a given information set plays a central role.222Throughout 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” (Granger, 1980, p.330).
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, Hecq et al. (2021) (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 Tibshirani (1996), with the post-double selection procedure of
Belloni et al. (2014) 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 (see e.g. Smeekes and Wijler,
2020, and the references therein), 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 (see Dickey and Fuller,
1979). 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 (Smeekes and Wilms, 2023). 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
(Cochrane, 1991) 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 Toda and Yamamoto (1995) who first
considered a simple lag augmentation of the system, which they showed provides
asymptotically normal estimators irrespective of potential unit roots and
cointegration (see also Dolado and Lütkepohl, 1996). The same approach has
also recently been used to develop inference on impulse response functions
(Inoue and Kilian, 2020, Montiel Olea and Plagborg-Møller, 2021) 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 2 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 3. In Section 4 we investigate finite-sample performance through a
simulation study, while Section 5 proposes a data-driven method to find a
sensible upper bound for the lag length in the VAR. Section 6 uses the
proposed testing framework to investigate the causes and effects of economic
uncertainty in the context of the FRED-MD dataset. Section 7 concludes. In
Appendix A 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 B is devoted to the
validation of high-level assumptions, while additional empirical results are
presented in Appendix C.

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}|x_{i}|^{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$. $|S|$ denotes the cardinality of the set
$S$. We use $\xrightarrow{p}$ and $\xrightarrow{d}$ to denote convergence in
probability and distribution, respectively.

##  2 Granger Causality Tests for Nonstationary High-Dimensional Data

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 2.1 first presents the model and
lag-augmented Granger causality test. Next, 2.2 sets out the post-double-
selection procedure.

###  2.1 Lag-Augmented Granger Causality Testing

Let $\bm{z}_{1},\ldots,\bm{z}_{T}$ be a $K$-dimensional multiple time series
process, where $\bm{z}_{t}=(y_{t},x_{t},\bm{w}_{t})^{{}^{\prime}}$. Here
$y_{t}$ is the series we would like to test being Granger caused, $x_{t}$ the
potentially Granger causing series and $\bm{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 $\bm{z}_{t}$ is generated by a VAR($p$) process as

$\bm{z}_{t}=\bm{A}_{1}\bm{z}_{t-1}+\cdots+\bm{A}_{p}\bm{z}_{t-{p}}+\bm{u}_{t},\qquad
t=p+1,\ldots,T\;,$ (1)

where $\bm{A}_{1},\ldots,\bm{A}_{p}$ are $K\times K$ parameter matrices and
$\bm{u}_{t}$ is a $K\times 1$ vector of error terms.

######  Assumption 1.

The VAR model in (1) satisfies:

  1. 1.

$\\{\bm{u}_{t}\\}_{t=1}^{T}$ is an mds with respect to the filtration
$\mathcal{F}_{t}=\sigma(\bm{z}_{t},\bm{z}_{t-1},\bm{z}_{t-2},\ldots)$
$\bm{u}_{t}$ such that $\mathbb{E}(\bm{u}_{t}|\mathcal{F}_{t-1})=\bm{0}$ for
all $t$; the $K\times K$ covariance matrix
$\bm{\varSigma}_{u}=\mathbb{E}(\bm{u}_{t}\bm{u}_{t}^{\prime})$ is positive
definite and $\mathbb{E}|\bm{u}_{t}|^{2+\delta}\leq\infty$, for $\delta>0$.

  2. 2.

The roots of $\det(\bm{I}_{K}-\sum_{j=1}^{p}\bm{A}_{j}z^{j})$ can either lie
on the unit disc or outside.

Note that Assumption 1(2) allows for the time series to have unit roots and be
cointegrated. We specifically allow elements of $\bm{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 3 and Appendix B.

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 $\bm{w}_{t}$. For the moment we
assume the lag-length $p$ in (1) to be known; we shall further elaborate on
data-driven ways to select $p$ in Section 5. 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
$\bm{X}_{-p}=(\bm{x}_{-1},\ldots,\bm{x}_{-p})$ denote the $T\times p$ matrix
containing the $p$ lags of the Granger causing variable $x_{t}$, where
$\bm{x}_{-j}$ is the vector containing the observations corresponding to its
$j$-th lag.333To 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 $\bm{Y}_{-p}$
containing the $p$ lags of the Granger caused variable $y_{t}$ and
$\bm{W}_{-p}$ for the conditioning set such that $\bm{W}_{-p}$ is a
$T\times(K-2)p$ matrix. Then, our equation of interest for the Granger
causality testing can be stated as

$\bm{y}=\bm{X}_{-p}\bm{\beta}+\bm{Y}_{-p}\bm{\delta}_{1}+\bm{W}_{-p}\bm{\delta}_{2}+\bm{u}=\bm{X}_{-p}\bm{\beta}+\bm{V}\bm{\delta}+\bm{u},$
(2)

where for notational simplicity we define $\bm{V}=(\bm{Y}_{-p},\bm{W}_{-p})$
as the matrix containing all control variables and
$\bm{\delta}=(\bm{\delta}_{1}^{\prime},\bm{\delta}_{2}^{\prime})^{\prime}$ its
coefficients.

Testing for no Granger causality is then equivalent to testing the following
null hypothesis:

$H_{0}:\bm{\beta}=\bm{0}\quad\text{against}\quad H_{1}:\bm{\beta}\neq\bm{0}.$
(3)

To account for the potential unit roots in the system, we follow the approach
pioneered by Toda and Yamamoto (1995) and Dolado and Lütkepohl (1996) 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

$\bm{y}=\bm{X}_{-p}\bm{\beta}+\sum_{j=1}^{d}\beta_{p+j}\bm{x}_{-(p+j)}+\bm{V}\bm{\delta}+\bm{u}=\bm{X}_{-(p+d)}\bm{\beta}_{+}+\bm{V}\bm{\delta}+\bm{u},$
(4)

where $\bm{X}_{-(p+d)}=(\bm{X}_{-p},\bm{x}_{-(p+1)},\ldots,\bm{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
$\bm{\beta}_{+}=(\bm{\beta}^{\prime},\bm{0}_{d}^{\prime})^{\prime}$, as we are
adding redundant variables. It should be clear that by testing whether the
first $p$ elements of $\bm{\beta}_{+}$ – those corresponding to $\bm{\beta}$
in (1) – 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)$:

$y_{t}=\beta x_{t-1}+u_{1,t},\qquad x_{t}=\rho
x_{t-1}+u_{2,t},\quad{\left\lvert\rho\right\rvert}\leq 1.$ (5)

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 (5) as

$y_{t}=\beta x_{t-1}+0x_{t-2}+u_{1,t}=\beta u_{2,t-1}+\beta\rho
x_{t-2}+u_{1,t},$ (6)

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
A.1 we show formally that not only the equivalence in (6) 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 Toda and Yamamoto
(1995) 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.

######  Remark 1.

For ease of exposition we confine our attention to the study of bivariate
Granger causality relations, conditional on a large information set. At the
cost of more involved notation and algebra, our approach can be extended to
Granger causality between multiple variables. In that case, lag augmentation
is needed for all Granger causing variables, which is only feasible if the
block of Granger causing variables is not too large. HMS provide details about
this in the stationary setting; the same approach can be adapted here.

######  Remark 2.

One important aspect of the current framework is that the lag-length $p$ of
the VAR is necessarily larger than, or at least equal to, the suspected
maximum order of integration $d$. It is therefore important to specify $p$
correctly in practice. We return to this issue in Section 5. One might also
worry that specifying $d$ too high, or having mixed orders of integration
among multiple Granger causing variables would lead to over-differencing,
i.e., moving average unit roots being introduced by differencing stationary
time series (see e.g., Chang and Dickey, 1994). This, however, does not happen
here as the additional lags of the Granger causing and Granger caused
variables are used “at convenience”; if they are not needed because the
variables are already stationary, inclusion of these redundant will at most
marginally decrease the power of the test given the small over-specification
of the lag length.

######  Remark 3.

Even if the variables in a particular dataset are thought to be at most of
order $I(1)$, it may still pay off to take $d=2$. The choice of $d=2$ is
supported by simulations reported later in Section 4. It is well known that
when one or more roots of the characteristic polynomial are close to unity,
the distribution of the least squares estimator becomes skewed, yielding
estimators which tend to underestimate the true autoregressive parameters (see
e.g. Fuller, 2009). This also causes difficulties in performing inference on
these parameters. By augmenting with $d=2$ lags, one can avoid any issues with
near unit roots. The simulations reported in Section 4 confirm that in the
presence of substantial autocorrelation beyond the first lag, augmenting with
$d=2$ lags is an easy way to improve the finite sample behaviour of the test.

######  Remark 4.

If the $I(d)$ series object of the Granger causality test are also actually
cointegrated, then the lag-augmentation does not serve any purpose as no
spurious relation would be estimated (Dolado and Lütkepohl, 1996). This causes
a slight loss in power which however is in practice minimal, as shown in
Section 4. If, however, the series in object are not cointegrated, then the
lag-augmentation becomes paramount to render the series difference-stationary
before the estimation and testing. Note that also all sorts of situations in
between are allowed: some variables may be cointegrated, others contain ‘pure’
stochastic trends, and others may be stationary. In fact, in our high-
dimensional setup where we apply the test to a large dataset containing many
variables, such mixed properties seem likely.

###  2.2 Inference after selection by the lasso

Assuming that $\bm{\delta}$ is sparse, in the sense that many elements of the
vector are equal to zero, one might consider estimating (4) 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 (Leeb and Pötscher,
2005), 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 Belloni et al.
(2014) 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 $\bm{X}_{-p,\backslash\\{j\\}}$ denote the matrix $\bm{X}_{-p}$ from which
the $j$-th column, corresponding to the $j$-th lag of $x_{t}$, has been
removed. Let $\bm{Z}_{-p}=(\bm{X}_{-p},\bm{Y}_{-p},\bm{W}_{-p})$ denote the
matrix of all (non-augmented) variables and
$\bm{Z}_{-p,\backslash\\{j\\}}=(\bm{X}_{-p,\backslash\\{j\\}},\bm{Y}_{-p},\bm{W}_{-p})$.
Then we consider the following regressions in the first step:

$\bm{y}=\bm{Z}_{-p}\bm{\eta}_{0}+\bm{u},\qquad\bm{x}_{-j}=\bm{Z}_{-p,\backslash\\{j\\}}\bm{\eta}_{j}+\bm{e}_{j},\quad
j=1,\ldots,p.$ (7)

In contrast to HMS, who follow Belloni et al. (2014) by only regressing
$\bm{y}$ and $\bm{x}_{-j}$ on the controls $\bm{Y}_{-p}$ and $\bm{W}_{-p}$, we
also add all other lags of $x_{t}$. This is done to avoid that the error terms
in (7) contain unit roots and the regressions become spurious.

######  Remark 5.

Note that in the situation where $p=d$, there may still be a risk of spurious
regression in the first step. E.g., consider $p=d=1$ where the regression of
$x_{t-1}$ on the remaining variables will not contain any lags (or leads) of
$x_{t-1}$. This yields a spurious regression if $x_{t}$ is a non-cointegrated
unit root process. We therefore recommend to always take $p\geq d+1$ in such
cases. If this is not desirable, e.g. if the sample size is too small to
sustain such a large $p$, it is also possible to augment the first-stage
regressions with an additional lag of $x_{t}$ to eliminate the possibility of
spurious regression. In the following we work under the condition that the
error terms $\bm{e}_{j}$ are $I(0)$, thereby implicitly assuming that one of
the two measures suggested above have been taken if needed.

Although the concept is tricky with integrated variables, the coefficients
$\bm{\eta}_{j}$ in (7) can be thought of as best linear predictors. Indeed, as
the error terms $\bm{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 2, we assume that the coefficients
$\bm{\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 (7).444We 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 $\bm{y}$ on
$\bm{X}_{-p}$, as they either have non-zero coefficients in (2) _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 $\bm{V}_{\mathcal{S}}$ denote the matrix consisting of only those columns
of $\bm{V}$ that corresponds to the selected variables in $\mathcal{S}$, and
$\bm{\delta}_{\mathcal{S}}$ the corresponding parameter vector. Our goal of
the PDS regression is then to recover the regression

$\bm{y}=\bm{X}_{-(p+d)}\bm{\beta}_{+}+\bm{V}_{\mathcal{S}}\bm{\delta}_{\mathcal{S}}+\bm{u},$
(8)

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 (see
Tibshirani, 1996) simultaneously performs variable selection and estimation of
the parameters in (7) by solving the following minimization problems

$\begin{split}&\hat{\bm{\eta}}_{0}=\underset{\bm{\eta}}{\arg\min}\bigg{(}T^{-1}{\left\lVert\bm{y}-\bm{Z}_{-p}\bm{\eta}\right\rVert}_{2}^{2}+\lambda{\left\lVert\bm{\eta}\right\rVert}_{1}\bigg{)},\\\
&\hat{\bm{\eta}}_{j}=\underset{\bm{\eta}}{\arg\min}\bigg{(}T^{-1}{\left\lVert\bm{x}_{-j}-\bm{Z}_{-p,\backslash\\{j\\}}\bm{\eta}\right\rVert}_{2}^{2}+\lambda{\left\lVert\bm{\eta}\right\rVert}_{1}\bigg{)},\quad
j=1,\ldots,p,\end{split}$ (9)

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 7 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 (7) using a penalized regression technique such as in (LABEL:eq:lasso). 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

$\bm{y}=\bm{X}_{-(p+d)}\bm{\beta}_{+}+\bm{V}_{\hat{\mathcal{S}}}\bm{\delta}_{\hat{\mathcal{S}}}+\bm{u},$
(10)

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 1 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.

Algorithm 1 Post-double selection lag-augmented (PDS-LA) Granger causality
test

  * [1]

Estimate the initial partial regressions in (7) by an appropriate sparsity-
inducing estimator such as the (adaptive) lasso, and let
$\hat{\bm{\eta}}_{0},\ldots,\hat{\bm{\eta}}_{p}$ denote the resulting
estimators. Collect the selected variables in the sets
$\hat{\mathcal{S}}_{j}=\\{m>p:\hat{\eta}_{m,j}\neq 0\\}$ for $j=0,1,\ldots,p$.

  * [2]

Let $\hat{\mathcal{S}}=\bigcup_{j=0}^{p}\hat{\mathcal{S}}_{j}$ denote the set
of all variables selected in at least one regression at Step [1]. Denote its
cardinality by $\hat{\mathcal{S}}={\left\lvert\hat{\mathcal{S}}\right\rvert}$.

Then, consider one of the following two tests:

PDS-LA-LM

  * [3]

Regress $\bm{y}$ on
$\bm{V}_{+,\hat{\mathcal{S}}}:=(\bm{x}_{-(p+1)},\ldots,\bm{x}_{-(p+d)},\bm{V}_{\hat{\mathcal{S}}})$
by OLS and obtain the residuals $\hat{\bm{\xi}}$.

  * [4]

Regress $\hat{\bm{\xi}}$ onto $\bm{V}_{+,\hat{\mathcal{S}}}$ plus the Granger
causing variables $\bm{X}_{-p}$, retaining the residuals $\hat{\bm{\nu}}$.
Then obtain
$R^{2}=1-\hat{\bm{\nu}}^{\prime}\hat{\bm{\nu}}/\hat{\bm{\xi}}^{\prime}\hat{\bm{\xi}}$

  * [5a]

Reject $H_{0}$ if $\text{LM}=TR^{2}>q_{\chi_{p}^{2}}(1-\alpha)$, where
$q_{\chi_{p}^{2}}(1-\alpha)$ is the $1-\alpha$ quantile of the $\chi^{2}$
distribution with $p$ degrees of freedom.

  * [5b]

Reject $H_{0}$ if
$\bigg{(}\frac{T-\hat{\mathcal{S}}-(p+d)}{p}\bigg{)}\bigg{(}\frac{R^{2}}{1-R^{2}}\bigg{)}>q_{F_{p,T-\hat{\mathcal{S}}-(p+d)}}(1-\alpha)$,
where $q_{F_{p,T-\hat{\mathcal{S}}-(p+d)}}$ is the $1-\alpha$ quantile of the
$F$ distribution with $p$ and $T-\hat{\mathcal{S}}-(p+d)$ degrees of freedom.

PDS-LA-Wald

  * [3]

Regress $\bm{y}$ on $\bm{X}_{-p}$ and
$\bm{V}_{+,\hat{\mathcal{S}}}:=(\bm{x}_{-(p+1)},\ldots,\bm{x}_{-(p+1)},\bm{V}_{\hat{\mathcal{S}}})$
by OLS, obtaining $\hat{\bm{\beta}}_{p}$, the estimated coefficients
corresponding to $\bm{X}_{-p}$, and the residuals $\hat{\bm{u}}_{+}$.

  * [4]

Reject $H_{0}$ if
$\text{Wald}=\hat{\bm{\beta}}_{p}^{\prime}\left(\widehat{\mathbb{V}ar}\hat{\bm{\beta}}_{p}\right)^{-1}\hat{\bm{\beta}}_{p}>q_{\chi_{p}^{2}}(1-\alpha)$,
where $\widehat{\mathbb{V}ar}\hat{\bm{\beta}}_{p}$ is the OLS standard error
of $\hat{\bm{\beta}}_{p}$.

######  Remark 6.

In Algorithm 1, the choice among Step [5a] or [5b] does not affect the finite
sample results of the test whenever the sample size $T$ is large enough. The
small sample correction in [5b] (see Kiviet, 1986) has a wider practical
applicability since [5a] suffers from size distortion in small samples,
therefore in Section 4 we always use [5b] for the Monte-Carlo simulations of
the PDS-LA-LM test. Unreported simulations show that the Wald test performs
very similarly to the LM test, if with slightly bigger size distortions.

######  Remark 7.

The algorithm designed in HMS employs a lower bound on the penalty to ensure
that in each selection regression at most $cT$ terms gets selected, for some
$0<c<1$. Similarly, we also employ a $c=0.5$ lower bound on the selected
variables in Algorithm 1. This ensures that in each equation the lasso does
not select too many variables, as this would render the union too large and
hence infeasible for post-least-squares estimation. Our PDS procedure does not
require consistent model selection but only consistency (see Assumption 2(3)
below). Mistakes are allowed to occur in the selection: variables might be
incorrectly included as long as the estimator remains sufficiently sparse and
consistency is guaranteed. Note that even with the lower bound it remains
possible that the lasso selects sufficiently distinct variables at every
selection step. Should such a case occur where the number of selected
variables $\hat{\mathcal{S}}$ is larger than the sample size, post-OLS would
be infeasible. One way to avoid this would be to impose an ad-hoc increase on
the tightness of the bound on the selected variables. Alternatively one could
switch to a sparser estimator, such as the adaptive lasso.

######  Remark 8.

Although it is not necessary for the theory to hold, we advocate to always
include the $p$ lags of $y_{t}$ in the second stage regression. Erroneously
omitting them might induce spurious regression, which is to be avoided. We
achieved this by including them without penalty in the first stage
regressions, such that they will be included in the active set by default.
Additionally, we do not penalise the lags of $x_{t}$ in the first-stage
regressions, to further reduce the probability of spurious regression and
improving finite sample behaviour.

##  3 Theoretical Properties

In this section we present the main theoretical result. We show the post-
selection, lag-augmented, least squares estimator $\hat{\bm{\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
$\bm{\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
$\bm{Z}_{-p}^{\prime}\bm{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 $\bm{z}_{t}$ can be written as the cointegrated
system

$\Delta\bm{z}_{t}=\bm{A}\bm{B}^{\prime}\bm{z}_{t-1}+\bm{u}_{t},$

where $\bm{u}_{t}$ is $I(0)$ and $\bm{A}$ and $\bm{B}$ are $n\times r$
matrices with $r<n$. Then define $\bm{\zeta}_{t}:=\bm{z}_{t}\bm{Q}$ where
$\bm{Q}:=(\bm{B},\bm{A}_{\perp})^{\prime}$ and $\bm{A}_{\perp}$ is the
$n\times n-r$ orthognal complement of $\bm{A}$, such that
$\bm{A}_{\perp}^{\prime}\bm{A}=\bm{0}$. Then we can write

$\Delta\bm{\zeta}_{t}=\bm{Q}\bm{A}\bm{B}^{\prime}\bm{Q}^{-1}\bm{\zeta}_{t-1}+\bm{Q}\bm{u}_{t}=\begin{bmatrix}\bm{B}^{\prime}\bm{A}&\bm{0}\\\
\bm{0}&\bm{0}\end{bmatrix}\bm{\zeta}_{t-1}+\bm{Q}\bm{u}_{t}.$

Now the first $r$ variables in $\bm{\zeta}_{t}$ are $I(0)$, while the
remaining ones are $I(1)$. Further details on constructing such a matrix
$\bm{Q}$ can be found in, e.g., Lütkepohl (2005, Chapter 7).

Given the existence of such a matrix $\bm{Q}$, we can without loss of
generality assume that we may partition $\bm{\zeta}_{t}$ as
$\bm{\zeta}_{t}=(\bm{\zeta}_{0,t},\bm{\zeta}_{1,t},\bm{\zeta}_{2,t})^{\prime}$,
where $\bm{\zeta}_{i,t}$ contains the $I(i)$ variables. We then consider a
scaling matrix

$\bm{D}_{T}:=\begin{bmatrix}T^{1/2}\bm{I}_{n_{0}}&\bm{0}&\bm{0}\\\
\bm{0}&T^{-1}\bm{I}_{n_{1}}&\bm{0}\\\
\bm{0}&\bm{0}&T^{-2}\bm{I}_{n_{2}}\end{bmatrix},$ (11)

where $n_{i}$ corresponds to the number of $I(i)$ variables in
$\bm{\zeta}_{t}$. Finally, we let $\bm{G}_{T}:=\bm{Q}\bm{D}_{T}$ as the final
‘scaling and rotation’ matrix.

######  Assumption 2.

Let $\delta_{T}$ and $\Delta_{T}$ denote sequences such that
$\delta_{T},\Delta_{T}\rightarrow 0$ as $T\rightarrow\infty$. For $\bm{D}_{T}$
as defined in (11), assume there exists a $Kp\times Kp$ matrix
$\bm{G}_{T,\bm{Z}}=\bm{Q}\bm{D}_{T}$ conformably with $\bm{Z}_{-p}$ such that
the following holds:

  1. 1.

Deviation Bound:  With probability at least $1-\Delta_{T}$ we have that the
empirical process satisfes
${\left\lVert\bm{G}_{T}^{-1\prime}\bm{Z}_{-p}^{\prime}\bm{u}\right\rVert}_{\infty}\leq\bar{\gamma}_{T}$,
for the deterministic sequence $\bar{\gamma}_{T}$.

  2. 2.

Boundedness: Let $\bm{\beta}$ in (2) be in the interior of a compact parameter
space $\mathbb{B}\subset\mathbb{R}^{K}$.

  3. 3.

Consistency: With probability $1-\Delta_{T}$, the following error bounds for
$\bm{\eta}_{j}$ and $\hat{\bm{\eta}}_{j}$, defined in (7) and
(LABEL:eq:lasso), hold:

$\begin{split}&{\left\lVert\bm{Z}_{-p}\left(\hat{\bm{\eta}}_{0}-\bm{\eta}_{0}\right)\right\rVert}_{2}^{2}\leq\delta_{T}^{2}T^{1/2},\qquad{\left\lVert\hat{\bm{\eta}}_{0}-\bm{\eta}_{0}\right\rVert}_{\infty}\leq\delta_{T}T^{-1/4},\\\
&{\left\lVert\bm{Z}_{-p,\backslash\\{j\\}}\left(\hat{\bm{\eta}}_{j}-\bm{\eta}_{j}\right)\right\rVert}_{2}^{2}\leq\delta_{T}^{2}T^{1/2},\quad
j=1,\ldots,p.\end{split}$

  4. 4.

Sparsity: Let $s={\left\lvert\mathcal{S}\right\rvert}$ and
$\hat{\mathcal{S}}={\left\lvert\hat{\mathcal{S}}\right\rvert}$ denote the
number of active variables in population and sample, respectively. Then with
probability at least $1-\Delta_{T}$, we have that
$\max(s,\hat{\mathcal{S}})\leq\bar{s}_{T}$ for the deterministic sequence
$\bar{s}_{T}$.

  5. 5.

Restricted Sparse Eigenvalues: For any $\bm{\eta}\in\mathbb{R}^{Kp}$ with
${\left\lVert\bm{\eta}\right\rVert}_{0}\leq\bar{s}_{T}$, we have with
probability at least $1-\Delta_{T}$

${\left\lVert\bm{\eta}\right\rVert}_{1}\leq\sqrt{\bar{s}_{T}}{\left\lVert\bm{Z}_{-p}\bm{G}_{T,\bm{Z}}^{-1}\bm{\eta}\right\rVert}_{2}/\kappa_{T,\min},$

for the deterministic sequence $\kappa_{T,\min}>0$.

  6. 6.

Rate Conditions: the deterministic sequences bounding sparsity
($\bar{s}_{T}$), thickness of empirical process tails ($\bar{\gamma}_{T}$) and
minimum eigenvalue ($\kappa_{T,\min}$) must satisfy

$T\frac{\bar{s}_{T}\bar{\gamma}_{T}}{\kappa_{T,\min}}\leq\delta_{T}.$ (12)

Condition (1) bounds the empirical process with high probability. Such
deviation bounds for non-stationary time series can be found in, among others,
Smeekes and Wijler (2021, Lemma A.3), Mei and Shi (2022, Proposition 1) and
Wijler (2022, Lemma 2). Condition (2) is standard and assumes compactness of
the parameter space of the vector $\bm{\beta}$ which in turn implies the
boundedness.

Condition (3) 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 (1), (4) and (5). For stationary time series many
results exist under a variety of settings; see e.g. Masini et al. (2022,
Theorem 1) for VAR models or Adamek et al. (2023, Corollary 1) for general
time series models. Nonstationary time series are treated in Smeekes and
Wijler (2021, Theorem 2), Mei and Shi (2022, Theorem 1 and 2) and Wijler
(2022, Theorem 2).

Condition (4) 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 Belloni et al.
(2014), where it is assumed that the exact sparse model is a (good)
approximation to the true DGP, or weak sparsity following Adamek et al.
(2023), where many non-zero but small coefficients are allowed.

Condition (5) 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., Adamek et al., 2023, Lemma
A.3 and A.5 and Masini et al., 2022, Lemma 3 and Proposition 2), time series
with unit roots require more care. Indeed, Smeekes and Wijler (2021, Theorem
B.2), Mei and Shi (2022, Lemma 2) and Wijler (2022, Theorem 1) show for unit
root regressors that the standard rate of $T^{-2}$ coming from the
$\bm{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 (6) 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.555For an illustration of the complexities of these relations, we
refer to Figure C.1 in Adamek et al. (2023) 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 1.

Let $\hat{\bm{\beta}}_{p}$ denote the OLS estimator of the first $p$ elements
of $\bm{\beta}_{+}$ (equal to $\bm{\beta}$) in (10), and
$\tilde{\bm{\beta}}_{p}$ the corresponding estimator in (8). Then uniformly
over a parameter space $\mathcal{B}$ for which Assumptions 1 and 2 holds for
all elements in $\mathcal{B}$, we have that

$\sqrt{T}(\hat{\bm{\beta}}_{+}-\bm{\beta})=\sqrt{T}(\tilde{\bm{\beta}}_{+}-\bm{\beta})+o_{p}(1).$

Theorem 1 establishes the asymptotic equivalence of the estimator in the
feasible second-stage regression (10) based on the estimated active set
$\hat{\mathcal{S}}$, and the estimator in the infeasible second-stage
regression (10) 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 2(4) – this assumption follows directly from well-known
results in the unit root and cointegration literature; see Appendix B. We here
state the assumption more generally to also accommodate sparsity that
increases with the sample size.

######  Assumption 3.

Let $\\{\bm{\zeta}_{t}\\}_{t=1}^{T}$ denote an $n$-dimensional process, where
$n=n_{T}$ may increase with $T$. Let $n_{0}$, $n_{1}$ and $n_{2}$ denote
integers such that $n=n_{0}+n_{1}+n_{2}$ and partition
$\bm{\zeta}_{t}=(\bm{\zeta}_{0,t},\bm{\zeta}_{1,t},\bm{\zeta}_{2,t})^{\prime}$
comformably with $(n_{0},n_{1},n_{2})$, and assume that
$\mathbb{E}[\bm{\zeta}_{t}u_{t}]=\bm{0}$ for all $t=1,\ldots,T$. Define

$\bm{D}_{T}:=\begin{bmatrix}T^{1/2}\bm{I}_{n_{0}}&\bm{0}&\bm{0}\\\
\bm{0}&T\bm{I}_{n_{1}}&\bm{0}\\\
\bm{0}&\bm{0}&T^{2}\bm{I}_{n_{2}}\end{bmatrix}=:\begin{bmatrix}\bm{D}_{T,0}&\bm{0}&\bm{0}\\\
\bm{0}&\bm{D}_{T,1}&\bm{0}\\\ \bm{0}&\bm{0}&\bm{D}_{T,2}\end{bmatrix},$

and let $\Delta_{T},\delta_{T}$ denote deterministic sequences such that
$\Delta_{T},\delta_{T}\to 0$ as $T\to\infty$. For varying subscripts ‘$\cdot$’
given below, let $\phi_{T,\cdot}$, $\kappa_{T,\cdot}$ and $\gamma_{T,\cdot}$
denote sequences of positive real numbers, such that with probability
$1-\Delta_{T}$ the following statements hold jointly:

  1. 1.

${\left\lVert\bm{D}_{T,i}^{-1}\sum_{t=1}^{T}\bm{\zeta}_{i,t}\bm{\zeta}_{j,t}^{\prime}\bm{D}_{T,j}^{-1}\right\rVert}_{2}\leq\phi_{T,ij}$
for $i,j=0,1,2$ and $i<j$;

  2. 2.

$\lambda_{\min}\left(\bm{D}_{T,0}^{-1}\sum_{t=1}^{T}\bm{\zeta}_{0,t}\bm{\zeta}_{0,t}^{\prime}\bm{D}_{T,0}^{-1}\right)\geq\kappa_{T,0}$
and
$\lambda_{\min}\left(\bm{D}_{T,-0}^{-1}\sum_{t=1}^{T}\bm{\zeta}_{-0,t}\bm{\zeta}_{-0,t}^{\prime}\bm{D}_{T,-0}^{-1}\right)\geq\kappa_{T,12}$,
where
$\bm{\zeta}_{-0,t}=(\bm{\zeta}_{1,t}^{\prime},\bm{\zeta}_{2,t}^{\prime})^{\prime}$
and $\bm{D}_{t,-0}=\text{diag}(\bm{D}_{T,1},\bm{D}_{T,2})$;

  3. 3.

${\left\lVert\bm{D}_{T,i}^{-1}\sum_{t=1}^{T}\bm{\zeta}_{i,t}u_{t}\right\rVert}_{2}\leq\gamma_{T,u,i}$
for $i=0,1,2$;

  4. 4.

${\left\lVert
T^{-1}\sum_{t=1}^{T}\left[\bm{\zeta}_{0,t}\bm{\zeta}_{0,t}^{\prime}-\mathbb{E}\left(\bm{\zeta}_{0,t}\bm{\zeta}_{0,t}^{\prime}\right)\right]\right\rVert}_{2}\leq\delta_{T}$.

Assume that the rates above are bounded as

  1. 1.

$\kappa_{T,12}-(\phi_{T,01}^{2}+\phi_{T,02}^{2})/\kappa_{T,0}\geq\mu_{T,(i)}$;

  2. 2.

$\phi_{T,01}+\phi_{T,02}+2(\phi_{T,01}^{2}+\phi_{T,02}^{2})/\kappa_{T,0}\leq\mu_{T,(ii)}$;

  3. 3.

$\gamma_{T,u,1}+\gamma_{T,u,2}+(\phi_{T,01}+\phi_{T,02})\gamma_{T,u,0}/\kappa_{T,0}\leq\mu_{T,(iii)}$;

for sequences $\mu_{T,(i)},\mu_{T,(ii)},\mu_{T,(iii)}$ that satisfy the
condition

$\mu_{T,(ii)}(\mu_{T,(ii)}+\mu_{T,(iii)})\leq\delta_{T}\mu_{T,(i)}.$ (13)

In addition, let $\bm{R}_{m}$ denote a deterministic $m\times n_{0}$ matrix
where $m<\infty$ is not depending on $T$. Then we have that

$T^{-1/2}\sum_{t=1}^{T}\bm{R}_{m}\bm{\zeta}_{0,t}u_{t}\xrightarrow{d}N(\bm{0},\bm{\varOmega}),$

where
$\bm{\varOmega}=\lim_{T\rightarrow\infty}T^{-1}\bm{R}_{m}\mathbb{E}(\bm{\zeta}_{0,t}u_{t}u_{t}^{\prime}\bm{\zeta}_{0,t}^{\prime})\bm{R}_{m}^{\prime}=\sigma_{u}^{2}\bm{R}_{m}\mathbb{E}(\bm{\zeta}_{0,t}\bm{\zeta}_{0,t}^{\prime})\bm{R}_{m}^{\prime}$.

Assumption 3 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 (13) – 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.
Hamilton (1994) for an overview. We show in Appendix B how Assumption 3 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.,
Zhang et al. (2019, Remark 3.4) and Smeekes and Wijler (2021, Theorem B.3).
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
2(1) 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.

######  Theorem 2.

Let Wald and LM be as defined in Algorithm 1. Assume that Assumption 3 holds
for $\bm{\zeta}_{t}=\bm{Q}_{\mathcal{S}}\bm{v}_{+,\mathcal{S},t}$, where
$\bm{v}_{+,\mathcal{S},t}=(x_{t-p-1},\ldots,x_{t-p-d},\bm{v}_{\mathcal{S},t})$.
Then, uniformly over a parameter space $\mathcal{B}$ on which Assumptions 1
and 2 holds for all elements in $\mathcal{B}$, we have that

$\displaystyle\text{LM},\text{Wald}\xrightarrow{d}\chi_{p}^{2},\qquad\text{as
}T\rightarrow\infty,$

under the null hypothesis that $\bm{\beta}=\bm{0}$.

######  Remark 9.

Though we focus on homoskedastic error terms for simplicity,
heteroskedasticity-robust versions of the test can easily be constructed and
shown to be valid, as this would only affect the low-dimensional part of our
results in Theorem 2. Standard techniques can therefore be used for
constructing heteroskedasticity-robust tests: for the Wald test, the OLS
standard errors can be replaced by Eicker-White standard errors, while the LM
test can be modified as in Wooldridge (1987). We refer to HMS, Algorithm 2 for
a full treatment.

##  4 Monte-Carlo Simulations

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 Kock and Callot (2015):

$\displaystyle\text{DGP1:}\qquad\Delta\bm{z}_{t}=\begin{bmatrix}0.5&0&\ldots&0\\\
0&0.5&\ldots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\ldots&0.5\end{bmatrix}\Delta\bm{z}_{t-1}+\bm{u}_{t},$
$\displaystyle\text{DGP2:}\qquad\Delta\bm{z}_{t}=\begin{bmatrix}(-1)^{|i-j|}a^{|i-j|+1}&\ldots&(-1)^{|i-j|}a^{|i-j|+1}\\\
(-1)^{|i-j|}a^{|i-j|+1}&\ldots&(-1)^{|i-j|}a^{|i-j|+1}\\\
\vdots&\ddots&\vdots\\\
(-1)^{|i-j|}a^{|i-j|+1}&\ldots&(-1)^{|i-j|}a^{|i-j|+1}\end{bmatrix}\Delta\bm{z}_{t-1}+\bm{u}_{t},$

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 4 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.

######  Remark 10.

As mentioned in Remark 7, in order to obtain the results for the size and
power when $T\leq Kp$ we need to impose a lower bound on the lasso penalty
$\lambda$ which guarantees to select at most $c\,T$ variables in each relevant
equation of the VAR, for some $0<c<1$. The bound should be set as strict as
the system requires and often there is not a universal constant $c$ that works
in all settings, therefore this choice needs to be adaptive. For instance, if
the lag length is $p=2$, this implies 3 selection steps plus $d$ augmented
lags of Granger causing variables to be added. In some cases this might lead
to too many variables being selected in order to perform least squares in the
second step. In these cases we tighten the bound using either $c=0.33$ or
$c=0.25$.

Simulation results for the PDS-LA-LM Granger causality test $T$ 50 100 200 500
1000 50 100 200 500 1000 DGP $\rho$ $K$ Size Power 1 0 10 10.4 7.2 5.5 4.2 5.5
23.3 41.0 74.2 99.4 100 20 12.5 9.5 7.3 6.3 5.9 19.7 40.1 74.6 99.5 100 50
10.7 7.1 7.2 5.4 4.7 17.7 38.9 66.5 99.2 100 100 9.9 8.4 6.8 6.1 5.3 16.0 34.1
68.1 98.5 100 2 10 9.0 6.2 4.7 4.2 5.4 19.7 36.0 69.3 99.2 100 20 13.3 8.1 5.6
5.5 6.1 17.3 36.9 69.4 98.5 100 50 8.6 7.3 6.9 5.9 4.7 18.2 34.4 63.9 98.8 100
100 8.6 8.2 7.5 6.4 5.4 13.9 31.4 62.4 98.0 100 1 .7 10 13.6 9.4 6.2 5.0 5.5
20.9 20.9 44.1 87.5 99.3 20 12.8 11.4 9.1 6.2 5.8 17.8 23.7 41.2 84.9 99.4 50
11.7 9.6 10.3 7.2 5.2 15.5 22.6 39.8 80.5 99.2 100 12.6 12.4 7.2 10.1 6.3 13.3
22.2 37.6 72.7 98.7 2 10 12.1 9.0 6.0 4.8 6.1 17.3 19.7 40.1 84.1 98.3 20 12.4
8.8 8.0 6.4 5.9 15.9 19.9 39.0 81.6 98.8 50 11.7 9.6 9.3 7.4 5.6 14.6 19.4
34.0 74.9 98.2 100 11.1 10.1 9.1 9.5 5.7 12.8 21.0 33.3 66.0 97.9

  * 1.

Notes: Size and Power for the different DGPs are reported for 1000
replications. $T=(50,100,200,500)$ is the time series length,
$K=(10,20,50,100)$ the number of variables in the system, the lag length is
fixed to $p=2$ and BIC is used to select the tuninig parameter for the lasso.
$\rho$ indicates the correlation employed to simulate the time series with the
Toeplitz covariance matrix.

##  5 Lag Length Selection

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$
(Zellner and Palm, 1974, Cubadda et al., 2009). On the other hand, Cubadda et
al. (2009) 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

$z_{i,t}=\sum_{j=1}^{p}\gamma_{j,p}z_{i,t-j}+\varepsilon_{i,p,t},\qquad
i=1,\ldots,K,$

by OLS. Then, letting
$\hat{\omega}_{i}=\frac{1}{T}\sum_{t=1}^{T}\hat{\varepsilon}_{i,p,t}^{2}$, we
set
$\hat{\bm{\varOmega}}_{p}=\text{diag}(\hat{\omega}_{1},\ldots,\hat{\omega}_{K})$
and choose the $p$ that minimizes

$\text{IC}^{*}(p)=\ln\left(\det\hat{\bm{\varOmega}}_{p}\right)+C_{T}\frac{pK}{T}=\sum_{i=1}^{N}\log\hat{\omega}_{i}+C_{T}\frac{pK}{T},$

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{\bm{\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 4 for the
combinations of $K=(10,20,50,100)$ and $T=(50,100,200,500,1000)$.666We
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$.

Figure 1: Lag length selection frequencies for AIC∗ and BIC∗

We show the results in Figure 1 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.

##  6 Empirical Application: Causes and Effects of Economic Uncertainty

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
(Ludvigson et al., 2021). 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 Bloom (2009) economic uncertainty has been at the
forefront of macroeconomic debate. Bloom (2009) 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 McCracken and Ng (2016) 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.777In 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 C We extend our analysis to a recent post-
Covid vintage with VIXCLSx included.

An alternative, very popular measure of uncertainty was proposed by Baker et
al. (2016), 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 (McCracken and Ng, 2016) 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 Smeekes and Wijler (2020) and
Smeekes and Wilms (2023), 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.888To run these analyses we used the authors R
package HDGCvar available at https://github.com/Marga8/HDGCvar.

For the Economic Policy Uncertainty index we gather the data directly from
Baker et al. (2016).999The data is freely available on the website of Scott R.
Baker, Nick Bloom and Steven J. Davis at
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.101010We also repeated the analysis using the News-based Policy
Uncertainty index as described in Baker et al. (2016) 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 C 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 3 and 3
we instead just investigate the predictive relations to and from VXOCLSx with
all the macroeconomic series of FRED-MD using PDS-LA-LM.

Figure 2: PDS-LA-LM, Causes of VXOCLSx, $\alpha=0.01$

Figure 3: PDS-LA-LM, Caused by VXOCLSx, $\alpha=0.01$

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.

Figure 4: PDS-LM, Causes of VXOCLSx, $\alpha=0.01$

Figure 5: PDS-LM, Caused by VXOCLSx, $\alpha=0.01$

The difference between Figure 3, 3 and Figure 5, 5 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 6 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%.

Figure 6:  $p$-values for lag-augmented (LA) and stationary transformed (ST)
Granger causality tests between FRED-MD and VXOCLSx. Full bars equal a
$p$-value of 0, empty bars 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.

Figure 7: PDS-LA-LM, Causes of US-EPU, $\alpha=0.01$

Figure 8: PDS-LA-LM, Caused by US-EPU, $\alpha=0.01$

In Figure 8 we only find five connections when it come to Granger causality
relations from FRED-MD series to US-EPU. Vice-versa, in Figure 8 there are 13
causal paths from US-EPU to the other macroeconomic series. The results are
quantitatively similar to FRED$\leftrightarrow$VXOCLSx in Figure 3, 3, 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 8 that are not occurring in Figure 3 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.

Figure 9: PDS-LM, Causes of US-EPU, $\alpha=0.01$

Figure 10: PDS-LM, Caused by US-EPU, $\alpha=0.01$

As before, in Figure 10, 10 we repeat the analysis, this time around however
we apply the recommended stationary transformations (ST) from FRED-MD and we
use Hecq et al. (2021) 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 11 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.

Figure 11:  $p$-values for lag-augmented (LA) and stationary transformed (ST)
Granger causality tests between FRED-MD and US-EPU. Full bars equal a
$p$-value of 0, empty bars a $p$-value above 10%

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 Ludvigson et al. (2021) 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.

##  7 Conclusion

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 Toda and
Yamamoto (1995) 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.

## References

  * Adamek et al. (2023) Adamek, R., Smeekes, S., and Wilms, I. (2023).  Lasso inference for high-dimensional time series.  Journal of Econometrics, 235(2):1114–1143.
  * Baker et al. (2016) Baker, S. R., Bloom, N., and Davis, S. J. (2016).  Measuring economic policy uncertainty.  Quarterly Journal of Economics, 131(4):1593–1636.
  * Belloni et al. (2014) Belloni, A., Chernozhukov, V., and Hansen, C. (2014).  Inference on treatment effects after selection among high-dimensional controls.  Review of Economic Studies, 81(2):608–650.
  * Billio et al. (2012) Billio, M., Getmansky, M., Lo, A. W., and Pelizzon, L. (2012).  Econometric measures of connectedness and systemic risk in the finance and insurance sectors.  Journal of Financial Economics, 104(3):535–559.
  * Bloom (2009) Bloom, N. (2009).  The impact of uncertainty shocks.  Econometrica, 77(3):623–685.
  * Chang and Dickey (1994) Chang, M. C. and Dickey, D. A. (1994).  Recognizing overdifferenced time series.  Journal of Time Series Analysis, 15(1):1–18.
  * Cochrane (1991) Cochrane, J. H. (1991).  A critique of the application of unit root tests.  Journal of Economic Dynamics and Control, 15(2):275–284.
  * Cubadda et al. (2009) Cubadda, G., Hecq, A., and Palm, F. C. (2009).  Studying co-movements in large multivariate data prior to multivariate modelling.  Journal of Econometrics, 148(1):25–35.
  * Dickey and Fuller (1979) Dickey, D. A. and Fuller, W. A. (1979).  Distribution of the estimators for autoregressive time series with a unit root.  Journal of the American Statistical Association, 74(366a):427–431.
  * Dolado and Lütkepohl (1996) Dolado, J. J. and Lütkepohl, H. (1996).  Making wald tests work for cointegrated var systems.  Econometric Reviews, 15(4):369–386.
  * Friston et al. (2013) Friston, K., Moran, R., and Seth, A. K. (2013).  Analysing connectivity with granger causality and dynamic causal modelling.  Current opinion in neurobiology, 23(2):172–178.
  * Fuller (2009) Fuller, W. A. (2009).  Introduction to statistical time series, volume 428.  John Wiley & Sons.
  * Granger (1969) Granger, C. W. (1969).  Investigating causal relations by econometric models and cross-spectral methods.  Econometrica, pages 424–438.
  * Granger (1980) Granger, C. W. (1980).  Testing for causality: a personal viewpoint.  Journal of Economic Dynamics and control, 2:329–352.
  * Hamilton (1994) Hamilton, J. D. (1994).  Time Series Analysis, volume 2.  Princeton University Press.
  * Hecq et al. (2021) Hecq, A., Margaritella, L., and Smeekes, S. (2021).  Granger Causality Testing in High-Dimensional VARs: A Post-Double-Selection Procedure.  Journal of Financial Econometrics, 21(3):915–958.
  * Inoue and Kilian (2020) Inoue, A. and Kilian, L. (2020).  The uniform validity of impulse response inference in autoregressions.  Journal of Econometrics, 215(2):450–472.
  * Johansen (1992) Johansen, S. (1992).  A representation of vector autoregressive processes integrated of order 2.  Econometric Theory, 8(2):188–202.
  * Johansen (1995) Johansen, S. (1995).  Likelihood-based inference in cointegrated vector autoregressive models.  OUP Oxford.
  * Kiviet (1986) Kiviet, J. F. (1986).  On the rigour of some misspecification tests for modelling dynamic relationships.  The Review of Economic Studies, 53(2):241–261.
  * Kock and Callot (2015) Kock, A. B. and Callot, L. (2015).  Oracle inequalities for high dimensional vector autoregressions.  Journal of Econometrics, 186(2):325–344.
  * Leeb and Pötscher (2005) Leeb, H. and Pötscher, B. M. (2005).  Model selection and inference: Facts and fiction.  Econometric Theory, 21:21–59.
  * Ludvigson et al. (2021) Ludvigson, S. C., Ma, S., and Ng, S. (2021).  Uncertainty and business cycles: exogenous impulse or endogenous response?  American Economic Journal: Macroeconomics, 13(4):369–410.
  * Lütkepohl (2005) Lütkepohl, H. (2005).  New Introduction to Multiple Time Series Analysis.  Springer Science & Business Media.
  * Masini et al. (2022) Masini, R. P., Medeiros, M. C., and Mendes, E. F. (2022).  Regularized estimation of high-dimensional vector autoregressions with weakly dependent innovations.  Journal of Time Series Analysis, 43(4):532–557.
  * McCracken and Ng (2016) McCracken, M. W. and Ng, S. (2016).  FRED-MD: A monthly database for macroeconomic research.  Journal of Business & Economic Statistics, 34(4):574–589.
  * Mei and Shi (2022) Mei, Z. and Shi, Z. (2022).  On LASSO for high dimensional predictive regression.  arXiv e-print 2212.07052.
  * Miao et al. (2023) Miao, K., Phillips, P. C., and Su, L. (2023).  High-dimensional vars with common factors.  Journal of Econometrics, 233(1):155–183.
  * Montiel Olea and Plagborg-Møller (2021) Montiel Olea, J. L. and Plagborg-Møller, M. (2021).  Local projection inference is simpler and more robust than you think.  Econometrica, 89(4):1789–1823.
  * Park and Phillips (1988) Park, J. Y. and Phillips, P. C. (1988).  Statistical inference in regressions with integrated processes: Part 1\.  Econometric Theory, 4(3):468–497.
  * Phillips and Durlauf (1986) Phillips, P. C. and Durlauf, S. N. (1986).  Multiple time series regression with integrated processes.  The Review of Economic Studies, 53(4):473–495.
  * Phillips and Hansen (1990) Phillips, P. C. and Hansen, B. E. (1990).  Statistical inference in instrumental variables regression with i (1) processes.  Review of Economic Studies, 57(1):99–125.
  * Phillips and Solo (1992) Phillips, P. C. and Solo, V. (1992).  Asymptotics for linear processes.  The Annals of Statistics, pages 971–1001.
  * Seth et al. (2015) Seth, A. K., Barrett, A. B., and Barnett, L. (2015).  Granger causality analysis in neuroscience and neuroimaging.  Journal of Neuroscience, 35(8):3293–3297.
  * Sims et al. (1990) Sims, C. A., Stock, J. H., Watson, M. W., et al. (1990).  Inference in linear time series models with some unit roots.  Econometrica, 58(1):113–144.
  * Smeekes and Wijler (2020) Smeekes, S. and Wijler, E. (2020).  Unit roots and cointegration.  In Fuleky, P., editor, Macroeconomic Forecasting in the Era of Big Data, volume 52 of Advanced Studies in Theoretical and Applied Econometrics, chapter 17, pages –541–584. Springer.
  * Smeekes and Wijler (2021) Smeekes, S. and Wijler, E. (2021).  An automated approach towards sparse single-equation cointegration modelling.  Journal of Econometrics, 221(1):247–276.
  * Smeekes and Wilms (2023) Smeekes, S. and Wilms, I. (2023).  bootur: An r package for bootstrap unit root tests.  Journal of Statistical Software, 106(12):1–39.
  * Stern and Kaufmann (2014) Stern, D. I. and Kaufmann, R. K. (2014).  Anthropogenic and natural causes of climate change.  Climatic Change, 122(1-2):257–269.
  * Stock and Watson (1989) Stock, J. H. and Watson, M. W. (1989).  Interpreting the evidence on money-income causality.  Journal of Econometrics, 40(1):161–181.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society - Series B, pages 267–288.
  * Toda and Yamamoto (1995) Toda, H. Y. and Yamamoto, T. (1995).  Statistical inference in vector autoregressions with possibly integrated processes.  Journal of Econometrics, 66(1-2):225–250.
  * Wijler (2022) Wijler, E. (2022).  A restricted eigenvalue condition for unit-root non-stationary data.  (2208.12990).
  * Wooldridge (1987) Wooldridge, J. M. (1987).  A regression based lagrange multiplier statistic that is robust in the presence of heteroskedasticity.
  * Zellner and Palm (1974) Zellner, A. and Palm, F. C. (1974).  Time series analysis and simultaneous equation econometric models.  Journal of Econometrics, 2(1):17–54.
  * Zhang et al. (2019) Zhang, R., Robinson, P., and Yao, Q. (2019).  Identifying cointegration by eigenanalysis.  Journal of the American Statistical Association, 114(526):916–927.

##  Appendix A Main Theoretical Results

###  A.1 Justification of Lag-Augmentation

By adding a “free” lag of $x_{t}$, we in essence allow the system to
“difference itself” into the correct order to remove the unit roots. We will
now make this argument more precise using some algebra manipulations similar
to Toda and Yamamoto (1995). We will focus on the bivariate case, but give
some pointers where needed to extend the approach to multiple variables of
interest. In the following, we will generally focus on the case where $d\leq
2$. While one can allow for higher orders $d$ in the estimation procedure,
this would typically not be required for most economic applications of
interest, where it is commonly accepted that the maximum order of integration
of economic data hardly ever exceeds 2. Hence, we confine our presentation up
to the case of $I(2)$ series.

Let us introduce the following $(p+d)\times(p+d)$ upper-triangular
transformation matrices:

$\bm{R}=\begin{bmatrix}1&1&1&\cdots&1\\\ 0&1&1&\cdots&1\\\ 0&0&1&\cdots&1\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\cdots&1\end{bmatrix},\quad\bm{P}_{d}=\bm{R}^{d}=\begin{bmatrix}1&2&3&\cdots&p+d\\\
0&1&2&\cdots&p+d-1\\\ 0&0&1&\cdots&p+d-2\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ 0&0&0&0&\cdots&1\end{bmatrix}.$

$\bm{P}_{d}$ can be considered as an “integration matrix” which integrates the
vector it is multiplied with to order $d$. Similarly, we can consider the
inverse matrix $\bm{P}_{d}^{-1}$ as a difference matrix, which takes
differences $d$ times. Now define the transformed data and parameters

$\begin{split}&\begin{bmatrix}\bm{\beta}_{+}^{*}\\\
\bm{\delta}\end{bmatrix}=\begin{bmatrix}\bm{P}_{d}&\bm{0}\\\
\bm{0}&\bm{I}_{p(K-1)}\end{bmatrix}\begin{bmatrix}\bm{\beta}_{+}\\\
\bm{\delta}_{2}\end{bmatrix},\qquad\begin{bmatrix}\bm{X}_{-(p+d)}^{*}\\\
\bm{V}\end{bmatrix}=\begin{bmatrix}\bm{X}_{-(p+d)}\\\
\bm{V}\end{bmatrix}\begin{bmatrix}\bm{P}_{d}&\bm{0}\\\
\bm{0}&\bm{I}_{p(K-1)}\end{bmatrix}^{-1}.\end{split}$ (A.1)

Then we can rewrite (4) as

$\displaystyle\bm{y}$
$\displaystyle=\bm{X}_{-(p+d)}\bm{P}_{d}^{-1}\bm{P}_{d}\bm{\beta}_{+}+\bm{V}\bm{P}_{d}^{-1}\bm{P}_{d}\bm{\delta}+\bm{u},=$
(A.2)
$\displaystyle=\bm{X}_{-(p+d)}^{*}\bm{\beta}_{+}^{*}+\bm{V}\bm{\delta}+\bm{u}.$

$\bm{X}_{-(p+d)}^{*}=\bm{P}_{d}^{-1}\bm{X}_{-(p+d)}$ now contains $d$-th order
differences for the first $p$ entries, which are $I(0)$ if the original data
are $I(d)$. Furthermore, performing the test of whether the first $p$ elements
of $\bm{\beta}_{+}$ are zero has a one-to-one correspondence to a test on the
first $p$ elements of $\bm{\beta}_{+}^{*}$. This then allows us to transform
the test of Granger non-causality on $I(d)$ into an equivalent test on $I(0)$
data, thus leading to a test with a standard limit distribution regardless of
the order of integration of the data.

We will now show that the null hypothesis of Granger non-causality in (A.2)
has a one-to-one correspondence to a test on the first $p$ elements of
$\bm{\beta}_{+}^{*}$, denoted as $\bm{\beta}_{p}^{*}$.

######  Lemma A.1.

Consider testing (3) in (2). This is equivalent to testing

$H_{0}:\bm{\beta}_{p}^{*}=\bm{0}\quad\text{against}\quad
H_{1}:\bm{\beta}_{p}^{*}\neq\bm{0}$ (A.3)

in (A.2)

###### Proof of Lemma A.1.

Define the selection matrix
$\bm{M}=\left(\begin{matrix}\bm{I}_{p}&\bm{0}_{p\times d}\\\ \bm{0}_{d\times
p}&\bm{0}_{d\times d}\end{matrix}\right)$. Let $\bm{R}_{p\times p}$ denote the
upper left $p\times p$ block of $\bm{R}$, and note that

$\bm{M}\bm{P}_{d}=\bm{M}\bm{R}^{d}=\left(\begin{matrix}\bm{R}_{p\times
p}&\bm{0}_{{p\times d}}\\\ \bm{0}_{d\times p}&\bm{0}_{d\times
d}\end{matrix}\right)^{d}=(\bm{M}\bm{R})^{d}.$

Therefore, we can write

$\left(\begin{matrix}\bm{\beta}_{p}^{*}\\\ \bm{0}_{d\times
1}\end{matrix}\right)=\bm{M}\bm{\beta}_{+}^{*}=\bm{M}\bm{P}_{d}\bm{\beta}_{+}=(\bm{M}\bm{R})^{d}\bm{\beta}_{+}=\left(\begin{matrix}\bm{R}_{p\times
p}^{d}\bm{\beta}\\\ \bm{0}_{d\times 1}\end{matrix}\right).$ (A.4)

As $\bm{R}_{p\times p}$ is invertible, it follows that the test of
$\bm{\beta}_{p}=\bm{0}\Longleftrightarrow\bm{\beta}_{p}^{*}=\bm{R}_{p\times
p}\bm{0}=\bm{0}$. ∎

The following examples clarify how the lag augmentation allows the interest
variables to be re-expressed as difference-stationary. Let $p=d=1$ and let
$\mathbb{x}_{t}^{*}$ be a row of $\bm{X}_{-(p+d)}^{*}$. Then

$\bm{R}=\bm{P}_{1}=\begin{bmatrix}1&1\\\
0&1\end{bmatrix};\;\bm{P}_{1}^{-1}=\begin{bmatrix}[r]1&-1\\\
0&1\end{bmatrix};\;\bm{P}_{1}^{-1}\mathbb{x}_{t}^{*^{\prime}}=\begin{bmatrix}[r]1&-1\\\
0&1\end{bmatrix}\begin{bmatrix}x_{t-1}\\\
x_{t-2}\end{bmatrix}=\begin{bmatrix}\Delta x_{t-1}\\\ x_{t-2}\end{bmatrix}.$

For $p=2,d=1$ we get

$\bm{R}=\bm{P}_{1}=\begin{bmatrix}1&1&1\\\ 0&1&1\\\
0&0&1\end{bmatrix};\;\bm{P}_{1}^{-1}=\begin{bmatrix}[r]1&-1&0\\\ 0&1&-1\\\
0&0&1\end{bmatrix};\;\bm{P}_{1}^{-1}\mathbb{x}_{t}^{*^{\prime}}=\begin{bmatrix}[r]1&-1&0\\\
0&1&-1\\\ 0&0&1\end{bmatrix}\begin{bmatrix}x_{t-1}\\\ x_{t-2}\\\
x_{t-3}\end{bmatrix}=\begin{bmatrix}\Delta x_{t-1}\\\ \Delta x_{t-2}\\\
x_{t-3}\end{bmatrix}.$

As a final example, let $p=d=2$:

$\displaystyle\bm{R}=\begin{bmatrix}1&1&1&1\\\ 0&1&1&1\\\ 0&0&1&1\\\
0&0&0&1\end{bmatrix};\;\bm{P}_{2}=\bm{R}\bm{R}=\begin{bmatrix}1&2&3&4\\\
0&1&2&3\\\ 0&0&1&2\\\
0&0&0&1\end{bmatrix};\;\bm{P}_{2}^{-1}=\begin{bmatrix}[r]1&-2&1&0\\\
0&1&-2&1\\\ 0&0&1&-2\\\ 0&0&0&1\end{bmatrix};\;$
$\displaystyle\bm{P}_{2}^{-1}\mathbb{x}_{t}^{*^{\prime}}=\begin{bmatrix}[r]1&-2&1&0\\\
0&1&-2&1\\\ 0&0&1&-2\\\ 0&0&0&1\end{bmatrix}\begin{bmatrix}x_{t-1}\\\
x_{t-2}\\\ x_{t-3}\\\ x_{t-4}\end{bmatrix}=\begin{bmatrix}\Delta^{2}x_{t-1}\\\
\Delta^{2}x_{t-2}\\\ x_{t-3}-2x_{t-4}\\\ x_{t-4}\end{bmatrix}.$

We now show the numerical equivalence of the Wald and LM test statistic for
testing (3) in (4) and (A.3) in (A.2). Let $\bm{W}_{S}$ denote a subset $S$ of
the columns of $\bm{W}_{-p}$. In the low-dimensional case $S$ may simply
contain all columns of $\bm{W}_{-p}$, whereas in the high-dimensional case $S$
is determined through some sort of (double) selection procedure. Let
$\bm{V}_{S}=(\bm{Y}_{-p},\bm{W}_{S})$ and
$\bm{\delta}_{S}=(\bm{\delta}_{1},\bm{\delta}_{2,S})$. Let the coefficients of
the regression

$\bm{Y}=\bm{X}_{-(p+d)}\bm{\beta}_{+}+\bm{V}_{S}\bm{\delta}_{S}+\bm{u},$ (A.5)

be denoted as
$(\hat{\bm{\beta}}_{+}^{\prime},\hat{\bm{\delta}}_{S}^{\prime})^{\prime}$ and
the coefficients of the regression

$\bm{Y}=\bm{X}_{-(p+d)}^{*}\bm{\beta}_{+}^{*}+\bm{V}_{S}\bm{\delta}_{S}+\bm{u},$
(A.6)

be denoted as
$(\hat{\bm{\beta}}_{+}^{*\prime},\hat{\bm{\delta}}_{S}^{*\prime})^{\prime}$.

######  Lemma A.2.

Let $\hat{\bm{\beta}}_{p}$ denote the first $p$ elements of
$\hat{\bm{\beta}}_{+}$, the OLS estimator of $\bm{\beta}_{+}$ in (A.5), and
$\hat{\bm{\beta}}_{p}^{*}$ denote the first $p$ elements of
$\hat{\bm{\beta}}_{+}^{*}$, the OLS estimator of $\bm{\beta}_{+}^{*}$ in
(A.6). Also define

$\hat{\bm{u}}_{+}=\bm{y}-\bm{X}_{-(p+d)}\hat{\bm{\beta}}_{+}-\bm{V}_{S}\hat{\bm{\delta}}_{S},\qquad\hat{\bm{u}}_{+}^{*}=\bm{y}-\bm{X}_{-(p+d)}^{*}\hat{\bm{\beta}}_{+}^{*}-\bm{V}_{S}\hat{\bm{\delta}}_{S}^{*},$

as the corresponding residuals. Then consider the Wald and LM statistics as
defined in Algorithm 1 (with the generic set $S$ replacing
$\hat{\mathcal{S}}$). Then $\hat{\bm{\beta}}_{p}^{*}=\bm{R}_{p\times
p}^{d}\hat{\bm{\beta}}_{p}$, $\text{Wald}=\text{Wald}^{*}$ and
$\text{LM}=\text{LM}^{*}$.

###### Proof of Lemma A.2.

Define the residual-maker
$\mathcal{M}(\bm{X}):=\bm{I}-\bm{X}(\bm{X}^{\prime}\bm{X})^{-1}\bm{X}^{\prime}$.
Then,

$\displaystyle\left(\begin{matrix}\hat{\bm{\beta}}_{p}^{*}\\\ \bm{0}_{d\times
1}\end{matrix}\right)$
$\displaystyle=\bm{M}\hat{\bm{\beta}}_{+}^{*}=\bm{M}\left(\bm{X}_{-(p+d)}^{*\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}^{*}\right)^{-1}\bm{X}_{-(p+d)}^{*\prime}\mathcal{M}(\bm{V}_{S})\bm{y}$
$\displaystyle=\bm{M}\bm{P}_{d}\left(\bm{X}_{-(p+d)}^{\prime}M(\bm{V}_{S})\bm{X}_{-(p+d)}\right)^{-1}\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{y}=\bm{M}\bm{P}_{d}\hat{\bm{\beta}}_{+}=\left(\begin{matrix}\bm{R}_{p\times
p}^{d}\hat{\bm{\beta}}_{p}\\\ \bm{0}_{d\times 1}\end{matrix}\right).$

It therefore follows that $\hat{\bm{\beta}}_{p}^{*}=\bm{R}_{p\times
p}^{d}\hat{\bm{\beta}}_{p}$.

Next, let $\bm{Z}_{+}:=(\bm{X}_{-(p+d)},\bm{V}_{S})$ and
$\bm{Z}_{+}^{*}:=(\bm{X}_{-(p+d)}^{*},\bm{V}_{S})$, and define

$\bm{P}_{d}^{*}=\left(\begin{matrix}\bm{P}_{d}&\bm{0}_{(p+d)\times(K-1)p}\\\
\bm{0}_{(K-1)p\times(p+d)}&\bm{I}_{(K-1)p}\end{matrix}\right).$

It then follows directly that $\bm{Z}_{+}^{*}=\bm{Z}_{+}\bm{P}_{d}^{*}$ and
consequently

$\begin{split}\hat{\bm{u}}_{+}^{*}=M(\bm{Z}_{+}^{*})\bm{y}=\left(\bm{I}-\bm{Z}_{+}\bm{P}_{d}^{*}(\bm{P}_{d}^{*\prime}\bm{Z}_{+}^{\prime}\bm{Z}_{+}\bm{P}_{d}^{*})^{-1}\bm{P}_{d}^{*\prime}\bm{Z}_{+}\right)\bm{y}=M(\bm{Z}_{+})\bm{y}=\hat{\bm{u}}_{+}.\end{split}$

Letting $\bm{M}_{-}=(\bm{I}_{p},\bm{0}_{p\times d})$, then it follows directly
from (A.4) that $\bm{M}_{-}\bm{P}_{d}=\bm{R}_{p\times p}^{d}\bm{M}_{-}$. The
test statistic can be written as

$\begin{split}\text{Wald}^{*}&=\bm{y}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}^{*}\left(\bm{X}_{-(p+d)}^{*\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}^{*}\right)^{-1}\bm{M}_{-}^{\prime}\left[\hat{\sigma}_{u}^{*2}\left(\bm{M}_{-}\bm{X}_{-(p+d)}^{*\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}^{*}\bm{M}_{-}^{\prime}\right)^{-1}\right]^{-1}\\\
&\quad\times\bm{M}_{-}\left(\bm{X}_{-(p+d)}^{*\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}^{*}\right)^{-1}\bm{X}_{-(p+d)}^{*\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}^{*}\bm{y}\\\
&=\hat{\sigma}_{u}^{-2}\bm{y}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\left(\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\right)^{-1}\bm{P}_{d}^{\prime}\bm{M}_{-}^{\prime}\left(\bm{M}_{-}\bm{P}_{d}^{-1\prime}\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\bm{P}_{d}^{-1}\bm{M}_{-}^{\prime}\right)\\\
&\quad\times\bm{M}_{-}\bm{P}_{d}\left(\bm{X}_{-(p+d)}^{\prime}M(\bm{V}_{S})\bm{X}_{-(p+d)}\right)^{-1}\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\bm{y}\\\
&=\hat{\sigma}_{u}^{-2}\bm{y}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\left(\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\right)^{-1}\bm{M}_{-}^{\prime}\bm{R}_{p\times
p}^{d\prime}\\\ &\quad\times\left(\bm{R}_{p\times
p}^{d\prime}\right)^{-1}\bm{M}_{-}\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\bm{M}_{-}^{\prime}\left(\bm{R}_{p\times
p}^{d}\right)^{-1}\\\ &\quad\times\bm{R}_{p\times
p}^{d}\bm{M}_{-}\left(\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\right)^{-1}\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\bm{y}\\\
&=\hat{\sigma}_{u}^{-2}\bm{y}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\left(\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\right)^{-1}\bm{M}_{-}^{\prime}\bm{M}_{-}\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\bm{M}_{-}^{\prime}\\\
&\quad\times\bm{M}_{-}\left(\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\right)^{-1}\bm{X}_{-(p+d)}^{\prime}\mathcal{M}(\bm{V}_{S})\bm{X}_{-(p+d)}\bm{y}=\hat{\bm{\beta}}_{p}^{\prime}\left(\widehat{\mathbb{V}ar}\hat{\bm{\beta}}_{p}\right)^{-1}\hat{\bm{\beta}}_{p}=\text{Wald}.\end{split}$

Let $\bm{V}_{+}=(\bm{x}_{-(p+1)},\ldots,\bm{x}_{-(p+d)},\bm{V})$ and
$\bm{V}_{+}^{*}=(\bm{x}_{-(p+1)}^{*},\ldots,\bm{x}_{-(p+d)}^{*},\bm{V})$.
Defining $\tilde{\sigma}_{u}^{2}=\hat{\bm{\xi}}^{\prime}\hat{\bm{\xi}}$, we
can write the LM test as

$\begin{split}\text{LM}&=\tilde{\sigma}_{u}^{-2}\bm{y}^{\prime}\mathcal{M}(\bm{V}_{+,S})\bm{X}_{-p}\left[\bm{X}_{-p}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{X}_{-p}\right]^{-1}\bm{X}_{-p}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{y}\\\
&=\text{Wald}\frac{\tilde{\sigma}_{u}^{2}}{\hat{\sigma}_{u}^{2}}.\end{split}$

As we have already shown that $\text{Wald}=\text{Wald}^{*}$ and
$\hat{\sigma}_{u}^{2}=\hat{\sigma}_{u}^{*2}$, it remains to show that
$\tilde{\sigma}_{u}^{2}=\tilde{\sigma}_{u}^{*2}$. From the definition of
$\bm{P}_{d}$ we have that $\bm{x}_{-(p+1)}^{*},\ldots,\bm{x}_{-(p+d)}^{*}$ are
a linear combination of only $\bm{x}_{-(p+1)},\ldots,\bm{x}_{-(p+d)}$.
Therefore,
$\mathcal{M}(\bm{V}_{+,S})\bm{y}=\mathcal{M}(\bm{V}_{+,S}^{*})\bm{y}$ from
which the result follows. This completes the proof of the lemma. ∎

###  A.2 Proofs of Main Results

###### Proof of Theorem 1.

Let $\mathcal{P}(\bm{A})=\bm{A}(\bm{A}^{\prime}\bm{A})^{-1}\bm{A}^{\prime}$
denote the projection on the space spanned by $\bm{A}$ and let
$\mathcal{M}(\bm{A})=I-\mathcal{P}(\bm{A})$ denote the corresponding residual-
maker matrix. For any matrix $\bm{A}$, let the norm
${\left\lVert\cdot\right\rVert}_{p}$ represent the induced $l_{p}$-matrix norm
${\left\lVert\bm{A}\right\rVert}_{p}=\sup_{x\neq
0}{\left\lVert\bm{A}x\right\rVert}_{p}/{\left\lVert x\right\rVert}_{p}$.

We are working on the model

$\bm{y}=\bm{X}_{-p}\bm{\beta}+\bm{V}\bm{\delta}+\bm{u}=\bm{X}_{-(p+d)}\bm{\beta}_{+}+\bm{V}\bm{\delta}+\bm{u}=\bm{X}_{-(p+d)}^{*}\bm{\beta}^{*}+\bm{V}\bm{\delta}+\bm{u}.$
(A.7)

Define the set $\hat{\mathcal{S}}_{j}^{+}$ as the union of $\hat{\mathcal{S}}$
and all lags of $\bm{x}$ except the $j$-th, such that the matrix
$\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}}=(\bm{x}_{-1},\ldots,\bm{x}_{-j+1},\bm{x}_{-j-1},\ldots,\bm{x}_{-(p+d)},\bm{V}_{\hat{\mathcal{S}}})$
contains all variables used in in the second-stage regression except
$\bm{x}_{-j}$. Similarly, define the set $\mathcal{S}_{j}^{+}$ as the union of
$\mathcal{S}$ and all lags of $\bm{x}$ except the $j$-th, such that the matrix
$\bm{Z}_{\mathcal{S}_{j}^{+}}=(\bm{x}_{-1},\ldots,\bm{x}_{-j+1},\bm{x}_{-j-1},\ldots,\bm{x}_{-(p+d)},\bm{V}_{\mathcal{S}})$
contains all variables used in an ‘oracle’ version of the second-stage
regression. For the $j$-th coefficient in $\hat{\bm{\beta}}_{+}$
($\bm{\beta}_{+}$), we can then write

$\begin{split}\sqrt{T}(\hat{\beta}_{j,+}-\beta_{j+})&=\underbrace{\left(T^{-2a}(\bm{x}_{-j}^{\prime}\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{x}_{-j})\right)^{-1}}_{A_{T}^{-1}}\underbrace{T^{1/2-2a}\bm{x}_{-j}^{\prime}\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\left[\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}_{0,-j}+\bm{u}\right]}_{B_{T}},\end{split}$

where $\bm{\eta}_{0,-j}$ is equal to $\bm{\eta}_{0}$ but with a zero inserted
for the coefficient of $\bm{x}_{-j}$, and $a=1/2$ if $\bm{x}$ is I(0), $a=1$
if $\bm{x}$ is I(1) and $a=2$ if $\bm{x}$ is I(2).

We are now going to show in turn that

$\displaystyle
A_{T,j}=T^{-2a}\bm{x}_{-j}^{\prime}\mathcal{M}(\bm{Z}_{\mathcal{S}_{j}^{+}})\bm{x}_{-j}+o_{p}(1),\qquad
j=1,\ldots,p+d,$ (A.8)

and

$\displaystyle
B_{T,j}=T^{1/2-2a}\bm{x}_{-j}^{\prime}\mathcal{M}(\bm{Z}_{\mathcal{S}_{j}^{+}})\left[\bm{V}\bm{\delta}+\bm{u}\right]+o_{p}(1),\qquad
j=1,\ldots,p+d.$ (A.9)

We write the first-stage regressions as

$\begin{split}\bm{y}&=\bm{Z}_{-p}\bm{\eta}_{0}+\bm{e}_{0},\qquad\bm{x}_{-j}=\bm{Z}_{-p}\bm{\eta}_{j}+\bm{e}_{j},\quad
j=1,\ldots,p.\end{split}$

Note that here for notational convenience, in a slight abuse of notation we
insert $\bm{x}_{j}$ on the right-hand side into $\bm{Z}_{-p}$ instead of
considering $\bm{Z}_{-p,\\{-j\\}}$ as we do in the actual regression. Thereby,
we implicitly increase the vector $\bm{\eta}_{j}$ with one entry for
$\bm{x}_{-j}$, for which we insert a zero.

We now prove (A.8). We have that

$\displaystyle A_{T,j}$
$\displaystyle=T^{-2a}\left(\bm{Z}_{-p}\bm{\eta}_{j}+\bm{e}_{j}\right)^{\prime}\left[\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})-\mathcal{M}(\bm{V}_{\mathcal{S}_{j}^{+}})\right]\left(\bm{Z}_{-p}\bm{\eta}_{j}+\bm{e}_{j}\right)$
$\displaystyle=\underbrace{T^{-2a}\bm{\eta}_{j}^{\prime}\bm{Z}_{-p}^{\prime}\left[\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})-\mathcal{M}(\bm{Z}_{\mathcal{S}_{j}^{+}})\right]\bm{Z}_{-p}\bm{\eta}_{j}}_{A_{T,j,1}}$
$\displaystyle\quad+\underbrace{T^{-2a}\bm{e}_{j}^{\prime}\left[\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})-\mathcal{M}(\bm{Z}_{\mathcal{S}_{j}^{+}})\right]\bm{e}_{j}}_{A_{T,j,2}}+\underbrace{2T^{-2a}\bm{\eta}_{j}^{\prime}\bm{Z}_{-p}^{\prime}\left[\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})-\mathcal{M}(\bm{Z}_{\mathcal{S}_{j}^{+}})\right]\bm{e}_{j}}_{A_{T,j,3}}.$

For $A_{T,j,1}$, observe how $\bm{\eta}_{j}$ is zero outside the active set
$\mathcal{S}_{j}$ which is in turn a subset of $\mathcal{S}_{j}^{+}$, such
that we may write
$\bm{Z}_{-p}\bm{\eta}_{j}=\bm{Z}_{\mathcal{S}_{j}^{+}}\bm{\eta}_{j,\mathcal{S}_{j}^{+}}$.
It then directly follows that
$\mathcal{M}(\bm{Z}_{\mathcal{S}_{j}^{+}})\bm{Z}_{-p}\bm{\eta}_{j}=\bm{0}$.
Hence, $A_{T,j,1}$ reduces to

$\displaystyle{\left\lvert A_{T,j,1}\right\rvert}$
$\displaystyle={\left\lvert\bm{\eta}_{j}^{\prime}\bm{Z}_{-p}^{\prime}\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{Z}_{-p}\bm{\eta}_{j}\right\rvert}={\left\lVert\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{Z}_{-p}\bm{\eta}_{j}\right\rVert}_{2}^{2}\leq{\left\lVert\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}})\bm{Z}_{-p}\bm{\eta}_{j}\right\rVert}_{2}^{2}$
$\displaystyle\leq\min_{\underline{\bm{\eta}}:\underline{\bm{\eta}}_{m}=0,m\notin\hat{\mathcal{S}}_{j}}{\left\lVert\bm{Z}_{-p}\bm{\eta}_{j}-\bm{Z}_{\hat{\mathcal{S}}_{j}}\bm{\eta}\right\rVert}_{2}^{2}\leq{\left\lVert\bm{Z}_{-p}(\bm{\eta}_{j}-\hat{\bm{\eta}}_{j})\right\rVert}_{2}^{2}\leq\delta_{T}^{2}T^{1/2},$

which follows from observing that
$\hat{\mathcal{S}}_{j}\subseteq\hat{\mathcal{S}}_{j}^{+}$, and that the
constraint in the minimization is satisfied given
$\hat{\mathcal{S}}_{j}=\\{m:\hat{\eta}_{j,m}=0\\}$ and where the last
inequality follows with probability $1-\Delta_{T}$ from Assumption 2(3).

For $A_{T,j,2}$, note that the difference between the residuals is equal to
the difference between the corresponding projections, such that we can
simplify the expression to

$\displaystyle{\left\lvert A_{T,j,2}\right\rvert}$
$\displaystyle={\left\lvert\bm{e}_{j}^{\prime}\left(\mathcal{P}(\bm{Z}_{\mathcal{S}_{j}^{+}})-\mathcal{P}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\right)\bm{e}_{j}\right\rvert}\leq\underbrace{{\left\lvert\bm{e}_{j}^{\prime}\left(\mathcal{P}(\bm{Z}_{\mathcal{S}_{j}^{+}}))\right)\bm{e}_{j}\right\rvert}}_{A_{T,j,2}^{1}}+\underbrace{{\left\lvert\bm{e}_{j}^{\prime}\left(\mathcal{P}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\right)\bm{e}_{j}\right\rvert}}_{A_{T,j,2}^{2}}.$

For $j=1,\ldots,p$, we then find that

$\displaystyle A_{T,j,2}^{1}$
$\displaystyle={\left\lvert\bm{e}_{j}^{\prime}\bm{Z}_{\mathcal{S}_{j}^{+}}\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\left[\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\bm{Z}_{\mathcal{S}_{j}^{+}}^{\prime}\bm{Z}_{\mathcal{S}_{j}^{+}}\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\right]^{-1}\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\bm{Z}_{\mathcal{S}_{j}^{+}}^{\prime}\bm{e}_{j}\right\rvert}$
$\displaystyle\overset{(I)}{\leq}\sqrt{\bar{s}_{T}}{\left\lVert\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\bm{Z}_{\mathcal{S}_{j}^{+}}^{\prime}\bm{e}_{j}\right\rVert}_{\infty}{\left\lVert\left[\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\bm{Z}_{\mathcal{S}_{j}^{+}}^{\prime}\bm{Z}_{\mathcal{S}_{j}^{+}}\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\right]^{-1}\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\bm{Z}_{\mathcal{S}_{j}^{+}}^{\prime}\bm{e}_{j}\right\rVert}_{2}$
$\displaystyle\overset{(II)}{\leq}\bar{s}_{T}{\left\lVert\bm{D}_{T,\bm{Z}}^{-1}\bm{Z}_{-p}^{\prime}\bm{e}_{j}\right\rVert}_{\infty}{\left\lVert\left[\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\bm{Z}_{\mathcal{S}_{j}^{+}}^{\prime}\bm{Z}_{\mathcal{S}_{j}^{+}}\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\right]^{-1}\right\rVert}_{2}{\left\lVert\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\bm{Z}_{\mathcal{S}_{j}^{+}}^{\prime}\bm{e}_{j}\right\rVert}_{\infty}$
$\displaystyle\overset{(III)}{\leq}\bar{s}_{T}\bm{\kappa}_{T,\min}^{-1}{\left\lVert\bm{D}_{T,\bm{Z}}^{-1}\bm{Z}_{-p}^{\prime}\bm{e}_{j}\right\rVert}_{\infty}^{2}\overset{(IV)}{\leq}\bar{s}_{T}\bar{\gamma}_{T}^{2}\bm{\kappa}_{T,\min}^{-1}\leq\delta_{T}^{2}T^{1/2},$

where (I) follows from the dual norm inequality and bounding the $\ell_{1}$ by
the $\ell_{2}$ norm at the expense of the $\sqrt{\bar{s}_{T}}$ term; (II)
follows from the fact that the columns of $\bm{D}_{T,\bm{Z}}^{-1}\bm{Z}_{-p}$
contain all columns of
$\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\bm{Z}_{\mathcal{S}_{j}^{+}}$, the
definition of the induced norm and bounding the $\ell_{2}$ norm by the
$\ell_{\infty}$ norm. (III) then follows from the eigenvalue condition in
Assumption 2(5) while (IV) invokes the empirical process bound of Assumption
2(1); the last inequality then follows from Assumption 2(6). For
$A_{T,j,2}^{2}$ we can follow the same steps as above to again find that
$A_{T,j,2}^{2}\leq\delta_{T}^{2}T^{1/2}$.

We finally consider $A_{T,j,3}$. By the same arguments as for $A_{T,j,1}$ we
get that $A_{T,j,3}$ reduces to
$A_{T,j,3}={\left\lvert\bm{\eta}_{j}^{\prime}\bm{Z}_{-p}^{\prime}\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{e}_{j}\right\rvert}$.
Let us define the noiseless least squares estimator

$\tilde{\bm{\eta}}_{j,\hat{\mathcal{S}}_{j}^{+}}=\operatorname*{arg\,min}_{\bm{\eta}:\eta_{m}=0,m\notin\hat{\mathcal{S}}_{j}^{+}}{\left\lVert\bm{Z}_{-p}\bm{\eta}_{j}-\bm{Z}_{-p}\bm{\eta}\right\rVert}_{2}^{2},\qquad
j=1,\ldots,p,$

such that
$\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{Z}_{-p}\bm{\eta}_{j}=\bm{Z}_{-p}(\bm{\eta}_{j}-\tilde{\bm{\eta}}_{j,\hat{\mathcal{S}}_{j}^{+}})$.
Then, with probability $1-\Delta_{T}$,

$\displaystyle A_{T,j,3}$
$\displaystyle={\left\lvert(\bm{\eta}_{j}-\tilde{\bm{\eta}}_{j,\hat{\mathcal{S}}_{j}^{+}})^{\prime}\bm{D}_{T,\bm{Z}}\bm{D}_{T,\bm{Z}}^{-1}\bm{Z}_{-p}^{\prime}\bm{e}_{j}\right\rvert}\leq{\left\lVert\bm{D}_{T,\bm{Z}}^{-1}\bm{Z}_{-p}^{\prime}\bm{e}_{j}\right\rVert}_{\infty}{\left\lVert\bm{D}_{T,\bm{Z}}\left(\tilde{\bm{\eta}}_{j,\hat{\mathcal{S}}_{j}^{+}}-\bm{\eta}_{j}\right)\right\rVert}_{1}$
$\displaystyle\leq\bar{\gamma}_{T}\bar{s}_{T}\bm{\kappa}_{T,\min}^{-1}{\left\lVert\bm{Z}_{-p}\left(\tilde{\bm{\eta}}_{j,\hat{\mathcal{S}}_{j}^{+}}-\bm{\eta}_{j}\right)\right\rVert}_{2}\leq\bar{\gamma}_{T}\bar{s}_{T}\bm{\kappa}_{T,\min}^{-1}{\left\lVert\bm{Z}_{-p}\left(\hat{\bm{\eta}}_{j}-\bm{\eta}_{j}\right)\right\rVert}_{2}\leq\delta_{T}^{2}T^{1/2}.$

The first inequality follows from the dual norm inequality. The second
inequality follows from combining the sparsity and the restricted eigenvalue
conditions in Assumption 2(4) and 2(5). In particular, letting
$\tilde{\bm{\eta}}_{j,\hat{\mathcal{S}}_{j}^{+}}^{\dagger}=\bm{D}_{T,\bm{Z}}\tilde{\bm{\eta}}_{\hat{\mathcal{S}}_{j}^{+}}$
and $\bm{\eta}_{j}^{\dagger}=\bm{D}_{T,\bm{Z}}\bm{\eta}_{j}$, we bound
${\left\lVert\tilde{\bm{\eta}}_{j,\hat{\mathcal{S}}_{j}^{+}}^{\dagger}-\bm{\eta}_{j}^{\dagger}\right\rVert}_{2}\leq\bar{s}_{T}{\left\lVert\bm{Z}_{-p}\bm{D}_{T,\bm{Z}}\left(\tilde{\bm{\eta}}_{j,\hat{\mathcal{S}}_{j}^{+}}^{\dagger}-\bm{\eta}_{j}^{\dagger}\right)\right\rVert}_{1}$
as $\bm{D}_{T,\bm{Z}}$ is diagonal and therefore does not affect the sparsity.
We then use that $\tilde{\bm{\eta}}_{j,\hat{\mathcal{S}}_{j}^{+}}$ was
constructed to minimize the $\ell_{2}$ norm of the residuals. Finally, the
last inequality follows the consistency Assumption 2(3) and the rate
Assumption 2(6).

Combining all results, we find that for all $j=1,\ldots,p$, we get that
${\left\lvert
A_{T,j}\right\rvert}\leq\delta_{T}^{2}T^{1/2-2a}\leq\delta_{T}^{2}T^{-1/2}$,
as sufficient for establishing (A.8).

Let us now prove (A.9). By re-arranging terms we have

$\displaystyle B_{T,j}$
$\displaystyle=T^{1/2-2a}\underbrace{\bm{x}_{-j}^{\prime}\left[\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})-\mathcal{M}(\bm{Z}_{\mathcal{S}_{j}^{+}})\right]\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}_{0,-j}}_{B_{T,j,1}}$
$\displaystyle\quad+T^{1/2-2a}\underbrace{\bm{x}_{-j}^{\prime}\left[\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})-\mathcal{M}(\bm{Z}_{\mathcal{S}_{j}^{+}})\right]\bm{u}}_{B_{T,j,2}}.$

We first deal with $B_{T,j,1}$. As $\bm{\eta}_{0}$ is only non-zero on
$\mathcal{S}$, it follows directly that
$\mathcal{M}(\bm{Z}_{\mathcal{S}_{j}^{+}})\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}_{0,-j}=\bm{0}$
and $B_{T,j,1}$ reduces to
$B_{T,j,1}=\bm{x}_{-j}^{\prime}\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{Z}_{-p}\bm{\eta}_{0}$.
As above for (A.8), using the projection
$\bm{x}_{-j}=\bm{Z}_{-p}\bm{\eta}_{j}+\bm{e}_{j}$ we then get the following
expression

$\displaystyle B_{T,j,1}$
$\displaystyle=\underbrace{\bm{\eta}_{j}^{\prime}\bm{Z}_{-p}^{\prime}\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}_{0,-j}}_{B_{T,j,1}^{1}}+\underbrace{\bm{e}_{j}^{\prime}\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}_{0,-j}}_{B_{T,j,1}^{2}}.$
(A.10)

$B_{T,j,1}^{1}$ then follows analogously to $A_{T,j,1}$:

$\displaystyle{\left\lvert B_{T,j,1}^{1}\right\rvert}$
$\displaystyle\leq\sqrt{{\left\lVert\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{Z}_{-p}\bm{\eta}_{j}\right\rVert}_{2}{\left\lVert\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}_{0,-j}\right\rVert}_{2}^{2}}$
$\displaystyle\leq\sqrt{{\left\lVert\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}})\bm{Z}_{-p}\bm{\eta}_{j}\right\rVert}_{2}{\left\lVert\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{0}})\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}_{0,-j}\right\rVert}_{2}^{2}}$
$\displaystyle\leq\sqrt{\min_{\bm{\eta}:\eta_{j,m}=0,m\notin\hat{\mathcal{S}}_{j}}{\left\lVert\bm{Z}_{-p}\bm{\eta}_{j}-\bm{Z}_{\hat{\mathcal{S}}_{j}}\bm{\eta}\right\rVert}_{2}^{2}\min_{\bm{\eta}:\eta_{0,m}=0,m\notin\hat{\mathcal{S}}_{0}}{\left\lVert\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}_{0,-j}-\bm{Z}_{\hat{\mathcal{S}}_{0}}\bm{\eta}\right\rVert}_{2}^{2}}$
$\displaystyle\leq{\left\lVert\bm{Z}_{-p}(\bm{\eta}_{j}-\hat{\bm{\eta}}_{j})\right\rVert}_{2}{\left\lVert\bm{Z}_{-p\backslash\\{-j\\}}(\bm{\eta}_{0,-j}-\hat{\bm{\eta}}_{0})\right\rVert}_{2}\leq\delta_{T}T^{1/4}{\left\lVert\bm{Z}_{-p\backslash\\{-j\\}}(\bm{\eta}_{0,-j}-\hat{\bm{\eta}}_{0})\right\rVert}_{2}.$

Furthermore,

$\displaystyle{\left\lVert\bm{Z}_{-p\backslash\\{-j\\}}(\bm{\eta}_{0,-j}-\hat{\bm{\eta}}_{0})\right\rVert}_{2}\leq{\left\lVert\bm{Z}_{-p}(\bm{\eta}_{0}-\hat{\bm{\eta}}_{0})\right\rVert}_{2}+{\left\lVert\bm{x}_{j}(\eta_{0,j}-\hat{\eta}_{0,j})\right\rVert}_{2}$
$\displaystyle\quad\leq\delta_{T}T^{1/4}+{\left\lvert\hat{\eta}_{0,j}-\eta_{0,j}\right\rvert}\sqrt{\sum_{t=p+d+1-j}^{T-j}x_{t-j}^{2}}\leq\delta_{T}T^{1/4}+\delta_{T}T^{-1/4}O_{p}(T^{a}),$

by the consistency of $\hat{\bm{\eta}}_{0}$ in Assumption 2(3) and since
$\sqrt{\sum_{t=p+d+1-j}^{T-j}x_{t-j}^{2}}=O_{p}(T^{a})$. We may then construct
a sequence, say $\eta_{T}$, such that $\eta_{T}\to\infty$ but
$\eta_{T}\delta_{T}\rightarrow 0$. It then follows that with probability
$1-\Delta_{T}$, we have that $\sqrt{\sum_{t=p+d+1-j}^{T-j}x_{t-j}^{2}}\leq
C\eta_{T}T^{a}$ and it then follows that
$B_{T,j,1}^{1}\leq\delta_{T}^{2}T^{1/2}+\delta_{T}^{2}\eta_{T}T^{a}\leq
2\delta_{T}T^{a}$.

$B_{T,j,1}^{2}$ can be treated analogously to $A_{T,j,3}$. For $j=1,\ldots,p$,
define the estimators

$\tilde{\bm{\eta}}_{0,\hat{\mathcal{S}}_{j}^{+}}=\operatorname*{arg\,min}_{\bm{\eta}:\eta_{m}=0,m\notin\hat{\mathcal{S}}_{j}^{+}}{\left\lVert\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}_{0,-j}-\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}\right\rVert}_{2}^{2},$

such that
$\mathcal{M}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\bm{Z}_{-p\backslash\\{-j\\}}\bm{\eta}_{0,-j}=\bm{Z}_{-p\backslash\\{-j\\}}(\bm{\eta}_{0,-j}-\tilde{\bm{\eta}}_{0,\hat{\mathcal{S}}_{j}^{+}})$.
Then, with probability $1-\Delta_{T}$,

$\displaystyle{\left\lvert B_{T,j,1}^{2}\right\rvert}$
$\displaystyle={\left\lvert(\bm{\eta}_{0,-j}-\tilde{\bm{\eta}}_{0,\hat{\mathcal{S}}_{j}^{+}})^{\prime}\bm{D}_{T,\bm{Z},-j}\bm{D}_{T,\bm{Z},-j}^{-1}\bm{Z}_{-p\backslash\\{-j\\}}^{\prime}\bm{e}_{j}\right\rvert}$
$\displaystyle\leq{\left\lVert\bm{D}_{T,\bm{Z},-j}^{-1}\bm{Z}_{-p\backslash\\{-j\\}}^{\prime}\bm{e}_{j}\right\rVert}_{\infty}{\left\lVert\bm{D}_{T,\bm{Z},-j}\left(\tilde{\bm{\eta}}_{0,\hat{\mathcal{S}}_{j}^{+}}-\bm{\eta}_{0,-j}\right)\right\rVert}_{1}$
$\displaystyle\leq\bar{\gamma}_{T}\bar{s}_{T}\bm{\kappa}_{T,\min}^{-1}{\left\lVert\bm{Z}_{-p}\left(\hat{\bm{\eta}}_{0}-\bm{\eta}_{0}\right)\right\rVert}_{2}\leq\delta_{T}^{2}T^{1/2}.$

For $B_{T,j,2}$, analogously to $A_{T,j,2}$ we have ${\left\lvert
B_{T,j,2}\right\rvert}\leq\underbrace{{\left\lvert\bm{e}_{j}^{\prime}\left(\mathcal{P}(\bm{Z}_{\mathcal{S}_{j}^{+}})\right)\bm{e}_{j}\right\rvert}}_{B_{T,j,2}^{1}}+\underbrace{{\left\lvert\bm{e}_{j}^{\prime}\left(\mathcal{P}(\bm{Z}_{\hat{\mathcal{S}}_{j}^{+}})\right)\bm{e}_{j}\right\rvert}}_{B_{T,j,2}^{2}}$.
Then

$\displaystyle B_{T,j,2}^{1}$
$\displaystyle\leq\bar{s}_{T}{\left\lVert\bm{D}_{T,\bm{Z}}^{-1}\bm{Z}_{-p}^{\prime}\bm{e}_{j}\right\rVert}_{\infty}{\left\lVert\bm{D}_{T,\bm{Z}}^{-1}\bm{Z}_{-p}^{\prime}\bm{u}\right\rVert}_{\infty}{\left\lVert\left[\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\bm{Z}_{\mathcal{S}_{j}^{+}}^{\prime}\bm{Z}_{\mathcal{S}_{j}^{+}}\bm{D}_{T,\mathcal{S}_{j}^{+}}^{-1}\right]^{-1}\right\rVert}_{2}$
$\displaystyle\leq\bar{s}_{T}\bm{\kappa}_{T,\min}^{-1}\bar{s}_{T}\bar{\gamma}_{T}^{2}\bm{\kappa}_{T,\min}^{-1}\leq\delta_{T}^{2}T^{1/2},$

and it directly follows that $B_{T,j,2}^{2}\leq\delta_{T}^{2}T^{1/2}$ as well.

It then follows directly that ${\left\lvert
B_{T,j}\right\rvert}\leq\delta_{T}T^{1/2-a}\leq\delta_{T}^{2}$, thereby
establishing (A.9) for all $j=1,\ldots,p$. ∎

###### Proof of Theorem 2.

We first show that the second stage regression delivers asymptotically normal
estimators due to the lag augmentation. Let
$\bm{V}_{+,\hat{\mathcal{S}}}=(\bm{x}_{-(p+1)},\ldots,\bm{x}_{-(p+d)},\bm{V}_{\hat{\mathcal{S}}})$
and
$\bm{V}_{+,\mathcal{S}}=(\bm{x}_{-(p+1)},\ldots,\bm{x}_{-(p+d)},\bm{V}_{\mathcal{S}})$
and similarly let $\bm{V}_{+,\hat{\mathcal{S}}}^{*}$ and
$\bm{V}_{+,\mathcal{S}}^{*}$ denote their counterparts with
$\bm{x}_{-(p+1)}^{*},\ldots,\bm{x}_{-(p+d)}^{*}$ replacing
$\bm{x}_{-(p+1)},\ldots,\bm{x}_{-(p+d)}$. By virtue of Lemmas A.1 and A.2 as
well as equations (A.8) and (A.9), we have that

$\begin{split}\sqrt{T}(\hat{\bm{\beta}}_{p}^{*}-\bm{\beta}_{p}^{*})&=\sqrt{T}\bm{R}_{p\times
p}(\hat{\bm{\beta}}_{p}-\bm{\beta})=\sqrt{T}\bm{R}_{p\times
p}\left(\bm{X}_{-p}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{X}_{-p}\right)^{-1}\bm{X}_{-p}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{u}\\\
&=\sqrt{T}\bm{R}_{p\times
p}\left(\bm{X}_{-p}^{\prime}\mathcal{M}(\bm{V}_{+,\mathcal{S}})\bm{X}_{-p}\right)^{-1}\bm{X}_{-p}^{\prime}\mathcal{M}(\bm{V}_{+,\mathcal{S}})\bm{u}+o_{p}(1)\\\
&=\big{(}\underbrace{T^{-1}\bm{X}_{-p}^{*\prime}\mathcal{M}(\bm{V}_{+,\mathcal{S}}^{*})\bm{X}_{-p}^{*}\big{)}^{-1}}_{\bm{A}_{T}}\underbrace{T^{-1/2}\bm{X}_{-p}^{*\prime}\mathcal{M}(\bm{V}_{+,\mathcal{S}}^{*})\bm{u}}_{\bm{b}_{T}}+o_{p}(1).\end{split}$
(A.11)

Let $\bm{S}=\bm{V}_{+,\mathcal{S}}\bm{Q}_{\mathcal{S}}$ and
$\bm{G}_{T,\bm{Q}}=\bm{Q}_{\mathcal{S}}\bm{D}_{T,\bm{S}}$, , where
$\bm{D}_{T,\bm{S}}$ is the matrix with $T^{1/2},T$ and $T^{2}$ on the diagonal
for the I(0), I(1) and I(2) elements, respectively. Then we can write

$\begin{split}\bm{A}_{T}&=T^{-1}\bm{X}_{-p}^{*\prime}\bm{X}_{-p}^{*}-T^{-1/2}\bm{X}_{-p}^{*\prime}\bm{V}_{+,\mathcal{S}}^{*}\bm{G}_{T,\bm{Q}}^{-1}\big{(}\bm{G}_{T,\bm{Q}}^{-1}\bm{V}_{+,\mathcal{S}}^{*\prime}\bm{V}_{+,\mathcal{S}}^{*}\bm{G}_{T,\bm{Q}}^{-1}\big{)}^{-1}\bm{G}_{T,\bm{Q}}^{-1}\bm{V}_{+,\mathcal{S}}^{*\prime}\bm{X}_{-p}^{*}T^{-1/2}\\\
&=T^{-1}\bm{X}_{-p}^{*\prime}\bm{X}_{-p}^{*}-T^{-1/2}\bm{X}_{-p}^{*\prime}\bm{S}\bm{D}_{T,\bm{S}}^{-1}\big{(}\bm{D}_{T,\bm{S}}^{-1}\bm{S}^{\prime}\bm{S}\bm{D}_{T,\bm{S}}^{-1}\big{)}^{-1}\bm{D}_{T,\bm{S}}^{-1}\bm{S}^{\prime}\bm{X}_{-p}^{*}T^{-1/2}.\end{split}$

Let $\bm{S}_{-0}=(\bm{S}_{1},\bm{S}_{2})$ and $\bm{D}_{T,\bm{S}_{-0}}$ the
corresponding submatrix of $\bm{D}_{T,\bm{S}}$. It follows from inverting
block matrices that

$\displaystyle\bm{A}_{T}$
$\displaystyle=\underbrace{T^{-1}\bm{X}_{-p}^{*\prime}\mathcal{M}(\bm{S}_{0})\bm{X}_{-p}^{*}}_{\bm{C}_{T,1}}-\underbrace{T^{-1/2}\bm{X}_{-p}^{*\prime}\mathcal{M}(\bm{S}_{0})\bm{S}_{-0}\bm{D}_{T,\bm{S}_{-0}}^{-1}}_{\bm{R}_{T,2}^{\prime}}$
$\displaystyle\quad\times\big{(}\underbrace{\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\mathcal{M}(\bm{S}_{0})\bm{S}_{-0}\bm{D}_{T,\bm{S}_{-0}}^{-1}}_{\bm{R}_{T,1}}\big{)}^{-1}\underbrace{\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\mathcal{M}(\bm{S}_{0})\bm{X}_{-p}^{*}T^{-1/2}}_{\bm{R}_{T,2}}$

Then note that

$\displaystyle\lambda_{\min}(\bm{R}_{T,1})$
$\displaystyle\geq\lambda_{\min}\left(\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\bm{S}_{-0}\bm{D}_{T,\bm{S}_{-0}}^{-1}\right)$
$\displaystyle\quad-\lambda_{\max}\left(\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\bm{S}_{0}\bm{D}_{T,\bm{S}_{0}}^{-1}\big{(}\bm{D}_{T,\bm{S}_{0}}^{-1}\bm{S}_{0}^{\prime}\bm{S}_{0}\bm{D}_{T,\bm{S}_{0}}^{-1}\big{)}^{-1}\bm{D}_{T,\bm{S}_{0}}^{-1}\bm{S}_{0}^{\prime}\bm{S}_{-0}\bm{D}_{T,\bm{S}_{-0}}^{-1}\right)$
$\displaystyle\geq\lambda_{\min}\left(\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\bm{S}_{-0}\bm{D}_{T,\bm{S}_{-0}}^{-1}\right)$
$\displaystyle\quad-{\left\lVert\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\bm{S}_{0}\bm{D}_{T,\bm{S}_{0}}^{-1}\right\rVert}_{2}^{2}/\lambda_{\min}(\bm{D}_{T,\bm{S}_{0}}^{-1}\bm{S}_{0}^{\prime}\bm{S}_{0}\bm{D}_{T,\bm{S}_{0}}^{-1}\big{)},$

where the final inequality uses that, for compatible matrices $\bm{A}$ and
$\bm{B}$, we have
$\lambda_{\max}(\bm{A}^{\prime}\bm{B}^{-1}\bm{A})\leq\lambda_{\max}(\bm{A}^{\prime}\bm{A})\lambda_{\max}(\bm{B}^{-1})={\left\lVert\bm{A}\right\rVert}_{2}^{2}/\lambda_{\min}(\bm{B})$.
It then follows directly from Assumption 3 that

$\displaystyle\lambda_{\min}\left(\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\bm{S}_{-0}\bm{D}_{T,\bm{S}_{-0}}^{-1}\right)\geq\kappa_{T,12},\quad{\left\lVert\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\bm{S}_{0}\bm{D}_{T,\bm{S}_{0}}^{-1}\right\rVert}_{2}^{2}\leq\phi_{T,01}^{2}+\phi_{T,02}^{2},$
$\displaystyle\lambda_{\min}(\bm{D}_{T,\bm{S}_{0}}^{-1}\bm{S}_{0}^{\prime}\bm{S}_{0}\bm{D}_{T,\bm{S}_{0}}^{-1}\big{)}\geq\kappa_{T,0}.$

We may then conclude that

$\lambda_{\min}(\bm{R}_{T,1})\geq\kappa_{T,12}-(\phi_{T,01}^{2}+\phi_{T,02}^{2})/\kappa_{T,0}\geq\kappa_{T,-0}.$

In similar fashion we can bound

$\displaystyle{\left\lVert\bm{R}_{T,2}\right\rVert}_{2}$
$\displaystyle\leq{\left\lVert\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\bm{X}_{-p}^{*}T^{-1/2}\right\rVert}_{2}+{\left\lVert\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\bm{S}_{0}\bm{D}_{T,\bm{S}_{0}}^{-1}\right\rVert}_{2}$
$\displaystyle\quad\times{\left\lVert\bm{D}_{T,\bm{S}_{0}}^{-1}\bm{S}_{0}^{\prime}\bm{X}_{-p}^{*}T^{-1/2}\right\rVert}_{2}/\lambda_{\min}(\bm{D}_{T,\bm{S}_{0}}^{-1}\bm{S}_{0}^{\prime}\bm{S}_{0}\bm{D}_{T,\bm{S}_{0}}^{-1}\big{)}$
$\displaystyle\leq\phi_{T,01}+\phi_{T,02}+2(\phi_{T,01}^{2}+\phi_{T,02}^{2})/\kappa_{T,0}\leq\phi_{T,-0}.$

Letting

$\Sigma_{xx}:=\mathbb{E}\left(T^{-1}\bm{X}_{-p}^{*\prime}\bm{X}_{-p}^{*}\right),\qquad\Sigma_{sx}:=\mathbb{E}\left(T^{-1}\bm{S}_{0}^{\prime}\bm{X}_{-p}^{*}\right),\qquad\Sigma_{ss}:=\mathbb{E}\left(T^{-1}\bm{S}_{0}^{\prime}\bm{S}_{0}\right),$
(A.12)

it follows directly from Assumption 3 that
${\left\lVert\bm{A}_{T,1}-\bm{\varSigma}_{x|w}\right\rVert}_{2}\leq\delta_{T}$
with probability $1-\Delta$, where

$\bm{\varSigma}_{x|s}:=\bm{\varSigma}_{xx}-\bm{\varSigma}_{sx}^{\prime}\bm{\varSigma}_{ss}^{-1}\bm{\varSigma}_{sx}.$
(A.13)

Putting everything together, we then find that
$\bm{A}_{T}=\bm{\varSigma}_{x|s}+o_{p}(1)$.

Next we consider the numerator

$\begin{split}\bm{b}_{T}&=T^{-1/2}\bm{X}_{-p}^{*\prime}\bm{u}-T^{-1/2}\bm{X}_{-p}^{*\prime}\bm{S}\bm{D}_{T,\bm{S}}^{-1}\big{(}\bm{D}_{T,\bm{S}}^{-1}\bm{S}^{\prime}\bm{S}\bm{D}_{T,\bm{S}}^{-1}\big{)}^{-1}\bm{D}_{T,\bm{S}}^{-1}\bm{S}^{\prime}\bm{u}\\\
&=\underbrace{T^{-1/2}\bm{X}_{-p}^{*\prime}\mathcal{M}(\bm{S}_{0})\bm{u}}_{\bm{b}_{T,1}}-\underbrace{T^{-1/2}\bm{X}_{-p}^{*\prime}\mathcal{M}(\bm{S}_{0})\bm{S}_{-0}\bm{D}_{T,\bm{S}_{-0}}^{-1}}_{\bm{R}_{T,2}}\\\
&\quad\times\big{(}\underbrace{\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\mathcal{M}(\bm{S}_{0})\bm{S}_{-0}\bm{D}_{T,\bm{S}_{-0}}^{-1}}_{\bm{R}_{T,1}}\big{)}^{-1}\underbrace{\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\mathcal{M}(\bm{S}_{0})\bm{u}}_{\bm{r}_{T,3}}\\\
\end{split}$

The terms $\bm{R}_{T,1}$ and $\bm{R}_{T,2}$ were already treated for the
denominator. For $\bm{r}_{T,3}$ we get

$\displaystyle{\left\lVert\bm{r}_{T,3}\right\rVert}_{2}$
$\displaystyle\leq{\left\lVert\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\bm{u}\right\rVert}_{2}+{\left\lVert\bm{D}_{T,\bm{S}_{-0}}^{-1}\bm{S}_{-0}^{\prime}\bm{S}_{0}\bm{D}_{T,\bm{S}_{0}}^{-1}\right\rVert}_{2}$
$\displaystyle\quad\times{\left\lVert\bm{D}_{T,\bm{S}_{0}}^{-1}\bm{S}_{0}^{\prime}\bm{u}\right\rVert}_{2}/\lambda_{\min}(\bm{D}_{T,\bm{S}_{0}}^{-1}\bm{S}_{0}^{\prime}\bm{S}_{0}\bm{D}_{T,\bm{S}_{0}}^{-1}\big{)}$
$\displaystyle\leq\gamma_{T,u,1}+\gamma_{T,u,2}+(\phi_{T,01}+\phi_{T,02})\gamma_{T,u,0}/\kappa_{T,0}\leq\gamma_{T,u}.$

Finally, it directly follows from Assumption 3 that

${\left\lVert\bm{b}_{T,1}-T^{-1/2}\bm{X}_{-p}^{*\prime}\bm{u}-\bm{\varSigma}_{wx}^{\prime}\bm{\varSigma}_{ww}^{-1}\bm{S}_{0}^{\prime}\bm{u}\right\rVert}_{2}\leq\delta_{T}.$

Then we can invoke the final statement of Assumption 3 to conclude that

$\bm{b}_{T,1}=T^{-1/2}\begin{bmatrix}\bm{I}_{p}&-\bm{\varSigma}_{sx}^{\prime}\bm{\varSigma}_{ss}^{-1}\end{bmatrix}\begin{bmatrix}\bm{X}_{-p}^{*\prime}\bm{u}\\\
\bm{S}_{0}^{\prime}\bm{u}\end{bmatrix}+o_{p}(1)\xrightarrow{d}N(\bm{0},\sigma^{2}\bm{\varSigma}_{x|s}).$

Plugging everything back into (A.11), it then directly follows from the
results above that

$\begin{split}\sqrt{T}(\hat{\bm{\beta}}_{p}^{*}-\bm{\beta}_{p}^{*})&\xrightarrow{d}N\left(\bm{0},\sigma_{u}^{2}\bm{\varSigma}_{x|s}^{-1}\right).\end{split}$

We finally consider the test statistics. For the Wald test we get that

$\begin{split}\text{Wald}^{*}&=\hat{\sigma}_{u}^{*-2}\underbrace{T^{-1/2}\bm{y}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}}^{*})\bm{X}_{-p}^{*}}_{\bm{F}_{T}^{\prime}}\big{(}\underbrace{T^{-1}\bm{X}_{-p}^{*\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}}^{*})\bm{X}_{-p}^{*}\big{)}^{-1}}_{\bm{G}_{T}}\underbrace{T^{-1/2}\bm{X}_{-p}^{*\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}}^{*})\bm{y}}_{\bm{F}_{T}},\end{split}$

where
$\hat{\sigma}_{u}^{*2}=\hat{\bm{u}}_{+}^{*\prime}\hat{\bm{u}}_{+}^{*}/T=\hat{\bm{u}}_{+}^{\prime}\hat{\bm{u}}_{+}/T$.
It follows directly from the results obtained in the proof of Theorem 1 as
well as the current proof that $\bm{G}_{T}=\bm{\varSigma}_{x|w}+o_{p}(1)$ and
$\bm{F}_{T}\xrightarrow{d}N(\bm{0},\sigma^{2}\bm{\varSigma}_{x|s})$.
Furthermore,

$\begin{split}\hat{\sigma}_{u}^{*2}&=T^{-1}\bm{u}^{\prime}\bm{u}-T^{-1}\bm{u}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{X}_{-p}(\hat{\bm{\beta}}_{p}-\bm{\beta})+T^{-1}\bm{u}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{V}\bm{\delta}\\\
&\quad-T^{-1}(\hat{\bm{\beta}}_{p}-\bm{\beta})^{\prime}\bm{X}_{-p}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{u}+T^{-1}(\hat{\bm{\beta}}_{p}-\bm{\beta})^{\prime}\bm{X}_{-p}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{X}_{-p}(\hat{\bm{\beta}}_{p}-\bm{\beta})\\\
&\quad+T^{-1}\bm{u}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{V}\bm{\delta}-T^{-1}\bm{\delta}^{\prime}\bm{V}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{u}\\\
&\quad+T^{-1}\bm{\delta}^{\prime}\bm{V}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{X}_{-p}(\hat{\bm{\beta}}_{p}-\bm{\beta})+T^{-1}\bm{\delta}^{\prime}\bm{V}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{V}\bm{\delta}.\end{split}$
(A.14)

All terms – except for the leading $T^{-1}\bm{u}^{\prime}\bm{u}$ term – can be
treated exactly as in the proof of Theorem 1 to show that we may replace
$\hat{\mathcal{S}}$ with $\mathcal{S}$, and as in the current proof to show
their asymptotic negligibility. It then follows directly from the weak law of
large numbers that $\hat{\sigma}_{u}^{*2}\xrightarrow{p}\sigma_{u}^{2}$. It
then immediately follows that $\text{Wald}\to\chi^{2}_{p}$.

The result for the LM test follows identically, except that we replace
$\hat{\sigma}_{u}^{*2}$ by
$\tilde{\sigma}_{u}^{*2}=\bm{y}^{\prime}\mathcal{M}(\bm{V}_{+,\hat{\mathcal{S}}})\bm{y}$.
As $\bm{\beta}=\bm{0}$ under the null hypothesis, we can decompose
$\tilde{\sigma}_{u}^{*2}$ as in (A.14) with the terms involving set to 0. It
then directly follows that
$\tilde{\sigma}_{u}^{*2}\xrightarrow{p}\sigma_{u}^{2}$ and
$\text{LM}\xrightarrow{d}\chi^{2}_{p}$. ∎

##  Appendix B Verification of Assumption 3

In this section we show how Assumption 3 is satisfied using standard results
from the unit root and cointegration literature. We will work under a set of
assumptions commonly encountered in the cointegration literature (Johansen,
1995, see e.g.).111111For simplicity we omit any deterministic components such
as intercepts or linear trend terms in the analysis below. With suitable
modifications however, the analysis goes through in a similar way if they are
included.

######  Assumption B.1.

Let $\bm{z}_{t}$ denote an $n$-dimensional process, generated as
$\bm{A}(L)\bm{z}_{t}=\bm{u}_{t}$, where
$\bm{A}(z)=\bm{I}-\sum_{j=1}^{p}\bm{A}_{j}z^{j}$. Then $\bm{A}(L)$ satisfies
the following conditions:

  * (a)

${\left\lvert\bm{A}(z)\right\rvert}=0$ implies $|z|>1$ or $z=1$;

  * (b)

Let $\bm{\varPi}=-\bm{A}(1)$, then we can write
$\Pi=\bm{\mathcal{A}}\bm{B}^{\prime}$, for some $\bm{\mathcal{A}}$ and
$\bm{B}$, where $\bm{\mathcal{A}}$ and $\bm{B}$ are $n\times r$ matrices of
rank $r$.

In addition, for an $n\times r$ matrix $\bm{C}$ let $\bm{C}_{\perp}$ be an
$n\times(n-r)$ matrix of rank $n-r$ such that
$\bm{C}^{\prime}\bm{C}_{\perp}=\bm{0}$, and define
$\bar{\bm{C}}_{\perp}=\bm{C}_{\perp}\left(\bm{C}_{\perp}^{\prime}\bm{C}_{\perp}\right)^{-1}$.
Then, one of the following conditions holds:

  * (c)

$\bm{\mathcal{A}}_{\perp}^{\prime}\tilde{\bm{A}}(1)\bm{B}_{\perp}$ is
nonsingular, where
$\tilde{\bm{A}}(z)=\bm{I}_{n}-\sum_{j=1}^{p-1}\tilde{\bm{A}}_{j}z^{j}$ and
$\tilde{\bm{A}}_{j}=\sum_{i=j+1}^{p}\bm{A}_{i}$;

  * (c′)

$\bar{\bm{\mathcal{A}}}_{\perp}^{\prime}\tilde{\bm{A}}(1)\bar{\bm{B}}_{\perp}=\bm{F}\bm{G}^{\prime}$,
where and $\bm{F}$ and $\bm{G}$ are $(n-r)\times r^{\prime}$ matrices of rank
$r^{\prime}$ ($0\leq r^{\prime}<n-r$).

The rationale for these assumptions comes from rewriting the VAR model as a
vector error correction model (VECM), such that

$\Delta\bm{z}_{t}=\bm{\varPi}\bm{z}_{t-1}+\sum_{j=1}^{p-1}\tilde{\bm{A}}_{j}\Delta
z_{t-j}++\bm{u}_{t},$ (B.1)

where $\tilde{\bm{A}}_{j}=\sum_{i=j+1}^{p}\bm{A}_{i}$ and
$\bm{\varPi}=\bm{A}(1)$. Assumption B.1(a) rules out explosive processes but
allows for unit roots; B.1(b) allows for a reduced rank in $\bm{\varPi}$,
i.e., cointegration is present. If B.1(c) holds, $\bm{z}_{t}$ is $I(1)$ and
linear combinations $\bm{B}^{\prime}\bm{z}_{t}$ are $I(0)$. If B.1(c’) holds,
$\bm{z}_{t}$ is $I(2)$ and appropriate linear combinations are $I(1)$ or even
$I(0)$. This can be seen by again rewriting the VECM as

$\Delta^{2}\bm{z}_{t}=\bm{\varPi}\bm{z}_{t-1}+\bm{\varPhi}\Delta\bm{z}_{t-1}+\sum_{j=1}^{p-2}\bm{A}^{*}_{j}\Delta^{2}\bm{z}_{t-j}+\bm{u}_{t},$
(B.2)

where $\bm{\varPhi}=\tilde{\bm{A}}(1)$ and
$\bm{A}_{j}^{*}=\sum_{i=j+1}^{p-1}\tilde{\bm{A}}_{i}$.121212Formally, we need
another assumption to rule out that the process is $I(3)$, see Johansen (1992,
Theorem 3) for details.

Under these assumptions, we can find appropriate rotation matrices $\bm{Q}$
such that $\bm{\zeta}_{t}=\bm{Q}\bm{z}_{t}$ can be partitioned into variables
“full-rank” $\bm{\zeta}_{0,t},\bm{\zeta}_{1,t}$ and $\bm{\zeta}_{2,t}$ that
are $I(0)$, $I(1)$ and $I(2)$ respectively (where some sets may be empty).
Limit results for appropriate products of variables with different orders of
integration are well-known and have been derived by a variety of papers –
often under weaker assumptions than Assumptions 1 and B.1 – see for example
Phillips and Durlauf (1986), Park and Phillips (1988), Sims et al. (1990),
Phillips and Solo (1992), or Hamilton (1994) for an overview. We therefore
state these results here without proof.

######  Lemma B.1.

Let $\bm{z}_{t}$ satisfy Assumptions 1 and B.1 and let
$\bm{\zeta}_{t}=\bm{Q}\bm{z}_{t}$. Let
$\tilde{\bm{\zeta}}_{t}=(\bm{\zeta}_{0,t},\Delta\bm{\zeta}_{1,t},\Delta^{3}\bm{\zeta}_{2,t})$
and define
$\bm{\varSigma}:=\mathbb{E}\tilde{\bm{\zeta}}_{t}\tilde{\bm{\zeta}}_{t}^{\prime}$,
and partition it conformably such that
$\bm{\varSigma}_{ij}=\mathbb{E}\tilde{\bm{\zeta}}_{i,t}\tilde{\bm{\zeta}}_{j,t}^{\prime}$.
Let
$\bm{\varLambda}=\sum_{j=1}^{\infty}\mathbb{E}\tilde{\bm{\zeta}}_{t}\tilde{\bm{\zeta}}_{t+j}^{\prime}$
and $\bm{\varOmega}=\bm{\varSigma}+\bm{\varLambda}+\bm{\varLambda}^{\prime}$,
where $\bm{\varOmega}_{ij}$ refers to the same submatrix as
$\bm{\varSigma}_{ij}$, and assume that both $\bm{\varSigma}$ and
$\bm{\varOmega}$ are full rank. Let $\bm{B}(\cdot)$ denote an $n$-dimensional
Brownian motion with covariance matrix $\bm{\varOmega}$, which can be
partitioned conformably with $\bm{\zeta}_{0,t},\bm{\zeta}_{1,t}$ and
$\bm{\zeta}_{2,t}$ as
$\bm{B}(\cdot)=(\bm{B}_{0}(\cdot)^{\prime},\bm{B}_{1}(\cdot)^{\prime},\bm{B}_{2}(\cdot)^{\prime})^{\prime}$,
$\bar{\bm{B}}_{i}(s)=\int_{0}^{s}\bm{B}_{i}(u)du$ and $B_{u}(\cdot)$ a
univariate Brownian motion independent of $\bm{B}(\cdot)$. Then we have the
following results (where convergence applies jointly if needed):

$\displaystyle(a)\quad
T^{-1}\sum_{t=1}^{T}\bm{\zeta}_{0,t}\bm{\zeta}_{0,t}^{\prime}\overset{p}{\to}\bm{\varSigma}_{00};$
$\displaystyle(b)\quad
T^{-1/2}\sum_{t=1}^{T}\bm{\zeta}_{0,t-1}u_{t}\overset{d}{\to}N(\bm{0},\sigma_{u}^{2}\bm{\varSigma}_{00});$
$\displaystyle(c)\quad
T^{-1}\sum_{t=1}^{T}\bm{\zeta}_{1,t-1}u_{t}\overset{d}{\to}\int_{0}^{1}\bm{B}_{1}(s)^{\prime}d\bm{B}_{u}(s);$
$\displaystyle(d)\quad
T^{-1}\sum_{t=1}^{T}\bm{\zeta}_{1,t}\bm{\zeta}_{0,t}^{\prime}\overset{d}{\to}\int_{0}^{1}\bm{B}_{1}(s)d\bm{B}_{0}(s)^{\prime}+\bm{\varSigma}_{10}+\bm{\varLambda}_{10};$
$\displaystyle(e)\quad
T^{-2}\sum_{t=1}^{T}\bm{\zeta}_{1,t}\bm{\zeta}_{1,t}^{\prime}\overset{d}{\to}\int_{0}^{1}\bm{B}_{1}(s)\bm{B}_{1}(s)^{\prime}ds;$
$\displaystyle(f)\quad