Uniform Inference for Cointegrated Vector Autoregressive Processes

Uniform Inference For Cointegrated Vector Autoregressive Processes

Introduction

Persistence, i.e., long term sensitivity to small shocks, appears to be a commonly occurring characteristic of many stochastic systems encountered in practice. Such processes are often modeled with cointegration where the persistence can be attributed to a number of shared stochastic trends (random walks). Due to their relative simplicity, cointegration models are widely applied. Cointegration in vector-valued autoregressive processes arises when the characteristic polynomial possesses at least one unit root. The asymptotic theory and, hence, inference is heavily reliant upon the fact that a fixed number of roots can be assumed exactly one and the rest stay sufficiently far away from one. This is the case even for simple regression methods such as ordinary least squares. In practice, such assumptions are overly restrictive and more flexible models are often needed for a better description of the empirical data. Slight deviations from the unit root assumption can severely deteriorate the results of the statistical analysis elliott1998robustness. Thus, the need arises for inference methods that are uniformly valid over a range of stationary and non-stationary behaviours.

So far, methods of inference with proven uniform guarantees only work in the univariate case. Bootstrap inference algorithms are presented in andrews1993exactly, hansen1999grid with uniform guarantees given in mikusheva2007uniform. Furthermore, mikusheva2007uniform provides a uniform asymptotic framework for one-dimensional autoregressive processes with potential unit roots, which served as an inspiration for much of the work in this paper.

The problem is well-understood in one dimension, while less progress has been made in multiple dimensions. The only methods with uniform guarantees known to the authors employ lag augmentation dolado1996making, toda1995statistical. The idea is simple and easy to apply, but lacks efficiency since it essentially involves overfitting the model. Other methods, seeking to avoid this problem, impose restrictive assumptions on the process such that the inference is no longer uniform, but only holds for specific configurations of parameters. The main approach is to assume that the roots are all of a similar proximity to one. In particular, the autoregressive matrix is modeled as a sequence of matrices that approach the identity matrix at some given rate \(k_n\). This is the setup in the instrumental variable methodology (IVX) developed in kostakis2015robust, magdalinos2020econometric, phillips2009econometric. This offers some flexibility in terms of how close the roots can be to one. However, the framework does not allow for processes with simultaneously different degrees of persistence.

Generic results on uniform convergence are provided in andrews2020generic. However, they hold only insofar as one can establish the right asymptotic distributions under appropriate drifting sequences. While progress in this direction has been made phillips2013predictive, phillips2015limit, there has thus far been no general answer to this problem. For example, these papers do not allow the process to have parts that remain stationary for all sample sizes, and only diagonal regression matrices are allowed. Furthermore, the theory in andrews2020generic does not provide a way to perform inference once the appropriate limits are established. The method they propose for the one-dimensional autoregressive process is not feasible for multiple dimensions.

A lot of focus has been given to predictive regression problems in which the predictive power of the past of one process on the future of another is assessed. Efficient tests are developed in campbell2006efficient, jansson2006optimal and, based on the ideas of uniform inference for univariate autoregressive processes, phillips2014confidence leverage these methods for uniformly valid inference. Unfortunately, most of the work only covers the bivariate case in which the regressor is a univariate autoregressive process so that the theory for one dimension directly applies. Another branch of research concerns inference on cointegrating relations robust to deviations from the unit root assumption duffy2023cointegration, franchi2017improved, but, again, specific parameterizations of the deviations limit the generality of these results.

Extending the theory in mikusheva2007uniform to multiple dimensions runs into several difficulties. Firstly, it is not clear what assumptions to put on the autoregressive matrix to ensure that the asymptotic results hold uniformly while still covering all relevant cases. In one dimension, the autoregressive parameter is a scalar and it is sufficient that it is real, bounded in norm by 1, and bounded away from -1 by some small \(\delta\). We need to extend this idea to matrices. Secondly, while many of the results on the asymptotics of the sample covariance matrices generalize nicely to multiple dimensions, the proofs are more involved. For example, mikusheva2007uniform uses Skorohod's embedding, which famously only works in one dimension, to prove that the errors can be assumed to be Gaussian. The third and perhaps most profound difficulty is that the multivariate setting allows for cointegrated systems (or almost cointegrated systems in the case where the roots are only close to unity). This gives rise to certain asymptotic discontinuities. In particular, it necessitates a proper normalization of the sample covariances. These problems extend to the inferential side where, additionally, computational challenges arise. Naively adapting, for example, the grid bootstrap approach of hansen1999grid would cause the computational complexity to explode.

The present work provides methods of inference for multivariate vector autoregressive processes that are proven to be valid uniformly over a set of parameters given in Assumption (ref), including processes that are cointegrated and with roots arbitrarily close to the unit circle. Assumptions (ref) and (ref) are probably slightly stronger than necessary, but they are easy to work with and provide a clear interpretation while still being much more general than anything we have been able to find in the literature. Most of the work dealing with similar problems assume either that \(\Gamma\) is a drifting sequence of diagonal matrices (see, e.g., phillips2013predictive, phillips2015limit, magdalinos2020econometric) or impose structural assumptions such that the problem can be reduced to studying only diagonal \(\Gamma\) (see, e.g., duffy2023cointegration where they assume that \(\Gamma\) is normal).

The main contributions are the following. First, we extend the asymptotic results of mikusheva2007uniform to vector autoregressive processes and where uniformity also holds over a family of martingale difference error processes. The main result is Theorem (ref), stating that the asymptotic distributions of the relevant sample covariances can be approximated by stochastic integrals of Ornstein-Uhlenbeck processes (a direct analog to the univariate case). Both the result and the proof are interesting in their own rights. As part of the proof we show that one can approximate the finite sample distribution of crucial statistics by replacing the general error terms with Gaussian errors. We can sample from this approximation at a comparably low computational cost which facilitates inference greatly.

Second, we provide several ways to construct confidence regions for the autoregressive parameter and show that these are uniformly valid. Confidence regions constructed with IVX and lag augmentation are uniformly valid (but IVX only for the entire autoregressive matrix).

