Asymptotics of Cointegration Tests for High-Dimensional VAR($k$)

Asymptotics of cointegration tests for high-dimensional VAR($k$)

\address[Anna Bykhovskaya]{Duke University} \email{[email removed]}

\address[Vadim Gorin]{University of California at Berkeley} \email{[email removed]}

\@startsection{section}{1} \z@{1.0\linespacing\@plus\linespacing}{0.5\linespacing} {\normalfont}{Introduction}

Starting with the pioneering work of sims, vector autoregressions (VARs) became a workhorse model in macroeconomics and other fields. Many key time series in macroeconomics and finance (e.g., consumption and output) are nonstationary, and the properties of VARs can be very different depending on whether one is dealing with a stationary or nonstationary series. Moreover, there is a further subdivision to be accounted for in the case of nonstationary series: it is important to understand whether the data are cointegrated---that is, whether there exists a stationary nontrivial linear combination within the considered series (e.g., the log of consumption minus the log of output is stationary while the series themselves have unit roots).

Classical tools for testing cointegration (see, e.g., johansen_book, maddala, and juselius) fail to achieve the desired finite sample performance when the number of time series, $N$, is large. Thus, they are not commonly used in such settings, and the design of proper tools to handle cointegration under a large $N$ remained an open problem for years (see, e.g., choi15). Recently onatski_ecta, onatski_joe and BG have opened a new avenue based on the “$T/N$ converging to a constant” asymptotic regime. However, the testing procedures of these texts cover only VAR($1$), while, in practice, researchers rarely confine themselves to VARs of order $1$, instead usually considering at least two lags. Indeed, as noted already in pagan, “most applications of Sims’ methodology have put the number of lags between four and ten.” Since pagan the lengths of available time series and computing power have only increased, thus allowing researchers to work with even more complex models. Hence, it is important to generalize and extend the above papers to a VAR($k$) setting, which is the main topic of our text.

Our paper analyses a family of tests for the absence of cointegration for nonstationary VAR($k$), such as the Johansen likelihood ratio (LR) test (johansen1988, johansen1991) and related Hotelling--Lawley and Pillai--Bartlett tests (see, e.g., gonzalo_pitarakis1995 and references therein) as $N$ and $T$ jointly and proportionally go to infinity. The shared feature of these tests is that their statistics are based on the squared sample canonical correlations between certain transformations of current changes and past levels of the data. The main contribution of our paper is in the asymptotic analysis of these canonical correlations. First, we show that for VAR($k$) with general $k$, under the null of no cointegration (and some additional technical conditions) the empirical distribution of the squared sample canonical correlations converges to the Wachter distribution. As a corollary, we deduce the first-order deterministic limits of the above test statistics. Second, we introduce a modification of the testing procedure and prove much more refined results in the modified setting. By computing the exact asymptotic behavior of the probability distributions of individual canonical correlations after proper recentering and rescaling, we are able to compute the critical values for the test of no cointegration with correct asymptotic size as $N$ and $T$ jointly and proportionally go to infinity.

We remark that there is a wide scope of literature devoted to the corrections of Johansen's LR test and its relatives (originally developed based on fixed $N$, large $T$ asymptotics) for large values of $N$ (see, e.g., Reinsel_Ahn, johansen_correction, Swensen, cav_et_all, and onatski_joe). The distinguishing feature of our work is that we are not trying to correct the finite $N$ asymptotic statements, which stop working for large $N$, by introducing various empirical adjustments. Instead, we develop a theoretical framework for working with the large $N$ case directly. One advantage is that our approach explains the general phenomenology and predicts the asymptotic behavior. As a result, the empirical or simulational adjustments for particular values of the parameters of the model are no longer needed.

To achieve the above, in our proofs we use the VAR($1$) results of BG as a cornerstone. The main technical work is devoted to producing recursive arguments, which reduce the VAR($k)$ behavior to that of VAR($k-1$) and eventually to VAR($1$). The central role is played by highly nontrivial projections from the group of orthogonal $T\times T$ matrices to the smaller subgroup of orthogonal $(T-N)\times (T-N)$ matrices. While such projections have previously been used in asymptotic representation theory, to the authors' knowledge, this is their first appearance in the econometrics or statistics context. Thus, many new properties of those projections need to be developed in our framework.