Third, we show how these confidence regions can be used to answer more general inference questions. The two main applications are confidence intervals for a single coordinate and predictive regression testing with a multivariate regressor. To the best of our knowledge, there have thus far been no attempt in the literature to deal with these applications in a uniform fashion (apart from lag augmentation). We run Monte Carlo experiments to verify the theoretical results and compare the finite sample properties of the different methods.

Our last contribution is the development of efficient algorithms to solve these inferential tasks. We show how the Evaluation-Approximation-Maximization (EAM) algorithm from kaido2019confidence circumvents the exploding computational cost inherent in algorithms relying on grid-like methods. Combined with the Gaussian approximation results, this is what makes inference possible in our two main applications.

The paper is structured as follows. Section 2 introduces notation and relevant concepts. In particular, it explains the concept of uniform convergence of random variables and presents vector autoregressive processes. Section 3 is devoted to presenting and proving the main asymptotic results. Section 4 deals with inference and shows how the results of Section 3 can be applied to obtain uniformly valid confidence regions. Furthermore, it contains a section on predictive regression, lag augmentation, and IVX. Section 5 contains the results of our Monte Carlo experiments. Finally, Section 6 concludes. The Appendix contains proofs and further technical details on martingale limit results, the Gaussian approximation, the simulation experiments as well as the EAM algorithm, and details on lag augmentation and IVX.

Preliminaries

This paper is concerned with vector autoregressive processes of order 1 (VAR(1) processes) fulfilling certain assumptions. Specifically, processes that may be integrated of order 1 and cointegrated, that is, processes for which the first difference is stationary and there exists some linear combinations of the coordinate processes that are stationary.

Notation: For a matrix \(A\in\mathbb{C}^{d\times d}\), \(A^T\) denotes its conjugate transpose and its trace is \(\text{tr}(A)\). We write \(\sigma_{max}(A)\) (\(\sigma_{min}(A)\)) for the largest (smallest) singular value of \(A\), and \(\lambda_{max}(A)\) (\(\lambda_{min}(A)\)) for the eigenvalue of \(A\) with the largest (smallest) magnitude. \(||A||\) is the Frobenius norm and \(||A||_2\) is the spectral norm, i.e., \(||A||=\sqrt{\text{tr}(A^T A)}\) and \(||A||_2 = \sqrt{\sigma_{max}(A)}\). For vectors, \(||\cdot||\) is the usual Euclidean norm. Define \(S_{d}\) to be the set of \(d\times d\) positive semidefinite matrices. We employ the usual big-\(O\) and little-\(o\) notation and use \(o_p\) to denote convergence in probability and \(O_p\) to denote boundedness in probability.

Uniform convergence of random variables

The definitions of uniform convergence in probability and in distribution are essentially the same as in kasy2019uniformity, lundborg2021conditional. Assume some background probability space, \((\Omega, \mathcal{F}, \mathbb{P})\), on which all future random variables are defined. For two random vectors, \(X\) and \(Y\), taking values in \((\mathbb{C}^d, \mathcal{B}(\mathbb{C}^d))\), we denote by \(P_{X}\) and \(P_{Y}\) the law of \(X\) and \(Y\) and write \(X \stackrel{\mathcal{L}}{=} Y\) if they are equal in law. Let \(BL_1\) be the space of functions \(f:\mathbb{C}^d\rightarrow [-1, 1]\) that are Lipschitz continuous with constant at most 1. Let \(\mathcal{P}(\mathbb{C}^d, \mathcal{B}(\mathbb{C}^d))\) be the set of probability measures on \((\mathbb{C}^d, \mathcal{B}(\mathbb{C}^d))\). The bounded Lipschitz metric on \(\mathcal{P}(\mathbb{C}^d, \mathcal{B}(\mathbb{C}^d))\) is given by \[ d_{BL}(\mu, \nu):=\sup_{f\in BL_1}\left|\int_{\mathbb{C}^d} f d\mu - \int_{\mathbb{C}^d} f d\nu\right|, \quad \mu, \nu\in \mathcal{P}(\mathbb{C}^d, \mathcal{B}(\mathbb{C}^d)). \] We use the shorthand \(d_{BL}(X, Y) = d_{BL}(P_X, P_Y)\) to denote the bounded Lipschitz metric between the laws of two random variables, \(X\) and \(Y\). It is well known that \(d_{BL}\) metrizes weak convergence which motivates the following definition of uniform convergence.

Uniform convergence could also be stated as convergence along all sub-sequences \(\theta_n\subset\Theta\) (see Definition 2 and Lemma 1 in kasy2019uniformity). Additionally, we allow the limiting distribution to be a sequence, since the results below are stated in terms of an approximating sequence of random variables. We obtain the conventional notion by letting \(Y_{n, \theta} = Y_\theta\).

Model

Consider some \(\Theta \subset \mathbb{R}^{d\times d}\times S_d\times\mathbb{R}_+\). For any \(\theta\in\Theta\) there exist \(\Gamma_\theta, \Sigma_\theta\), and \(c_\theta\) such that \(\theta = (\Gamma_\theta, \Sigma_\theta, c_\theta)\). Let \(N_\theta\in\{1,...,d\}\) denote the number of distinct eigenvalues of \(\Gamma_\theta\) and \(\lambda_\theta\in\mathbb{C}^{N_\theta}\) the corresponding vector of ordered eigenvalues, that is, \(|\lambda_{\theta, 1}|\ge |\lambda_{\theta, 2}| \ge ... \ge |\lambda_{\theta, N_\theta}|\) with multiplicities \(m_{\theta, 1},..., m_{\theta, N}\in\{1,..., d\}\). Where this does not cause confusion, we omit the subscript \(\theta\). Let \((X_{t, \theta})_{t\in\mathbb{N}, \theta\in\Theta}\) and \((\epsilon_{t, \theta})_{t\in\mathbb{N}, \theta\in\Theta}\) be two families of \(\mathbb{R}^d\)-valued stochastic processes and \((\mathcal{F}_{t,\theta})_{t\in\mathbb{N}}\) the filtration generated by \((\epsilon_{t, \theta})_{t\in\mathbb{N}}\).

These assumptions ensure that \(X_{t, \theta}\) is a VAR(1) process started at \(0\) with model parameters given by the index \(\theta\). Assumptions (ref) and (ref) imply that \(X_{t, \theta}\) is adapted to \(\mathcal{F}_{t, \theta}\).

For future reference let us define the following set of \(d\times d\) matrices. For a given \(\delta > 0\), let \(\mathcal{J}_d(\delta)\subset\mathbb{C}^{d\times d}\) be the set of upper triangular matrices such that every \(J\in\mathcal{J}_d(\delta)\) can be decomposed as \(J=D+N\) with \(D\) diagonal such that \(|D_{11}|\ge |D_{22}|\ge\dots\ge |D_{dd}|\) and \(N\) equal to 0 everywhere except on the super-diagonal where it satisfies \(N_{i, i+1}\in\{0, 1\}\) if \(|D_{ii}|\le \delta\) and 0 otherwise for \(i=1,\dots, d-1\). In other words, every \(J\in\mathcal{J}_d(\delta)\) can be written as a block diagonal matrix where the upper left block is diagonal and contains all eigenvalues greater than \(\delta\) while the lower right block can have ones on the super-diagonal and has eigenvalues less than \(\delta\). We call the matrices in \(\mathcal{J}_d(\delta)\) Jordan-like.

Asymptotic Properties

The key building blocks for inference are the two covariance matrices

Obviously, \(S_{XX}\) and \(S_{X\epsilon}\) are families of stochastic processes depending on \(n\) and \(\theta\), which we suppress to avoid cluttering up the notation, but the dependence should be kept in mind. We first need to determine what happens to \(S_{XX}\) and \(S_{X\epsilon}\) when \(n\) goes to infinity and for varying \(\theta\). We cannot hope to say anything uniformly without further assumptions on \(\Theta\). The following assumptions are sufficiently general to cover a wide range of behaviours while still allowing for uniform asymptotic results.

Assumption (ref) is a moment condition on the error process which is needed for some of the triangular array martingale difference limit results. It is implied by the other conditions if the errors are i.i.d. Gaussian. Assumption (ref) states that \(||\Sigma_\theta||\) and \(||\Sigma_\theta^{-1}||\) are uniformly bounded for any matrix norm. In particular, \(\Sigma_\theta\) is of full rank. This is a natural condition when considering uniform convergence. The important assumptions are (ref) and (ref). Both assumptions have clear interpretations and are sufficient if we want to limit our attention to processes that are at most integrated of order 1 and without seasonal cointegration. In particular, to avoid higher orders of integration, we must restrict all eigenvalues to have magnitude less than 1 (ensured by Assumption (ref), see Fig. (ref)). For eigenvalues with magnitude 1, the corresponding Jordan block must be scalar archontakis1998alternative (ensured by Assumption (ref)). Note that the matrices \(J_\theta\) in Assumption (ref) are not required to be Jordan matrices so that the assumption allows, for example, for matrices of the form \[ \Gamma =

\] for \(\lambda, \lambda' \in \mathbb{R}\) arbitrarily close together as long as \(|\lambda|, |\lambda'|\le 1 - \alpha\).

To avoid seasonal cointegration, eigenvalues with magnitude 1 are restricted to be exactly equal to 1. Specifically, as the magnitude of an eigenvalue approaches one, the eigenvalue itself approaches 1. See Fig. (ref) for the region of the complex plane satisfying Assumption (ref) for some small \(\alpha \in (0, 1)\). Since $\alpha$ can be chosen freely, the parameter space can be made arbitrarily close to the unit circle.

To state our main result, define \(M_i = \sum_{j\le i} m_j\) and write \(i_k = \min\{i \ge 1 | M_i \ge k\}\). The limiting behaviour of the \(k\)'th coordinate of \(X_{t, \theta}\) depends on how close the eigenvalue, \(\lambda_{i_k}\), is to 1. As in the univariate case mikusheva2007uniform, we unify the range of asymptotics with an Ornstein-Uhlenbeck process. For any \(\theta\), we let \(C_n(\theta)\) be the \(d\times d\) diagonal matrix whose \(i\)'th diagonal block is \(n\log(|\lambda_i|)I_{m_i}\) with the convention that \(\log(0)=-\infty\). For convenience, we sometimes suppress the dependence on \(\theta\) and \(n\) and just write \(C\). In what follows, uniform convergence of random matrices means uniform convergence of the vectorization of these matrices so that Definition (ref) applies directly.

We prove the uniform results for the special case of $F_{\theta}$ being the identity, i.e., under the assumption that \(\Gamma\in \mathcal{J}_d(1 - \delta)\). Technically, this allows for complex-valued \(\Gamma\), thus, the proofs are more general than the real-valued case. Also, it is not hard to generalize to any $F_{\theta}$ fulfilling Assumption (ref). Indeed, for \(X_{t, \theta}\) generated by \(\Gamma\in \mathbb{R}^{d\times d}\), there exist \(F\in \mathbb{C}^{d\times d}\) and \(J\in J_d(1 - \delta)\) such that \(F\) is invertible with \(F^{-1}JF = \Gamma\). The transformed process \(\tilde{X}_{t, \theta} = FX_{t, \theta}\) is then of the required form with parameters \(\tilde{\theta}=(J, F\Sigma F^T, ||F||^{2+\delta}c)\). Assuming that \(F\) is uniformly invertible and bounded in norm then ensures that \(\tilde{\theta}\) satisfies Assumption (ref).