The rest of the paper is organized as follows. Section (ref) describes our setting and provides the first asymptotic results. Section (ref) constructs our modified test and computes its asymptotics, while Section (ref) presents supporting Monte Carlo simulations. Section (ref) illustrates our theoretical findings on S&P100 data and on the prices of cryptocurrencies. Finally, Section (ref) concludes. All proofs are in Sections (ref)--(ref). The accompanying R package is available at the Github \url{https://github.com/eszter-kiss/Largevars}.

\@startsection{section}{1} \z@{1.0\linespacing\@plus\linespacing}{0.5\linespacing} {\normalfont}{First-order asymptotics of sample canonical correlations}

We consider an $N$-dimensional vector autoregressive process of order $k$, VAR($k$), based on a sequence of i.i.d. mean zero Gaussian\footnote{We expect that all our results continue to hold for non-normally distributed errors as long as they have enough moments (cf.\ such distribution-independence results in other high-dimensional models, as in ErdosYau, Tao_Vu, HanPanYang, and FanYang).} errors $\{\varepsilon_t\}$ with nondegenerate covariance matrix $\Lambda$. That is, written in the error correction form,

where $\Delta X_t:=X_t-X_{t-1}$, $D_t$ is a $d_D$-dimensional vector of deterministic terms, such as a constant, a trend or seasonality (extra explanatory variables are also allowed as long as they are observed), and $\Gamma_1,\ldots,\Gamma_{k-1},\,\Pi,\,\Phi$ are unknown parameters. The process is initialized at fixed $X_{1-k},\ldots,X_0$. We do not impose any restrictions on $\Lambda$; thus, we allow for arbitrary correlations across coordinates of $X_t$. In contrast, many previous approaches rely on specific properties of the covariance matrix $\Lambda$; see, e.g., Breitung_Pesarann_2008, Bai_Ng_2008 and ZhangPanGao_2018.

We are interested in the behavior of the squared sample canonical correlations between transformed past levels (lags) and changes (first differences) of the data $X_t$. As shown in johansen1988, johansen1991 (see also Anderson), the correlations are related to whether the process is cointegrated. To be more specific, they appear in the likelihood ratio test for the presence and rank of the cointegration. Let us formally define these correlations. Here and below $^*$ denotes matrix transposition.

Johansen's LR statistic for testing the hypothesis $\mathrm{rank}(\Pi)\le r_1$ (at most $r_1$ cointegrating relationships) versus the alternative $\mathrm{rank}(\Pi)\in (r_1,r_2]$ (between $r_1$ and $r_2$ cointegrating relationships) with $r_2>r_1$ has the form\footnote{We omit a usual scaling factor of $T$ for the statistics (ref), (ref), and (ref).}

the Pillai--Bartlett statistic is

and Hotelling--Lawley statistic is

See gonzalo_pitarakis1995 for a discussion and many references about these statistics.

In Theorem (ref) we show that the empirical measure of eigenvalues of $\mathcal C$ converges (weakly in probability) to the Wachter distribution. The theorem generalizes the results of onatski_ecta from the VAR($1$) to the VAR($k$) setting.\footnote{While onatski_ecta allows the data to be VAR($k$) under restriction (ref), they construct the matrix $\mathcal C$ involved in the statistical testing procedures as if the data were VAR($1$). That is, the formal procedure is based on a misspecified VAR($1$) setting, while we use the true VAR($k$) procedure. As illustrated in Section (ref), using an underspecified VAR can lead to severe size distortions.}

Imposing assumption (ref) can be viewed as a dimension reduction (cf. sparsity assumption). Approximating data with low-rank matrices is a widely used and powerful technique in data science, in machine learning applications such as recommender systems (e.g., movie preference recognition), and in computational mathematics. We refer the reader to lowrk_highdim for theoretical explanations of the suitability of low-rank models and many references to situations in which they are very efficient. In our particular context, the number of unknown parameters in the VAR model (ref) is proportional to $N^2$, and we have access to $N T$ observations. Since $N^2$ and $NT$ are of the same order in the asymptotic regime (ref), the model (ref) can overfit the data. We view the rank restriction (ref) as a natural way to avoid overfitting\footnote{An alternative way to introduce a low-rank assumption into the VAR model is, instead of using the error correction form in (ref), to rewrite the evolution as $X_t=\sum\limits_{i=1}^{k} A_i X_{t-i}+\Phi D_t+\varepsilon_t$. Then, in the spirit of factor models, one can impose the low-rank assumption on $A_i,\,i=1,\ldots,k$. Notice that $A_i=\Gamma_i-\Gamma_{i-1}$ for $i=2,\ldots,k-1$, so that low-rank assumptions on higher-order lags in this and our setting are related. This alternative low-rank restriction complements ours via the rank of $\Pi$: in our setting, the number of cointegrating relationships grows sublinearly in $N$, while in the alternative setting, this number is close to $N$. While none of our theorems directly cover the factor setting, numeric simulations in Section (ref) indicate that the tests that we develop remain useful.} (cf. the discussion in JASA_low_rank_VAR and Wang_Tsay_low_rank_VAR). Section (ref) illustrates that the results obtained under this assumption are consistent with the behavior of large-dimensional financial datasets. Another setting in which we can expect (ref) to be satisfied is when there are a few special coordinates in $X_t$, e.g., some macroeconomic indicators, that mostly drive the behavior of the entire vector $X_t$. This would correspond to the case where the columns of $\Gamma_i$ corresponding to those indicators are nonzero while the other columns are zero.

Figure (ref) illustrates Theorem (ref) for independent standard normal errors and $k=2,\,N=150,\,T=1500$. The parameters are $\Phi D_t=1_N$, $\Gamma_1=0.95E_{12}$, $\Pi=-0.1E_{\cdot1}$, where $1_N$ is an $N\times 1$-column matrix of ones, $E_{12}$ is an $N\times N$ matrix with one at the intersection of the 1st row and the 2nd column and zeros everywhere else, and $E_{\cdot1}$ is an $N\times N$ matrix with ones in the first column and zeros everywhere else. Thus, all the matrices have rank one. The parameters of the Wachter distribution are $\mathfrak p=2,\,\mathfrak q=T/N-k=8$. The single separated (rightmost) eigenvalue corresponds to $\mathrm{rank}(\Pi)=1$, i.e., one cointegrating relationship. Generally, we expect that if there are $r$ separated eigenvalues and the value of $k$ in Procedure (ref) is correctly specified, then there are at least $r$ cointegrating relationships.

The proof of Theorem (ref) is based on treating the setting of (ref) as a small-rank perturbation of a more restrictive setting analyzed in Section (ref). The small-rank assumption in (ref) is crucial for the validity of the theorem, and we expect that the asymptotic behavior changes in situations when (ref) fails. This expectation is supported by the Monte Carlo experiment in Section (ref). In contrast, the assumption that $T>(k+1)N$ can potentially be relaxed. When $T<(k+1)N$, the matrix $\mathcal C$ has deterministic eigenvalues equal to $1$, which should be taken into account in $N\to\infty$ asymptotics. This case can be also addressed by our methods, but we do not continue in this direction.

An important corollary of Theorem (ref) is that it provides the asymptotic behavior of various tests constructed from eigenvalues of $\mathcal C$.

Note that, when $\rho_1=0$, for the Johansen LR and Hotelling--Lawley (HW) statistics, we obtain inequalities rather than equalities. This is because of the singularity of $\ln(1-\lambda)$ and $\frac{\lambda}{1-\lambda}$ at $\lambda=1$: while Eq. (ref) controls average behavior, it does not control individual eigenvalues. Hence, the largest eigenvalue $\lambda_1$ can be arbitrarily close to $1$, so that the LR and HW statistics reach large negative and positive values, respectively. However, for similar statistics in the modified setting of the next section the inequalities turn into equalities (as can be proven by combining Theorem (ref) and Proposition (ref)).

For commonly used tests, one often takes $r_2=N$ (i.e., $H_1:\mathrm{rank}(\Pi)\leq N$), in which case $F^{-1}(\rho_2)=0$. We also remark that the integrals in Corollary (ref) can be explicitly computed in many situations. For instance, the one appearing in the asymptotic of the Pillai--Bartlett statistic for $\rho_1=0$, $\rho_2=1$ is $$ \int_{0}^{1} x\, \mu_{2, \tau-k}(x)\, \mathrm d x=\frac{2}{\tau+2-k}. $$

There are several applications of Theorem (ref) and Corollary (ref):

To draw further economical and statistical conclusions and to develop precise statistical tests and their critical values, one needs to go beyond the first-order asymptotic results of Theorem (ref) and Corollary (ref). In the next section we introduce relevant modifications and develop appropriate second-order asymptotics.

\@startsection{section}{1} \z@{1.0\linespacing\@plus\linespacing}{0.5\linespacing} {\normalfont}{Cointegration test: Second-order asymptotics}

In the regime of $N$ and $T$ growing simultaneously and proportionally, the first-order asymptotics of tests based on the squared sample canonical correlations are given in Corollary (ref). To perform testing and be able to reject at a given significance level, we need to be more precise and find a centered limit, which would be a random variable rather than a constant. To do this, we need to impose additional conditions on $\Gamma_i$, $D_t$, and $\Phi$ in Eq. (ref). Let us first describe the modified procedure and then state the asymptotic results.

Test

We restrict our attention to the case $D_t=1$, i.e.,

The null hypothesis of no cointegration is $H_0:\:{\rm rank}(\Pi)=0$ or $\Pi\equiv0$. The complement to $H_0$ is ${\rm rank}(\Pi)>0$. However, to design our test we use an alternative hypothesis: $$ H(r):\quad {\rm rank}(\Pi)\in[1,r]. $$ As in BG, our test is based on a modification of the Johansen LR test. The Johansen LR test for the original $H_0$ (i.e., $\Pi\equiv0$) versus $H(r)$ is

where $\lambda_1,\lambda_2,\ldots$ are defined in Procedure (ref). Let us describe how our modified test proceeds.

An alternative way to write residuals $\tilde{R}_{i},\,i=0,k$ is via an orthogonal projector: Let $\mathcal W$ be a linear subspace of dimension $N(k-1)+1$ in $T$-dimensional vector space, spanned by vector $(1,1,\dots,1)$ and all rows of matrices $(\Delta X) (L_c^{i})^*$, $1\le i\le (k-1)$, where $L_c$ is a cyclic version of the conventional lag operator and $L_c^{i}$ is its $i$th power, that is, the cyclic lag applied $i$ times. The cyclic lag operator $L_c$ maps a vector $(x_1,x_2,\dots,x_T)$ to $(x_{T},x_1,x_2,\dots,x_{T-1})$. Let $P_{\bot \mathcal W}$ denote the projector on orthogonal complement to $\mathcal W$. Then,

Second-order asymptotics

In this section we show that, under additional restrictions, the eigenvalues $\tilde{\lambda}_i,\,{i=1,\ldots,N}$ are very close (up to $N^{-1+\epsilon}$ for arbitrary $\epsilon>0$) to a known random matrix distribution. From this result we deduce our main theorem (Theorem (ref)), which gives the large $N,T$ limit of the test statistic $LR_{N,T}(r)$ in Eq. (ref). Before we formally state the results, let us define the relevant random matrix distributions.

Definitions

The Jacobi ensemble is a generalization of the Beta distribution to the space of square matrices (when $N=1$, we obtain the Beta distribution). It plays a prominent role in statistics; e.g., it appears in canonical correlation analysis for independent data sets and in multivariate analysis of variance (see, e.g., Muirhead_book).

The marginals of the Airy$_1$ point process can be calculated via various methods (see, e.g., forrest for more details).

Theorems

The null $H_0$ for (ref) is not a point hypothesis, as it does not specify $\Gamma_i,\,i=1,\ldots,k-1$. A simplifying procedure when we are faced with such a composite space of the maintained hypothesis is to assume some fixed values of the parameters as a proxy for the null hypothesis. Along these lines, for the next theorems we are going to introduce additional restrictions and specify the values of $\Gamma_i,\, i=1,\ldots,k-1$. Thus, our model is going to be fully specified (up to a constant $\mu$, which will disappear in the testing procedure). We proceed to implement this approach in testing the hypothesis of no cointegration and introduce the restricted $\widehat H_0$\footnote{We discuss the consequences of using $\widehat H_0$ for testing the null $H_0$ after Theorem (ref).}:

In other words, under $\widehat H_0$ the data generating process turns into

where $\mu$ is an (unknown) $N$-dimensional vector.

The proof of Theorem (ref) relies on two steps. First, we modify our matrix $\tilde{\mathcal C}$ a bit, which leads to a surprising appearance of the Jacobi ensemble, as shown in Section (ref). Second, in Section (ref) we show that the distance between the original model and the modified one becomes small as $N\to\infty$.

Combining Theorem (ref) with known asymptotic results for the Jacobi ensemble, which we recall in Proposition (ref) in Section (ref), we derive the asymptotics of (ref) in the following theorem.

Theorem (ref) gives us the basis of cointegration testing in the large $N,T$ setting. Treating $\widehat H_0$ as a proxy for $H_0$, we can use our asymptotic results to test high-dimensional VARs for the presence of cointegration. Formally, to perform testing, one first needs to calculate the statistic $LR_{N,T}(r)$ following Procedure (ref). We recommend using small\footnote{In Theorem (ref) $r$ is kept fixed as $N$ and $T$ grow. The role of this choice and the motivations for sticking to it are discussed in detail in BG.} values of $r$, such as $r=1,2$, or $3$. Then, one needs to calculate $\frac{LR_{N,T}(r)- r \cdot c_1(N,T)}{ N^{-2/3} c_2(N,T)}$, as in Theorem (ref), and compare the result with quantiles of the sum of Airy$_1$, $\sum_{i=1}^r \mathfrak a_i$. If the rescaled statistic is larger than the $\alpha$ quantile, we reject the null of no cointegration at the $(1-\alpha)$ level. We report the quantiles for $r=1,2,3$ in Table (ref). See also vignette_largevars for more detailed tables for $r=1,\ldots,10$.

Note that, although the asymptotic result (ref) is shown under the restrictions $\widehat{H}_0$, we believe that it extends well beyond $\widehat H_0$: the same asymptotic results and testing procedures continue to hold in many situations with nonzero $\Gamma_i$ in Eq. (ref). While we do not have a full rigorous proof, we expect the following to be true:

For the data generating process (ref), assume that the ranks of all $\Gamma_i$ are bounded, as are the norms of all the matrices and vectors involved in the specification of the process (see Section (ref) for more details). Then conclusion (ref) of Theorem (ref) should continue to hold.

We collect extensive evidence supporting this statement. In Section (ref) we report results from Monte Carlo simulations consistent with it. Further, in Section (ref) we present a precise mathematical conjecture in this direction and give a heuristic argument for its validity. The intuition is that generic small-rank matrices are negligible relative to the scale of the rest of the process and, thus, their addition does not change the asymptotics. For this intuition to hold, it is important to correctly specify the parameter $k$ in the procedure to be equal to (or greater than) its true value. Otherwise (i.e., if we do not regress on the relevant $\Delta X_{t-i}$ in the procedure), the presence of $\Gamma_i$ can have an effect similar to that of the presence of nonzero $\Pi$: it leads to the appearance of special highly correlated linear combinations of rows of $\tilde R_0$ and $\tilde R_k$, which changes the behavior of the largest canonical correlations $\tilde \lambda_i$; see also the simulations in Section (ref).

Theorem (ref) means that under $\widehat H_0$ the largest eigenvalues $\tilde \lambda_i$ are close to $\lambda_+$, which is the right point of the support of the Wachter distribution in Eq. (ref). The relevance of this theorem for cointegration testing stems from the fact that we expect some of the eigenvalues to be much larger than $\lambda_+$ when cointegrating relationships are present. As an illustration, see Figure (ref), where $\Pi$ of rank $1$ leads to the largest eigenvalue being to the right of $\lambda_+$ and separated from the other eigenvalues. The separation is due to the small rank of $\Pi$. However, even if the rank of $\Pi$ is large, we expect the largest eigenvalue to be significantly larger than $\lambda_+$, and, thus, the test remains relevant (see Section (ref)). Providing rigorous results on the consistency of the test is an important task for future research. We present the first result in this direction in Corollary (ref) in Section (ref), where we produce a lower bound on the power of the test against a particular “one cointegrating relationship” alternative and show that the power tends to $1$ as $T/N$ tends to infinity.

\@startsection{section}{1} \z@{1.0\linespacing\@plus\linespacing}{0.5\linespacing} {\normalfont}{Monte Carlo simulations}

Size

We refer to BG for the finite sample size performance of our test for $k=1$. The results for VAR($k$) are similar, and we do not show them in much detail here. For illustration purposes and to represent the comparative statics, Table (ref) reports the empirical size for $T=522,\,N=92$ (those numbers correspond to our empirical example in Section (ref)) for tests based on VAR($k$), $k=1,2,3,4$ procedures.\footnote{Depending on the assumed order of autoregression, we have different numbers of regressors in Procedure (ref).} We can see that the numbers are close to the desired $5\%$ and, for the same $N$ and $T$, a lower order of VAR leads to slightly better results.

$\mathbf{H_0}$ vs. $\mathbf{\widehat H_0}$

An important aspect of our analysis for $k>1$ is the introduction of the additional restrictions $\widehat H_0$ maintained under the null. We would like to check whether Theorem (ref) can hold under the less restrictive $H_0$ instead of $\widehat H_0$. Some theoretical results in this direction are provided in Section (ref). Here we complement them with Monte Carlo simulations. For $N=100$, $T=500$ we simulate the data based on i.i.d. $\mathcal{N}(0,1)$ errors $\varepsilon_{it}$, zero $\Pi$, and nonzero $\Gamma_i$ (i.e., this corresponds to $H_0$ but not $\widehat H_0$). We then compare the density of the test based on the largest eigenvalue $\tilde{\lambda}_1$ ($r=1$ case of Theorem (ref) with statistic $\frac{\ln(1-\tilde{\lambda}_1)-c_1(N,T)}{N^{-2/3} c_2(N,T)}$), with the density of the first coordinate of the Airy$_1$ point process, $\mathfrak a_1$. If the densities coincide, then it means that we can still use the asymptotics from Theorem (ref) to test the null of no cointegration.

Let $E_{ij}$ be a matrix with $1$ at the cell $(i,j)$ and $0$s everywhere else and let $E_{\cdot j}$ be a matrix with $1$s filling the entire column $j$ and $0$s everywhere else. In the first two experiments we take $k=2$. We set $\Gamma_1=0.95E_{11}$ in the first one, which guarantees stationarity of $\Delta X_t$ but allows for strong time correlations in the first coordinate via the $0.95$ factor. In this case the rank of $\Gamma_1$ is $1$. In the second experiment we consider an asymmetric matrix $\Gamma_1=0.1E_{\cdot1}+0.95E_{12}$, which has a close to $1$ singular value because of the $0.95$ factor; the rank of $\Gamma_1$ is $2$ in this case. The results are illustrated in Figure (ref). In the third experiment, we take $k=3$, $\Gamma_1=E_{11}$, and $\Gamma_2=-\tfrac2{9}E_{11}$, so that both matrices are of rank $1$. The value $-\tfrac{2}{9}$ guarantees stationarity of $\Delta X_t$, since $1-z+\tfrac{2}{9}z^2=\left(1-\tfrac1{3}z\right)\left(1-\tfrac{2}{3}z\right)$. The result is shown in Figure (ref). We interpret the outcomes of these three experiments as a strong argument toward the validity of an analogue of Theorem (ref) well beyond the $\widehat H_0$ setting.\footnote{The minor mismatches between densities as in Figures (ref) and (ref) should be expected even under $\widehat H_0$. Theorem (ref) (after multiplication of the result by $N^{2/3}$, as in Eq.\ (ref)) predicts errors of at least ${\rm const}\cdot N^{-1/3}$ in the approximations under $\widehat H_0$.}

Order of VAR

It is essential for the experiments in the last paragraph that the data generating process is VAR(2) and that the procedure we use also corresponds to VAR(2), i.e., $k=2$ in the notations of Sections (ref) and (ref). As illustrated in Figure (ref), using a larger $k$ would lead to similar results, while incorrectly using a VAR(1) procedure when the data generating process is VAR(2) would imply wrong centering and scaling. Moreover, one can spot in Figure (ref) that underestimation of the order of the VAR can be misinterpreted as a presence of cointegration\footnote{Related simulations for $\Gamma_1=\theta E_{11}$ are also reported in BG: For small $\theta$, such as $\theta=0.5$, the VAR(1) procedure still performs well. However, as $\theta$ grows to $1$, the performance quickly deteriorates ($\theta=0.95$ in Figure (ref)).} (largest eigenvalue separated from the rest leading to the large value of the test statistic). However, as we increase the order, the largest eigenvalue becomes inseparable from the rest, and no sign of false cointegration remains present. Thus, practitioners are encouraged to experiment with the order of the VAR to make sure that they are detecting cointegration and not simply using the wrong model.

Note that one should be careful if using classical information criteria for estimating the order $k$ of a VAR in our situation. They are known to be unreliable in high-dimensional settings and may underestimate $k$ (see, e.g., the simulations in gonzalo_pitarakis2002). A possible approach to choosing $k$ is to look sequentially at histograms of eigenvalues at $k=1,2,\dots$. If the outlier eigenvalues larger than $\lambda_+$ exist for the procedures with all $k$ and perhaps move closer to $\lambda_+$ as $k$ grows (corresponding to a decrease in the power of the test), then this is a strong indication of the presence of cointegration. On the other hand, if there is a sharp transition---i.e., outlier eigenvalues are present when we use the VAR($k'$) procedure for $k'<k$ and abruptly disappear at $k=k'$---then this is an indication that the true model is VAR($k$) without cointegration (see Figure (ref)).

All of the above reinforces the importance of using the VAR($k$) rather than the VAR($1$) procedure.

Small ranks

To illustrate the importance of small ranks (e.g., (ref) in Theorem (ref)), we also redo the same procedure for a matrix $\Gamma_1$ of full rank and set $\Gamma_1$ to be $0.95I_N$, where $I_N$ is an $N\times N$ identity matrix. The result is shown in Figure (ref). While the shape and the scale (corresponding to $N^{2/3}$ rescaling in Theorem (ref)) of the distribution remain similar, the location changes. Thus, the small-rank restriction of Eq.\ (ref) is important not only in the context of Theorem (ref) but also for correct centering in possible generalizations of Theorem (ref).

Power

Finally, we simulate the process based on $\Pi\neq0$ to assess the power of our cointegration testing procedure. We refer the reader to BG for many simulations in the $k=1$ case and do not repeat similar experiments here.

For the first experiment, we use $k=2$, $\Gamma_1=0$, and $\Pi=-0.95E_{11}$. Figure (ref) shows the results of the simulation. Two curves are separated; moreover, the black straight line (test distribution) is flatter than the blue dashed curve ($\mathfrak a_1$). First, the separation of the curves is in line with the usefulness of Theorem (ref) in cointegration hypothesis testing, since the test statistic was designed to distinguish between $\Pi=0$ and $\Pi\neq 0$. Second, the distinct variances are due to the fact that (as we expect from a comparison with results on spiked random matrices in the literature; see, e.g., BBP) under the alternative ($\Pi\neq0$) the test needs to be scaled differently: Instead of $N^{2/3}$ rescaling one should use $N^{1/2}$. This result is in line with the power analysis of our test in Section (ref).

For the second experiment, we use $k=2$, $N=150$, $T=1500$, $\Gamma_1=0.95E_{12}$, and $\Pi=-0.8I_\rho X_{t-2}$, where $I_\rho$ is a matrix with ones on the first $\rho$ diagonal elements and zeros elsewhere. The histograms of eigenvalues $\tilde{\lambda}_1,\ldots,\tilde{\lambda}_N$ for $\rho=5,50, 150$ are shown in Figure (ref). As $\rho$ becomes large, we are no longer in the framework of Theorem (ref), and the result about the convergence to the Wachter distribution does not apply. Nevertheless, we observe that the largest eigenvalues are significantly larger than $\lambda_+$ of Theorem (ref). Recall that our testing procedure is based on comparing the largest eigenvalues\footnote{More precisely, logarithms of 1 minus eigenvalues vs. $\log(1-\lambda_+)$.} with $\lambda_+$. Thus, this leads to the conclusion that our test is useful for all values of $\rho\in [1,N]$: the test rejects $H_0$ of no cointegration (i.e., the $\Pi=0$ hypothesis) at a very high statistical significance level.

\@startsection{section}{1} \z@{1.0\linespacing\@plus\linespacing}{0.5\linespacing} {\normalfont}{Empirical illustrations}

S$\&$P$\mathbf{100}$

We illustrate our asymptotic theorems on the S$\&$P100 data. We use logarithms of weekly prices of assets in the S$\&$P100 over ten years (January 1, 2010, to January 1, 2020), which gives us $522$ observations across time. More detailed description of the variables can be found in BG.

For the S$\&$P100 data set we use Procedure (ref) to calculate $\tilde{\lambda}_1,\ldots,\tilde{\lambda}_N$. We do this for various choices of $k$. The case $k=1$ (VAR(1)) corresponds to BG. The results are shown in Figure (ref).

We see a striking match between the histograms and Wachter densities for all $k=1,2,3,4$, which is an indication that the setting of Theorem (ref) is a proper modeling for the S$\&$P data. We do not see any outliers in the largest eigenvalues, which would appear if there were cointegration. Indeed, our test statistics based on Theorem (ref) are $-0.28,\,-0.71,\,-1.07,\,-3.84$ for $k=1,2,3,4$, respectively, while the $5\%$ and $10\%$ critical values are $0.97$ and $0.44$. Because the former numbers are smaller than the latter numbers, we do not reject the “no cointegration” hypothesis.

Cryptocurrencies

In this subsection we redo the calculations for cryptocurrencies instead of S$\&$P stocks. We use the data from crypto_coint_paper (25 series from the Github repository). Logarithms of daily prices for two years (from October $5$, $2017$, to October $4$, $2019$) are shown in Figure (ref). The results of Procedure (ref) used to calculate $\tilde{\lambda}_1,\ldots,\tilde{\lambda}_N$ are shown in Figure (ref).

Similarly to the S$\&$P example in the previous subsection, we see a match between the eigenvalues and the Wachter distribution.\footnote{The Wachter distribution depends on the order of the VAR, $k$, and on the ratio $T/N$. Thus, the orange curves in Figures (ref) and (ref) have different shapes and supports.} However, there is a major difference between Figures (ref) and (ref): The latter has around 3 eigenvalues to the right of the support of the orange curve (Wachter distribution). This is an indication of the presence of approximately $3$ cointegrating relationships. This is reinforced by our test, which has p-values below $0.01$ for all four choices of the order of VAR($k$) ($k=1,2,3,4$).

The difference in results for traditional stocks and cryptocurrencies can be explained by the fact that the cryptocurrency market is still very inefficient and, thus, has numerous trading possibilities. The presence of cointegration can be one such inefficiency.

\@startsection{section}{1} \z@{1.0\linespacing\@plus\linespacing}{0.5\linespacing} {\normalfont}{Conclusion}

High-dimensional data are becoming increasingly widespread in economics and other sciences. Thus, appropriate machinery for handling such data is needed. We believe that the use of random matrix theory is inevitable for the development of the area: As soon as dimensions are high, random matrices start to contribute. Along these lines, in our paper the central role is played by random matrix objects: the Wachter distribution, Airy$_1$ point process, and Jacobi ensemble.

The present paper focused on nonstationary high-dimensional VARs and presented the asymptotic limit of the Johansen LR test for cointegration and its modifications. Because the limit is nonrandom, the appropriate second-order statistic was derived, and a new test for the presence of cointegration was proposed. The new test builds upon the Johansen LR, while having some extra modifications. This new test is suitable for a vector autoregression of order $k$ with an intercept.

The main focus of the present paper is the null of no cointegration. The next essential step is to be able to test whether the cointegration rank is $r$ for $r>0$, i.e., to find the true rank of cointegration. Heuristics for finding the correct value of $r$ can already be seen in our simulations and data sets (cf.\ Figures (ref), (ref), and (ref)): When there are no cointegrations, all eigenvalues (squared canonical correlations) are to the left of the end-point of the support of the Wachter distribution. In contrast, we expect each cointegrating relationship to lead to an eigenvalue between the right end-point of the support of the Wachter distribution and $1$. Identifying the exact conditions under which this heuristic is correct represents an important problem for future research.

\@startsection{section}{1} \z@{1.0\linespacing\@plus\linespacing}{0.5\linespacing} {\normalfont}{Appendix 1: Proofs}

This appendix contains the proofs of Theorems (ref), (ref), and (ref) from the main text.

First, in Section (ref) we collect known statements about the asymptotics of the Jacobi ensemble of Definition (ref), which will be used in our subsequent proofs.

Second, in Theorem (ref) of Section (ref) we introduce a novel random matrix model for the Jacobi ensemble. Our proof of Theorem (ref) proceeds through certain intricate inductive computations of large-dimensional matrix integrals.

Third, in Section (ref) we connect the matrix model of Section (ref) to the cointegration setting: for that we use the rotational symmetry of the Gaussian law to express the squared sample canonical correlations solving (ref) under the hypothesis $\widehat H_0$ in terms of a certain deterministic orthogonal matrix. Replacing this deterministic matrix by a uniformly random one, we arrive at the Jacobi ensemble of Theorem (ref). We proceed by bounding the error in this replacement, which relies on the rigidity estimate for orthogonal matrices (ref), but needs special care due to various matrix inversions involved in our procedures. Eventually, we arrive at Theorem (ref). This theorem is our main technical result. Combining Theorem (ref) with Proposition (ref) from Section (ref), we finish the proof of Theorem (ref) from the main text.

Finally, in Section (ref) we prove Theorem (ref) by combining Theorem (ref) with Proposition (ref) of Section (ref) and general statements about small rank perturbations.

Asymptotic of Jacobi ensemble

In this section we review the asymptotic results for the Jacobi ensemble $\mathbf J(N; p,q)$ introduced in the Definition (ref) as $N\to\infty$.

We assume that as $N\to\infty$, also $p,q\to\infty$, in such a way that

where $\mathfrak p$ and $\mathfrak q$ are two parameters, which stay bounded away from $1$ and from $\infty$ as $N\to\infty$.\footnote{Johnstone_Jacobi suggests to use $\frac{p-1}{N}$ and $\frac{q-1}{N}$ instead of $\frac{p}{N}$ and $\frac{q}{N}$, respectively, in order to improve the speed of convergence. However, we found in BG that for the tests in VAR($1$) case the usefulness of this correction depends on the exact value of the ratio $T/N$ and we are not going to pursue this direction here.} We further define the equilibrium measure $\mu_{\mathfrak p,\mathfrak q}$ of the Jacobi ensemble through:

where the support $[\lambda_-,\lambda_+]$ of the measure is defined via

One can check that $0<\lambda_-<\lambda_+<1$ for every $\mathfrak p,\mathfrak q>1$. Further, define

and note that $$ \mu_{\mathfrak p,\mathfrak q}(x-\lambda_{\pm})\approx \frac{c_\pm}{\pi} \sqrt{|x-\lambda_{\pm}|}, \text{ as } x\to\lambda_{\pm}\quad \text{ inside }\quad [\lambda_-,\lambda_+], $$ where the normalization $\frac{1}{\pi} \sqrt{|x-\lambda_{\pm}|}$ was chosen to match the behavior of the Wigner semicircle law $\frac{1}{2\pi} \sqrt{4-x^2}$ near edges $\pm 2$.

A new model for the Jacobi ensemble

The Jacobi ensemble appearing in Theorem (ref) originates in the following computation of exact distribution. In addition to real symmetric matrices ($\beta=1$ in the usual random matrix notations) it also covers the case of complex Hermitian matrices ($\beta=2$).

The $k=1$ case of Theorem (ref) is established in BG. The proof of Theorem (ref) uses the following three auxiliary ingredients.

A perturbation of the Jacobi ensemble.

In this section we use Theorem (ref) to prove Theorem (ref).

Recall the cyclic shift\footnote{Note that in BG we expressed all the operators in terms of $F=L_c^{-1}$ rather than $L_c$.} operator $L_c$ acting in $T$--dimensional space. Let $V$ be the $(T-1)$--dimensional space orthogonal to the vector $(1,1,\dots,1)$, i.e., $V=\{(x_1,\dots,x_T)\mid x_1+\dots+x_T=0\}$. Note that $V$ is an invariant space for $L_c$ and let $L_V$ denote the restriction of $L_c$ on the subspace $V$.

Take a uniformly-random orthogonal (or unitary if $\beta=2$) operator $\tilde O$ acting in $(T-1)$--dimensional space $V$ and define an operator $\tilde L$ acting in $V$: $$ \tilde L=- \tilde O L_V \tilde O^*. $$

Comparing Proposition (ref) with Theorem (ref) and Remark (ref) one notices that the differences are in restricting on the subspace $V$ (hence, decreasing the dimension by $1$) and in replacement $O\leftrightarrow \tilde L$.

The next proposition explains the effect of the replacement $O\leftrightarrow \tilde L$ on the eigenvalues of the product of projectors in Theorem (ref) and Proposition (ref). We need to introduce some additional notations.

Choose positive integers $k$, $N$, and $\mathcal T$, such that $\mathcal T\ge (k+1) N$ and an arbitrary $N$--dimensional subspace $\mathcal V$ in $\mathcal T$--dimensional space. Let $$f^{k,N,\mathcal T;\mathcal V}: SO(\mathcal T)\to \{0\le x_1\le x_2\le\dots \le x_N\le 1 \}$$ be a map from the group $SO(\mathcal T)$ of orthogonal $\mathcal T\times \mathcal T$ matrices of determinant $1$ to $N$--tuples of reals on $[0,1]$ interval, defined by the following procedure: Take $O\in SO(\mathcal T)$. Let $P$ be the orthogonal projector on the space orthogonal to $O \mathcal V$, $O^2 \mathcal V$,\dots, $O^{k-1} \mathcal V$. Let $P_1$ be the projector on the subspace $P \mathcal V$ and $P_2$ be the projector on the subspace $P O^{k}(I_{\mathcal T}+O)^{-1} \mathcal V$. Then $f^{k,N,\mathcal T;\mathcal V}$ maps $O$ to $N$ largest eigenvalues of $P_1P_2P_1$.

We also need three norms:

Proposition (ref) claims a continuity of map $f^{k,N,\mathcal T;\mathcal V}$. The proof needs care because of the inversions in the definition of the map $f^{k,N,\mathcal T;\mathcal V}$.

The proof of Proposition (ref) relies on three lemmas which we prove later in this section. For these lemmas we write matrices $O_1$ and $O_2$ of Proposition (ref) in the block forms according to the splitting $\mathcal T=(\mathcal T-N)+N$:

Small rank perturbations

In this section we prove Theorem (ref) by combining Theorem (ref) with Proposition (ref) and general statements about small rank perturbations. The key step of the proof is the following observation:

\@startsection{section}{1} \z@{1.0\linespacing\@plus\linespacing}{0.5\linespacing} {\normalfont}{Appendix 2. Discussion of asymptotics under $H_0$ and $H_1$}

The goal of this section is to discuss the asymptotics of the test statistic $LR_{N,T}(r)$ of (ref) under various data generating processes (ref) generalizing $\widehat H_0$ of (ref) and Theorem (ref).

Beyond $\widehat H_0$

We start by working under a slightly more restrictive assumption than (ref) of Theorem (ref). Let $\|A\|_2$ be the spectral norm of a matrix $A$.

We do not expect the conditions in Conjecture (ref) to be optimal. For instance, the Gaussianity assumption can likely be relaxed, as the simulations of BG indicate, and it is plausible that $\mathrm{rank}(\Gamma_i)<C$ condition can be replaced with slow growth of $\mathrm{rank}(\Gamma_i)$, as in Theorem (ref). Nevertheless, we wanted to record Conjecture (ref) in the present form, as a precise statement to be addressed in the future work. We are not giving a proof of Conjecture (ref) here: the required mathematical apparatus does not exist so far. Instead, we are going to provide a heuristic argument for its validity based on our recent results in BG_CCA in a related, yet different setting.

BG_CCA studied the following general setting: let $\mathbf U$ and $\mathbf V$ be two random linear subspaces in $S$--dimensional space with $\dim(\mathbf U)=K$, $\dim(\mathbf V)=M$ and all three numbers $K,M,S$ assumed to be growing to infinity. In addition, suppose that there are $\mathbbm q$ special vectors $\mathbf u_1,\dots,\mathbf u_{\mathbbm q}$ inside $\mathbf U$ and other $\mathbbm q$ special vectors $\mathbf v_1,\dots,\mathbf v_{\mathbbm q}$ inside $\mathbf V$, where $\mathbbm q$ is assumed to stay finite as other parameters grow. We directly observe $\mathbf U$ and $\mathbf V$, but not $\mathbf u_1,\dots,\mathbf u_{\mathbbm q}$ or $\mathbf v_1,\dots,\mathbf v_{\mathbbm q}$. Can we reconstruct $\mathbf u_1,\dots,\mathbf u_{\mathbbm q}$, $\mathbf v_1,\dots,\mathbf v_{\mathbbm q}$, or at least identify their presence by looking at the squared sample canonical correlations between $\mathbf U$ and $\mathbf V$ and corresponding vectors?

The connection to our cointegration tests comes from taking as $\mathbf U$ the space spanned by the $N$ rows of $\tilde{R}_0$, as defined after (ref), and as $\mathbf V$ the space spanned by the rows of $\tilde R_k$. The value ${\mathbbm q}$ corresponds to the cointegration rank and $\mathbf v_1,\dots,\mathbf v_{\mathbbm q}$ correspond to the cointegrating relationships.

While any finite $\mathbbm q$ can be analyzed in a similar fashion, let us stick to $\mathbbm q=1$ case for simplicity, so that we have a single vector $\mathbf u\in\mathbf U$ and another vector $\mathbf v\in\mathbf V$. The most important quantity is the sample squared correlation coefficient $r^2$ between vectors $\mathbf u$ and $\mathbf v$. It turns out that if $r^2$ is large (i.e., close to $1$, because $0\le r^2 \le 1$), then the largest canonical correlation between $\mathbf U$ and $\mathbf V$ is clearly separated from the rest (reminiscent of Figure (ref)) and the corresponding eigenvectors can be used to extract information on $\mathbf u$ and $\mathbf v$. On the other hand, if $r^2$ is small, then the histogram of the canonical correlations does not have such a spiked eigenvalue and all the information about $\mathbf u$ and $\mathbf v$ is washed out. BG_CCA proved the existence of $r^2_{\rm{critical}}\in(0,1)$ separating the above two regimes for a variety of settings for the data generating process for $\mathbf U$, $\mathbf V$, $\mathbf u$, and $\mathbf v$, see also bao2019canonical,yang2022limiting. However, the results of BG_CCA do not address the setting relevant to cointegration and further new ideas would be necessary to find the value of $r^2_{\rm{critical}}$ for cointegration or rigorously prove its existence. Nevertheless, because the cointegration testing is also based on canonical correlations, one expects that the same phenomenology is true for it and, therefore, there should be the following dichotomy:

We now present heuristics in favor of Conjecture (ref) based on this dichotomy. Some of the technical details are omitted as we try to express the key ideas instead.

We now come back to Lemma (ref).

Power

We proceed to our next computation, supplementing Conjecture (ref). This time we would like to explain what changes in the asymptotics, if the data generating process satisfies the alternative $H_1$, rather than the null-hypothesis $H_0$. For the clarity of the exposition, we only concentrate on one particular instance of $k=1$ case here and deal with the $N$-dimensional data generating process

$\varepsilon_t \thicksim\text{i.i.d.}~\mathcal{N}(0,\Lambda)$, $E_{11}$ is the matrix with $1$ in top-left corner and $0$s everywhere else. $\theta$ is a real parameter, we set $\beta=1+\theta$, which implies that the first coordinate of $X_t$ is a scalar process $y_t$ solving

Note that $\xi_t$ are i.i.d.\ $\mathcal N(0,\sigma^2)$ random variables for some constant $\sigma^2$. Because in the notations of (ref) the matrix $\Pi$ now has rank $1$, one hopes that our test statistic of Section (ref) under (ref) behaves significantly differently than under $H_0$ of Conjecture (ref). This would imply that the no-cointegration test based on Theorem (ref) has high power against the rank one alternative (ref). Let us prove that this is indeed true for large values of the ratio $T/N$.

As a corollary, we deduce that our cointegration test has a significant power against rank one stationary alternative, and this power tends to $1$ as $T/N\to\infty$. Here is a precise statement:

Our approach to the proof of Corollary (ref) gives a lower bound on $p(\alpha,\tau)$. Finding exact formulas for the power under $H_1$ of this corollary and under other alternatives remains an important open problem for the future research.