We prove Theorem (ref) in several steps. The main idea is to split \(\Theta\) into overlapping regions and prove that Theorem (ref) holds in each region. Consider \[ R_{n, 0} := \left\{\theta\in\Theta : |\lambda_1| \le 1 - \frac{\log n}{n}\right\}, \quad R_{n, d}:= \left\{\theta\in\Theta : |\lambda_N| \ge 1 - n^{-\eta}\right\}, \] where \(\eta\in(0, 1)\) is to be specified later. The two regions correspond to the stationary and local-to-unity (non-stationary) regimes, respectively, in the univariate case. Different asymptotics arise depending on the region and, in particular, on how fast the eigenvalues converge to unity. Throughout the rest of this section we assume that Assumptions (ref) and (ref) hold with \(F_\theta = I\).

Non-stationary asymptotics

In this section we consider sequences of parameters in the non-stationary region \(R_{n, d}\). For simplicity we assume throughout this subsection that \(\epsilon_{t, \theta}\) is i.i.d. Gaussian with mean 0 and covariance \(\Sigma_\theta\). In Appendix (ref) it is argued why this is not a restriction. Indeed, all the relevant sample moments can be approximated by Gaussian counterparts.

Let \((W_{t})_{t\in[0,1]}\) be a standard \(d\)-dimensional Brownian motion. Since \(\Sigma^{\frac{1}{2}}(W_{t/n} - W_{(t-1)/n})\stackrel{\mathcal{L}}{=} \epsilon_{t, \theta}/\sqrt{n}\) for all \(n\in\mathbb{N}\), \(0\le t \le n\) and \(\theta\in\Theta\), we get

where $H$ is defined in (ref) and \(f(t, s, n, \theta) = \sqrt{n}H^{-\frac{1}{2}}\Gamma^{\lfloor nt \rfloor - \lfloor ns \rfloor - 1} \Sigma^{\frac{1}{2}}\mathbf{1}\{s \le \lfloor nt\rfloor / n\}\). We then see that the following Lemma is a direct consequence of Lemma (ref) in Appendix (ref).

Stationary asymptotics

In this section we consider sequences of parameters in the stationary region \(R_{n, 0}\). We first show that the classical asymptotic theory for stationary VAR(1) processes applies. Since \(R_{n, 0}\) and \(R_{n, d}\) overlap, we then prove that the right hand sides of (ref) and (ref) converge to the standard stationary limiting distributions for the diagonal entries \(C_{ii}\) going to \(-\infty\).

The standard stationary theory in multiple dimensions mimics the univariate case. We follow the same strategy as in phillips2007limit, but allowing for multiple dimensions and a family of error processes \(\epsilon_{t, \theta}\). In this regime, we find that, when properly normalized, \(S_{XX}\) converges in probability to the identity matrix and \(\text{vec}(S_{Xe})\) converges in distribution to a \(d^2\)-dimensional standard Gaussian.

For the special case \(s=1\), equation (ref) shows that \(S_{XX}\) converges in probability to the identity matrix. By Theorem (ref) and Proposition 8 in the supplementary material for lundborg2021conditional, the proof of Theorem (ref) in the stationary regime is complete if we can show that, for any \((\theta_n\in R_{n, 0})_{n\in\mathbb{N}}\),

We emphasize that \(G\) and \(C\) in eqs. (ref)-(ref) are functions of \(\theta_n\) and therefore they are sequences of matrices. In particular, \(C_{ii}\le n\log(1-\log(n)/n)\rightarrow -\infty\) for \(n\rightarrow \infty\) and \(1\le i \le d\). Eqs. (ref)-(ref) are therefore a consequence of the following.

Mixed asymptotics

So far, all the eigenvalues were in the same regime. We have yet to explore what happens when there are eigenvalues in both regimes. We call this case the mixed regime. Define for \(1\le k \le d-1\) and some fixed \(\gamma\in (0, 1-\eta)\) the sets \[ R_{n, k} = \left\{\theta\in\Theta : M_{i_k} = k, |\lambda_{i_k}| \ge 1 - n^{-\eta - \gamma}, |\lambda_{i_{k+1}}|\le 1 - n^{-\eta - \gamma}\right\}. \] Since \(1-n^{-\eta} \le 1 - n^{-\eta - \gamma} \le 1 - \log(n)/n\), then for \(\theta\in R_{n, k}\) there are at least \(k\) coordinates in the non-stationary regime and \(d-k\) coordinates in the stationary regime (and some might be in both). Furthermore, for any \(n\in \mathbb{N}\), \(\Theta = \bigcup_{0\le k\le d} R_{n, k}.\) Thus, showing that (ref) and (ref) hold uniformly over \(R_{n, k}\) for any fixed \(1\le k\le d-1\) completes the proof of Theorem (ref). This is the content of the following lemma proved in Appendix (ref).

Uniform Inference

Having established the asymptotic properties of \(S_{XX}\) and \(S_{X\epsilon}\), we now seek to develop uniformly valid methods of inference. We focus on two important cases in which uniformly valid inference has so far proven challenging: predictive regression testing and coordinate confidence intervals. Inference in these settings can be hard even from a point-wise perspective since the presence of exact unit roots makes it problematic to construct test statistics with standard asymptotic distributions.

It is not trivial to conduct inference on \(\Gamma\) even in lieu of Theorem (ref). The main problem is the presence of the nuisance parameter \(C_n(\theta)\) in (ref)-(ref), which cannot be uniformly consistently estimated. Indeed, for a sequence \(\Gamma_n=I - C/n\) where the real part of the eigenvalues of \(C\in \mathbb{R}^{d\times d}\) are all strictly negative, the problem is essentially equivalent to estimating \(C = n(I - \Gamma)\). But it is well known that, in this setting, \(\Gamma\) can only be estimated at rate \(O(n^{-1})\). One way to solve this is by the use of test inversion or so-called grid bootstrap methods, which have been widely applied in the unitary case (see hansen1999grid, mikusheva2007uniform for grid bootstrap and campbell2006efficient, phillips2014confidence for an application to predictive regression). While this is fairly easy in one dimension, adapting these methods to vector autoregressive processes is prohibitive since the computational complexity quickly explodes. We now present an approach, which keeps the computational burden to a minimum. We omit the dependence on \(\theta\) in the subscript of all random variables.

Consider the least squares estimator, \(\hat{\Gamma}\), given by \[ \hat{\Gamma} = \frac{1}{n}\sum_{t=1}^n X_{t}X_{t-1}^T\left(\frac{1}{n}\sum_{t=1}^n X_{t-1}X_{t-1}^T\right)^{-1} = \Gamma + S_{X \epsilon}^T S_{XX}^{-1}. \] It follows from Theorem (ref) that \(\hat{\Gamma}\) is a uniformly consistent estimator of \(\Gamma\) with a rate of convergence depending on the proximity of the eigenvalues of \(\Gamma\) to one. Indeed, since \(\sqrt{n}(\hat{\Gamma} - \Gamma)=(\sqrt{n} S_{X\epsilon})^T S_{XX}^{-1}\), we find that \(n(\hat{\Gamma} - \Gamma)S_{XX}(\hat{\Gamma}-\Gamma)^T = O_p(1)\) which implies that \(\sqrt{n}(\hat{\Gamma} - \Gamma) = O_p(1)\). We define a uniformly consistent estimator of \(\Sigma\) by averaging the squared residuals, i.e., with \(\hat{\epsilon}_t = X_t - \hat{\Gamma}X_{t-1}\) we define \(\hat{\Sigma} = S_{\hat{\epsilon}\hat{\epsilon}}\) analogously to \(S_{XX}\) with \(X_{t-1}\) replaced by \(\hat{\epsilon}_t\). Another consequence of Theorem (ref) is a uniform approximation of the \(t^2\)-statistic. In particular,

uniformly over \(\Theta\) where \(C = C_n(\theta)\) is given in Section (ref) and \(\hat{J}_{t, C}\) is defined analogously to \(J_{t, C}\) but with \(\Sigma\) replaced by the consistent estimator \(\hat{\Sigma}\). Thus, for a fixed significance level, \(\alpha\in (0, 1)\), a uniformly valid \(100(1-\alpha)\%\) confidence region for \(\Gamma\) can be constructed by test-inversion. Let \(q_{n, \Gamma}(1-\alpha)\) denote the \(100(1-\alpha)\%\) quantile\footnote{Note that the quantile depends on \(\Gamma\) (which is also made explicit in the notation) since the asymptotic distribution depends on \(C\) which is a function of \(\Gamma\) and \(n\).} of \(t^2_\Gamma\) and define the confidence region

However, the distribution of \(t^2_\Gamma\) is non-standard and, therefore, computing the quantiles \(q_{n, \Gamma}\) requires extensive simulations and can be quite expensive.

Another approach relies on the Gaussian approximations (Appendix (ref)). It is similar to Andrew's Method (see mikusheva2007uniform), and similar in spirit to grid bootstrap. It was originally suggested by andrews1993exactly but has so far only been applied in the univariate case. For a given \(\Gamma\), define the \(VAR(1)\) process, \((Y_t)_{t\in \mathbb{N}}\), by \[ Y_t = \Gamma Y_{t-1} + e_t, \quad Y_0 = 0 \] where \(e_t\sim \mathcal{N}(0, \hat{\Sigma})\) i.i.d. Let \[ \tilde{t}^2_\Gamma = \text{tr}\left(n\hat{\Sigma}^{-\frac{1}{2}}S_{eY} S_{YY}^{-1} S_{eY}^T\hat{\Sigma}^{-\frac{1}{2}}\right) \] and denote by \(\tilde{q}_{n, \Gamma}(1-\alpha)\) the \(100(1-\alpha)\%\) quantile of \(\tilde{t}^2_\Gamma\). A confidence region for \(\Gamma\) is then obtained by

The distribution of \(\tilde{t}^2_\Gamma\) is still non-standard, but the quantiles \(\tilde{q}_{n, \Gamma}\) are much easier to compute by simulation. The following Theorem states that both confidence regions are uniformly asymptotically valid over the parameter space \(\Theta\). A proof can be found in Appendix (ref).

Predictive Regression

One application is robust inference in the predictive regression model. For a deeper discussion of why uniformly valid inference methods are important in this setting see campbell2006efficient, elliott1994inference, which cover the case of a univariate regressor, but the same considerations hold more generally. To fix ideas, consider \[ \Theta_P = \left\{\theta \in \Theta : \Gamma_{j1}= 0 \; \forall j=1,..., d\right\}. \] If \(X_t\) is a VAR(1)-process satisfying Assumption (ref) parameterized by \(\theta\in \Theta_P\), then we can split \(X_t=(Y_t, \tilde{X}_t^T)^T\) and \(\epsilon_t = (\rho_t, \tilde{\epsilon}_t^T)^T\) into their first coordinate and their last \(d-1\) coordinates such that

where \(\gamma^T = (\Gamma_{1j})_{2\le j\le d}\) and \(\tilde{\Gamma} = (\Gamma_{ij})_{2\le i,j\le d}\). The parameter of interest is \(\gamma\). The standard approach is to compute the least squares estimator \(\hat{\gamma}\) and base inference on the \(t^2\)-statistic. Unfortunately, we encounter the same issues as described above. To see this, let \(\Sigma_Y = \Sigma_{11}\), \(\Sigma_X = (\Sigma_{i, j})_{2\le i,j\le d}\), \(\Sigma_{YX} = (\Sigma_{1, j})_{2\le j\le d}\), and \(\Sigma_{XY}=\Sigma_{YX}^T\) and define \(\delta = \Sigma_{X}^{-1}\Sigma_{XY}\). Then, adopting previous notation, Theorem (ref) yields \[ \hat{t}_\gamma^2 \rightarrow_w \left(\int_0^1 \tilde{J}_{t,C}dB_{1,t}\right)^T\left(\int_0^1 \tilde{J}_{t, C}\tilde{J}_{t, C}^T dt\right)^{-1}\int_0^1 \tilde{J}_{t, C} dB_{1,t} =: t_\gamma^2 \] uniformly over \(\Theta_P\) where \(\tilde{J}_{t, C}\) consists of the last \(d-1\) coordinates of \(\hat{J}_{t, C}\) and \(B_{1, t}\) is the first coordinate and \(B_{2, t}\) the last \(d-1\) coordinates of \(W_t \hat{\Sigma}^{\frac{1}{2}}\). Since \(B'_{1,t} = (\Sigma_Y - \delta^T \Sigma_X\delta)^{-\frac{1}{2}}(B_{1, t} - \delta B_{2, t})\) is a standard \((d-1)\)-dimensional Brownian motion independent of \(B_{2, t}\) satisfying \(B_{1, t} = \delta B_{2, t} + (\Sigma_Y - \delta^T \Sigma_X\delta)^{\frac{1}{2}}B'_{1, t}\), we find that

where \(Z_{\tilde{\Gamma}}=(\int \tilde{J}_{t, C}\tilde{J}_{t, C}^Tdt)^{-\frac{1}{2}}\int \tilde{J}_{t, C}dB_{2, t}\), and \(Z\) is a \((d-1)\)-dimensional standard normal vector independent of \(Z_{\tilde{\Gamma}}\). The nuisance parameter, \(C=C_n(\theta)\), is therefore also present in the distribution of \(t^2_\gamma\) via \(Z_{\tilde{\Gamma}}\) necessitating alternative methods of inference. Using the results of Theorem (ref), we adopt the univariate Bonferroni approach of campbell2006efficient to obtain uniformly asymptotically valid confidence intervals. Say we want to find a confidence region for \(\gamma\) with significance level \(\alpha\in (0, 1)\). For \(\alpha_1, \alpha_2\in (0, 1)\) with \(\alpha_1 + \alpha_2 = \alpha\), the construction proceeds as follows: First construct a \(100(1-\alpha_1)\%\) confidence region for \(\tilde{\Gamma}\) using, e.g., either \(CR(\alpha_1) = CR_a(\alpha_1)\) or \(CR(\alpha_1) = CR_b(\alpha_1)\) (suitably modified for \(d-1\) dimensions). Then, for each \(\tilde{\Gamma}\in CR(\alpha_1)\), let \(CR_{\gamma|\tilde{\Gamma}}(\alpha_2)\) be a \(100(1-\alpha_2)\%\) confidence region for \(\gamma\) given \(\tilde{\Gamma}\). A confidence region not depending on \(\tilde{\Gamma}\) and with coverage of at least \(100(1-\alpha)\%\) is then obtained via a Bonferroni correction:

Let \(\hat{\gamma}_{\tilde{\Gamma}}\) be the estimator obtained by regressing \(Y_{\tilde{\Gamma}, t} = Y_t - \hat{\Sigma}_{YX}\hat{\Sigma}_X^{-1}(\tilde{X}_t - \tilde{\Gamma}\tilde{X}_{t-1})\) on \(\tilde{X}_{t-1}\) with standard error \(\hat{\sigma}^2_Y = \hat{\Sigma}_Y - \hat{\Sigma}_{YX}\hat{\Sigma}_X^{-1}\hat{\Sigma}_{XY}\). A choice for \(CR_{\gamma|\tilde{\Gamma}}(\alpha_2)\) is then given by

with \(q_{d-1, 1-\alpha_2}\) denoting the \(1-\alpha_2\) quantile of the \(\chi^2_{d-1}\) distribution and \(\hat{t}^2_{\gamma|\tilde{\Gamma}}\) the usual \(t^2\)-statistic for the estimator \(\hat{\gamma}_{\tilde{\Gamma}}\) evaluated at \(\gamma\). A proof of the following is given in Appendix (ref).

Lag agumentation

Other approaches have been suggested to be robust against deviations from exact unit root assumptions. The first one is the lag-augmented VAR methodology proposed by dolado1996making, toda1995statistical. In our setup of VAR(1) processes this approach regresses \(X_t\) on \(X_{t-1}\) and the additional augmented lag \(X_{t-2}\) upon which standard inference methodology is valid. Say, for example, that we are interested in testing the hypothesis \(H_0:A\text{ vec}(\Gamma) = b\) for \(A\in\mathbb{R}^{k\times d^2}\) of rank \(k\le d^2\) and \(b\in\mathbb{R}^k\). Let \(\hat{\Pi}_{LA} \in \mathbb{R}^{d\times 2d}\) denote the least squares estimator in the lag-augmented regression of \(X_{t}\) on \(\bar{X}_t = (X_{t-1}^T, X_{t-2}^T)^T\). Denote \(D = (I_d, 0)^T\) and define \(\hat{\Gamma}_{LA} = \hat{\Pi}_{LA}D\) along with the Wald-statistic

where \(\hat{\sigma}^2_{LA, A} = A\hat{\Sigma}_{LA}A^T\) and \(\hat{\Sigma}_{LA} = \hat{\Sigma}^{-1}\otimes \hat{\Sigma}\). Then, under the null, \(\hat{t}^2_{LA, A, b}\) converges in distribution to \(\chi^2_{k}\) uniformly over \(\Theta\) allowing for construction of uniformly valid confidence intervals and tests (see Appendix (ref)). For example, a \((1-\alpha)100\%\) confidence interval for \(\Gamma_{i j}\) is given by

where \(z_{1-\alpha/2}\) is the \(1-\alpha/2\) standard normal quantile.

IVX

Another approach, known as IVX, that deals specifically with the potential presence of unit roots uses endogenously constructed instrumental variables to slow down the rate of convergence of the estimator enough to ensure mixed Gaussian limiting distributions. It was first suggested by phillips2009econometric and later extended in kostakis2015robust, magdalinos2020econometric. The most general framework considered so far was proposed in magdalinos2020econometric. However, they make the crucial assumption that all roots converge to unity at the same speed. While this allows for easy construction of confidence intervals of general linear functions of \(\Gamma\) and simplifies the theory somewhat, this is a significant restriction. In particular, it does not yield uniform guarantees as the ones discussed in this paper. This excludes, for example, the mixed regime discussed in Section (ref) covering cases where parts of the process are stationary and others exhibit random walk behaviour. In this section we detail how one may achieve truly uniform results. This comes at the cost of less general confidence regions, which is essentially because we need to employ different normalizations depending on how close the different roots are to unity. This is akin to using \(\hat{t}^2_{\Gamma}\) for inference.

The idea of IVX is simple. We achieve Gaussian asymptotics by constructing an endogenously generated instrument that lies in the stationary regime and then perform IV-regression. For some fixed \(\beta\in(1/2, 1)\), we define, for each \(n\in\mathbb{N}\), the instrument \((Z_t)_{t\in\mathbb{N}}\) by \[ Z_t = (1-n^{-\beta})Z_{t-1} + \Delta X_t, \quad Z_0 = 0, \] where we have suppressed the dependence on \(n\) in the notation. For each \(n\in\mathbb{N}\), \(Z_t\) is a VAR(1) process with the error terms given by \(\Delta X_t\) and the sequence of coefficients, \((1-n^{-\beta})I\), fall inside the stationary regime. The IVX estimator is the IV estimator of regressing \(X_t\) on \(X_{t-1}\) and the instrumental variable \(Z_{t-1}\), i.e., \[ \hat{\Gamma}_{IV} = \sum_{t=1}^n X_tZ_{t-1}^T\left(\sum_{t=1}^n X_{t-1}Z_{t-1}^T\right)^{-1}. \] The corresponding \(t^2\)-statistic for testing the null \(H_0:\Gamma = \Gamma_0\) is then given by \[ \hat{t}^2_{IV, \Gamma_0} = \text{tr}\left(n\hat{\Sigma}^{-\frac{1}{2}}\left(\hat{\Gamma}_{IV} - \Gamma_0\right)S_{XZ}S_{ZZ}^{-1}S_{ZX}\left(\hat{\Gamma}_{IV} - \Gamma_0\right)^T\hat{\Sigma}^{-\frac{1}{2}}\right), \] where \(S_{XZ}\) and \(S_{ZZ}\) are defined analogously to \(S_{XX}\). It turns out that inference based on this statistic is standard. In particular, \(\hat{t}^2_{IV, \Gamma}\) has the standard asymptotic \(\chi^2_{d^2}\) distribution uniformly over the parameter space \(\Theta\) (see Appendix (ref) for more details) and therefore, for fixed \(\alpha\in(0, 1)\), a confidence region with asymptotic uniform level is given by \[ CR_{IV}(\alpha) = \left\{\Gamma : \hat{t}^2_{IV, \Gamma}\le q_{d^2, 1-\alpha}\right\}. \]

Simulations

In this section we investigate the finite sample properties of the methods described in the preceding section. First, we consider the problem of constructing a confidence interval for \(\Gamma_{ij}\). Testing whether \(X_{j, t}\) Granger causes \(X_{i, t}\), which in this case amounts to testing the null \(H_0: \Gamma_{ij} = 0\), is then equivalent to checking if 0 is contained in the confidence interval. Since there is nothing special about the choice of \(i\) and \(j\) we choose to focus on \(\Gamma_{11}\) for simplicity. The second problem we consider is that of testing \(H_0:\gamma = 0\) in the predictive regression model.

Confidence intervals

Throughout we fix the significance level at \(\alpha = 0.05\). We compare three different ways of constructing confidence intervals for \(\Gamma_{11}\). The first is lag-augmentation yielding the confidence interval \(CI_{LA}(\alpha)=CI_{LA, 11}(\alpha)\), equation (ref). The other two methods first compute a confidence region for the entire matrix \(\Gamma\) and then find a confidence interval for \(\Gamma_{11}\) by projecting the confidence region onto the first coordinate. That is, for a given confidence region of \(\Gamma\) with \((1 - \alpha)100\%\) coverage, \(CR(\alpha)\), we obtain the projected confidence interval \[ CI(\alpha)= \left(\inf_{\Gamma\in CR(\alpha)}\Gamma_{11}, \sup_{\Gamma\in CR(\alpha)}\Gamma_{11}\right). \] We let \(CI_b(\alpha)\) (respectively \(CI_{IV}(\alpha)\)) be the confidence interval obtained by projecting \(CR_b(\alpha)\) (respectively \(CR_{IV}(\alpha)\)). Of these, \(CI_b\) is by far the most costly to compute as the dimension increases. The constraint in the optimization problem is costly to evaluate due to the need for simulations to compute the critical value \(\tilde{q}_{n, \Gamma}(1-\alpha)\) at each \(\Gamma\). This issue can, however, be partly resolved by the use of the EAM-algorithm kaido2019confidence. See Appendix (ref) for details. \(CI_{IV}\) also involves an optimization problem, but the constraint is cheap to evaluate since the critical value in \(CR_{IV}\) is fixed and standard. \(CI_{IV}\) can be computed by standard solvers. We let \(\beta=0.9\) in the IVX regression.

To verify that the confidence intervals are truly uniform, we look at choices of \(\Gamma\) with eigenvalues in different regimes. In particular, for a fixed dimension \(d\) and sample size \(n\), we consider \(\Gamma\in\mathbb{R}^{d\times d}\) with eigenvalues \(\lambda_1(\Gamma) = 1\) and \(\lambda_i(\Gamma) = 1 - (1/n)^{1/(i-1)}\) for \(i=2,...,d\). For each simulation, we draw a new set of random eigenvectors and the errors are i.i.d. Gaussian with non-diagonal covariance matrix. For a detailed explanation of the setup, see Appendix (ref). The results are recorded in Table (ref).

All three confidence intervals have coverage greater than 0.95 in every case. As expected, both \(CI_b\) and \(CI_{IV}\) are conservative with practically a 100% coverage. This is already apparent in 3 dimensions. Despite the loss of efficiency, however, both yield shorter intervals than \(CI_{LA}\) in 3 dimensions for all three sample sizes. This can most likely be attributed to \(\Gamma\) having multiple roots close to unity, implying that the lag-augmented estimator converges at a rate slower than the IVX and the LS estimators. This advantage more or less vanishes in 4 dimensions and in 5 dimensions \(CI_{LA}\) is the clear winner. Intuitively, the dimension of the confidence regions is quadratic in \(d\) and the loss suffered by projection methods therefore quickly sets in. Interestingly, this phenomenon seems less pronounced for higher sample sizes. Another key result in Table (ref) is that \(CI_b\) is wider than \(CI_{IV}\) in every case. This is counter to the fact that the least squares estimator should be more efficient. One possible explanation is that \(\Gamma\) has roots that converge to 1 at a slower rate than \((1 / n)^{\beta}\) limiting the loss of efficiency of the IVX estimator. Another possible explanation is that finite sample behaviour of \(\hat{t}^2_{IV}\) is different from the asymptotic \(\chi^2\)-distribution for the sample sizes considered here, resulting in a confidence region for \(\Gamma\) with slightly lower coverage but still with conservative coverage when projected onto the first coordinate.

Predictive regression testing

Fix \(\alpha = 0.1\). We compare three methods to test \(H_0:\gamma = 0\) against the alternative \(H_A:\gamma\neq 0\) in the predictive regression model. The first one uses lag-augmentation and is based on the test-statistic in equation (ref). We denote this test by \(\phi_{LA}\). The other two tests employ the Bonferroni strategy described in Section (ref). In particular, for \(\alpha_1 = 0.05\), \(\alpha_2=0.05\), and \(CR(\alpha_1)\) a confidence region for \(\tilde{\Gamma}\) with uniform asymptotic level \(\alpha_1\), we define the test \(\phi_n = \mathbf{1}(\inf_{\tilde{\Gamma}\in CR(\alpha_1)}\hat{\sigma}^{-2}_Y\hat{t}^2_{0|\tilde{\Gamma}} \le q_{d-1, 1-\alpha_2})\). We consider \(CR(\alpha_1) = CR_b(\alpha_1)\) and \(CR(\alpha_1) = CR_{IV}(\alpha_1) \) denoting the corresponding tests by \(\phi_b\) and \(\phi_{IV}\).\footnote{Here, \(CR_b(\alpha_1)\) and \(CR_{IV}(\alpha_1)\) are confidence regions for \(\tilde{\Gamma}\) and not the entire matrix \(\Gamma\).} By Lemma (ref), the tests will have uniform asymptotic level. Computing the two latter test statistics involves an optimization problem. As in the case of the projection confidence intervals, it is much costlier to compute \(\phi_b\) and we employ the EAM-algorithm described in Appendix (ref) as a practically feasible solution.

We perform two sets of simulation experiments for the three tests. To verify that the uniform guarantees hold, we consider a sequence of \(\tilde{\Gamma}\) in the mixed regime. Throughout, we let \(d=4, 5, 6\) and \(\tilde{\Gamma}\in\mathbb{R}^{(d-1)\times (d-1)}\) is chosen as above. We also consider the case \(\tilde{\Gamma}=I\) so that \(\tilde{X}_t\) is a random walk, i.e., \(\tilde{\Gamma}\) is in the non-stationary regime. First we investigate the size properties of the three tests for different sample sizes. The results are depicted in Figure (ref). Evidently, \(\phi_b\) and \(\phi_{IV}\) quickly achieve a rejection rate well below the significance level for all three dimensions and in both regimes. This is in line with the theory since the Bonferroni corrections inherent in these tests result in tests with conservative sizes. For \(\phi_{LA}\) the asymptotics take a little longer to set in. At around \(n=150\) it achieves the correct size for 3 and 4 dimensions, but convergence is slower in 5 dimensions.

To compare the power of the three tests at finite sample sizes we consider a sequence of alternatives increasingly closer to \(0\). In particular, we let \(\gamma = \delta\mathbbm{1}\) for \(\delta\in\{0.005, 0.01,..., 0.1\}\) and fix the sample size at \(n=100\). In both regimes and for all three choices of \(d\), \(\phi_b\) and \(\phi_{IV}\) vastly outperform \(\phi_{LA}\) correctly rejecting the null around 90% of the time in the mixed regime for \(\delta=0.04\) compared to a rejection rate of only around 18% for \(\phi_{LA}\). It looks as though \(\phi_b\) is slightly better than \(\phi_{IV}\) especially as the dimension increases, but the two are close overall. Interestingly, both tests seem to fare better in the mixed regime than in the stationary regime although we would expect the LS and the IVX estimator to converge at a slower rate in the mixed regime. This might be related to the Bonferroni correction and the shape of the confidence regions in the different regimes. The power might be improved upon by using more efficient corrections than Bonferroni, see e.g. elliott2015,jansson2006optimal. The performance of \(\phi_{LA}\) does not depend much on the regime, but it does seem to slightly improve with the size of the dimension. The latter observation also holds for the other two tests and is probably a reflection of the fact that the alternative is easier to detect in larger dimensions.

Conclusion

We proved two major uniform asymptotic results for the sample covariance matrices of VAR(1) processes with potential unit roots. First we showed that \(S_{XX}\) and \(S_{X\epsilon}\) can be uniformly approximated by their Gaussian counterparts \(S_{YY}\) and \(S_{Y\rho}\). This result was used to derive another uniform approximation involving integrals of Ornstein-Uhlenbeck processes. While uniform asymptotic results akin to those presented here have been given in the literature for specific sequences of \(\Gamma\), this is the first time anything has been proven that is truly uniform over the parameter space \(\Theta\).

As an application of the uniform approximation results, we showed how to construct confidence regions for \(\Gamma\) with uniform asymptotic level. Similarly, we proved that the IVX methodology and lag augmentation also lead to uniformly valid inference if done properly.

Funding. The authors gratefully acknowledge financial support from Novo Nordisk Foundation through Grant NNF20OC0062958 and from Independent Research Fund Denmark | Natural Sciences through Grant 9040-00215B.