Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
81,393 characters · 28 sections · 80 citation commands
Cointegration in large VARs
The importance of cointegration in economics stems from the seminal papers granger1981 and engle_granger1987. For example, as they show, monthly rates on 1-month and 20-year treasury bonds are cointegrated, which means that they are both non-stationary, but they have a stationary linear combination. A lot of other variables in macroeconomics and finance such as price level, consumption, output, trade flows, interest rates and so on are non-stationary, and, thus, are potentially subject to cointegration. When dealing with non-stationary time series, it is always a question whether one should work with levels or with differences. For multivariate settings such as vector autoregression (VAR), the choice of model would depend on whether the series is cointegrated or not.
There are several ways to test for the presence of cointegration (see, e.g., maddala for the detailed description of various methods). One popular approach relies on checking whether the residuals from regressing one of the coordinates on the remaining ones are stationary. It is based on engle_granger1987 and was later extended in phillips_ouliaris1990. Another widely used technique is due to S{\o}ren Johansen johansen1988, johansen1991\footnote{A related approach was also proposed in stock_watson1988.}. This approach assumes VAR structure and relies on the likelihood ratio. It tests a null hypothesis of at most $\rho$ cointegrating relationships versus an alternative of between $\rho$ and $r>\rho$ cointegrating relationships. The Johansen test turns out to be related to the eigenvalues of some random matrix and has a non-standard asymptotic distribution.
Neither of the approaches is commonly used in the analysis of large systems. However, in many situations the data turns out to have both cross-sectional and time dimensions being large. Natural examples are given, e.g., by financial data (stock prices, exchange rates, etc.) or by monthly data on trade and investments between countries (where the number of pairs of countries is large). The main difficulty with applying the above cointegration testing approaches to high-dimensional settings is that both of them require $N$, the cross-sectional dimension of a time series, to be fixed and small. For the former approach, the larger is $N$, the more regressions we need to run and interpreting their results becomes ambiguous. For the latter approach, it turns out, the asymptotic theory stops providing a good approximation for moderate values of $T$. The test starts to over-reject the null of fewer cointegrations in favor of the alternative (see, for example, ho_sorensen1996, gonzalo_pitarakis). For instance, the simulations reported in Table 1 of gonzalo_pitarakis indicate that the empirical rejection rate based on the 95% asymptotic critical value for the no-cointegration hypothesis (constructed using the asymptotics of johansen1988, johansen1991) is $20\%$ at $N=5$, $T=30$. It eventually decreases to the desired $5\%$ as $T$ grows, but even at $T=400$ the empirical size is still $6\%$.
After size distortions became clear econometricians developed various procedures for correcting over-rejection. Popular methods include finite sample correction (e.g., Reinsel_Ahn, johansen_correction) and bootstrap (e.g., Swensen, cav_et_all). Such modifications help to restore correct size for moderate values of $T$ (e.g., $T=50$), when $N$ is of a smaller order (e.g., $N=5$). Yet, the question of larger $N$ remained open.
A recent ground-breaking paper onatski_ecta shows that the over-rejection when testing the null of no cointegration can be mathematically explained by considering an alternative asymptotics in which both $N$ and $T$ go to infinity jointly and proportionally. When $N$ is large such asymptotic regime turns out to better suit the finite sample properties of the data. While onatski_ecta point this out, they do not provide alternative testing procedures.
The observation that large VARs might have analyzable joint limits opens a new area of research. That is, one can try to develop various sophisticated joint asymptotics and derive, as a result, appropriate ways to test for cointegration in settings where both $N$ and $T$ are large. This, however, requires new tools rooted in the random matrix theory. In our paper we propose such cointegration tests and introduce asymptotic theorems which make testing possible. Let us describe our main results.
By modifying the Johansen likelihood ratio (LR) test, we come up with a way to examine the presence of cointegration in a time series when the cross-sectional dimension, $N$, and the time dimension, $T$, are of the same order. We consider a vector autoregression (VAR) of order $1$ in the error correction representation
where $\Delta X_t:=X_t-X_{t-1}$, errors $\{\varepsilon_t\}$ are Gaussian with $N\times N$ covariance matrix $\Lambda$, and $\Pi$, $\mu$ are unknown parameters.\footnote{Section (ref) discusses possible generalizations of our approach to VAR($k$) for general $k$.} We do not restrict $\Lambda$ to be diagonal, which means that any cross-sectional correlation structure is allowed. Such cross-sectional heteroscedasticity assumption may be important in applications, where one expects variables, e.g., countries, to be correlated. In contrast, many previous approaches relied on specific forms of covariance $\Lambda$, see, e.g., Breitung_Pesarann_2008, Bai_Ng_2008 for reviews and ZhangPanGao_2018 for recent developments.
Our model allows for an arbitrary linear trend in $X_t$, which is a desirable feature when one goes to applications. Yet, it also encompasses as a special case a model without a trend ($\mu=0$). The way we deal with the data turns out to be invariant to $\mu$, and we do not propose any different procedure even if a researcher has a prior knowledge that there is no trend. (We discuss this in more detail in Section (ref).)
We are interested in whether the $N$-dimensional non-stationary process $X_t$ is cointegrated or not. That is, whether there is a non-zero vector $\beta$ such that $\beta'X_t$ is trend stationary. Our main contributions lie in the construction of the appropriate test, analysis of its asymptotic distribution, and computation of critical values (some of them are reported in Table (ref) of Section (ref)).
Our procedure for cointegration testing relies on the following steps. First, we de-trend $X_t$ by subtracting $\frac{t}{T}\left(X_T-X_0\right)$ from $X_t$. Then we follow Johansen's procedure and regress both first differences, $\Delta X_t$, and lagged de-trended $X_t$ on a constant. That is, we de-mean the data. Then we calculate the squared sample canonical correlations between the residuals from those regressions. This corresponds to the eigenvalues from the modification of the Johansen LR statistics obtained by using the lagged de-trended version of $X_t$ and not de-trended first differences $\Delta X_t$.
The de-trending procedure can be interpreted in the following way: we take the first differences of the observed variables, de-mean them, and re-sum back. This is equivalent (up to a constant $X_0$ which disappears after we further de-mean the sums) to our de-trending. Such interpretation is related to ERS which analyzes unit root testing in the presence of a trend $d_t$. Under the null hypothesis of a unit root and when the trend $d_t$ is linear, ERS suggests to de-mean the first differences.
As discussed in phillips_detrending, de-trending plays an important role in improving the performance of cointegration tests. In our procedure there is an additional reason why de-trending matters. It leads to an unexpected connection with the Jacobi ensemble, a familiar object in random matrix theory, but a novel one in econometrics.\footnote{We remark that de-trending is implicitly present in the proofs of onatski_ecta,onatski_joe, yet in our work it becomes a significant ingredient, rather than a technical detail.} Eventually, this connection is the main technical tool leading to our asymptotic results and construction of the test.
We show that after proper rescaling our test under the null hypothesis of no cointegration converges to the sum of the first $r$ elements of the Airy$_1$ point process (we formally define that process in Section (ref)). Airy$_1$ process is a known object in the random matrix theory and its marginal distributions can be computed in various ways. We also present in Section (ref) a simulation study of the speed of convergence, which supports the limiting results even for moderate values of $N$ and $T$. In Table (ref) of that section we demonstrate the significant improvement over the finite sample behavior of the Johansen test and its corrections reported in gonzalo_pitarakis.
Along with the size we also report power simulations in Section (ref). The power depends on the choice of the alternative and we perform several experiments with random rank $1$ matrix $\Pi$ and random initial condition $X_0$ of varying magnitudes. In many cases the power is very close to $100\%$. In particular, we found high power even for moderate $N,T$, e.g., $N\approx30,\, T\geq150$. In general, the results are encouraging and we see quite large power envelopes.
We also present a small empirical illustration of our testing procedure on weekly S$\&$P100 log-prices over $10$ years. This gives us approximately $500$ observations across time and $T/N\approx5$, corresponding to high power and close to $5\%$ sample rejection rate in simulations. Log-prices are known to have a unit root, and we do not find any strong evidence that they are cointegrated.
Let us briefly indicate technical aspects of our proofs and their mathematical novelty. The key observation that we make is that a small perturbation of the matrix arising in the modified Johansen test has an explicit distribution of its eigenvalues. This distribution is called Jacobi ensemble, and its usual appearances in statistics include sample canonical correlations for two sets of independent data (as opposed to highly dependent in our case, see the discussion in Section (ref)) and multivariate analysis of variance (MANOVA), see, e.g., Muirhead_book. Our method has two main components which are new, as far as we are aware. First, our perturbation of the model is based on a replacement of a certain permutation in the matrix formulation of the modified Johansen test by a uniformly random orthogonal matrix. The second ingredient is a challenging computation of matrix integrals\footnote{Our arguments have some similarities with proofs in Hua.} leading to the identity of the law of the perturbed matrix with the Jacobi ensemble.
Put it otherwise, we discover an exactly-solvable\footnote{A random model is called exactly-solvable or integrable (see borodin_gorin_review for an overview) if there exist explicit formulas for the expectations of non-trivial random variables describing the system. Such formulas provide a basis for asymptotic analysis of the systems of interest. This is in contrast to generic models where explicit formulas are rarely available.} model in a small neighborhood of the (random) matrix of the Johansen cointegration test. Let us center our attention on this exactly-solvable model. It can be used as an initial point for perturbative arguments leading to asymptotic theorems for our and possibly other modifications of the Johansen test. In addition, our exactly-solvable model is not isolated, but rather it is a representative of a whole class of similar cases. We expect that our approach works in several other situation in which various test statistics in the vector autoregression context can be understood through (other) instances of the Jacobi ensemble. Justifications of this point of view are contained in Sections (ref) ,(ref), and (ref).
Section (ref) describes the setting and the main objects of interest. Section (ref) presents asymptotic results, while Section (ref) gives a sketch of their proofs. Section (ref) shows supporting Monte Carlo simulations and Section (ref) applies our test to S&P$100$. Additional results and extensions are presented in Section (ref). Finally, Section (ref) concludes. All proofs, unless otherwise noted, are in Appendix.
We consider an $N$-dimensional vector autoregressive process of order $1$, VAR($1$), based on a sequence of i.i.d. mean zero Gaussian errors $\{\varepsilon_t\}$ with non-degenerate covariance matrix $\Lambda$. That is,
where $\Delta X_t:=X_t-X_{t-1}$ and $\Pi,\,\mu$ are unknown parameters. We do not impose any restrictions on $\Lambda$, thus, we allow for arbitrary correlations across coordinates of $X_t$. The process is initialized at fixed $X_0$. We outline possible extensions to autoregressions of higher order, VAR($k$), in Section (ref).
We are interested in analyzing whether $X_t$ is cointegrated. That is, whether there exists a non-zero $N\times r_0$ matrix $\beta$ such that $\beta'X_t$ is (trend) stationary.
If there exists a $N\times r_0$ matrix $\beta$ of rank $r_0$ such that $\beta'X_t$ is (trend) stationary, but there does not exist a $N\times(r_0+1)$ matrix $\tilde{\beta}$ of rank $r_0+1$ such that $\tilde{\beta}'X_t$ is (trend) stationary, then we say that $X_t$ is cointegrated of order $r_0$. As shown in engle_granger1987, johansen_book the necessary condition for $X_t$ to be cointegrated of order $r_0$ is ${\rm rank}(\Pi)=r_0$ and, thus, there exist two $N\times r_0$ matrices $\alpha$ and $\beta$ of rank $r_0$ such that $\Pi=\alpha\beta'$. For the sufficiency one also needs to require an extra non-degeneracy condition involving $\alpha$, $\beta$, and $\Gamma_i$. This additional requirement is used to rule out the $I(2)$ processes, i.e., such that their second differences are stationary, while the first differences are not.
We test the null hypothesis $H_0$ of no cointegration, which is equivalent to ${\rm rank}(\Pi)=0$ or $\Pi\equiv 0$. The complement to the null is ${\rm rank}(\Pi)>0$. Yet, in order to design out test we set ${\rm rank}(\Pi)\in[1,r]$, where $r$ is a fixed finite number, to be our alternative hypothesis $H_1$. Thus, our alternative, in line with Granger representation theorem, can be interpreted as having at most $r$ cointegrating relationships. While the test has different power depending on $r$ and the true data-generating process, when working with real data for which the true value of ${\rm rank}(\Pi)\equiv r_0$ is unknown one can use any $r$ in order to reject $H_0$.\footnote{In practice we recommend using reasonably small values of $r$, since our asymptotic theorems are valid in the regime of fixed $r$ and large $N,T$.} Formally, we test
Our testing procedure is a modification of the Johansen likelihood ratio (LR) test johansen1988, johansen1991.\footnote{Johansen LR test is based on the maximization of the Gaussian likelihood.} We focus on the small $r$ regime, which differs from the classical Johansen LR test, where $r=N$ is more commonly used as an alternative. In that sense, our approach is close to the maximum eigenvalue ($\lambda_{\max}$) test. That turns out to be important in the large $N,T$ setting, see Section (ref). Our approach consists of several steps:
Step 1. De-trend\footnote{We discuss the role of de-trending in Section (ref).} the data and define
Note that we also do a time shift. This shift is in line with Johansen test, where lags and first differences are regressed on the observables.
Step 2. De-mean the data. That is, regress de-trended lags, $\tilde X_t$, and first differences $\Delta X_t$ on a constant. Define the residuals from those regressions as
Step 3. Calculate the squared sample canonical correlations between $R_0$ and $R_1$. That is, define
where here and below the notation $X^{\ast}$ means transpose of the matrix $X$ (transpose-conjugate whenever complex matrices appear). Then calculate $N$ eigenvalues $\lambda_1\geq\ldots\geq\lambda_N$ of the matrix $S_{10}S^{-1}_{00}S_{01}S^{-1}_{11}$. The eigenvalues solve the equation
Step 4. Form the test statistic
The subscript $N,T$ in (ref) indicates that we modify Johansen LR test to develop the large $N,T$ asymptotics. This statistic after centering and rescaling will be compared with appropriate critical values to decide whether one can reject $H_0$ or not (see Theorem (ref)).
Throughout the proofs and extensions, we will be using an alternative way to write the residuals and matrices $S_{ij}$. For this let us define the demeaning operator $\mathcal P$. It is a linear operator in a $T$--dimensional space defined by its matrix
$\mathcal P$ is an orthogonal projector on the space orthogonal to the vector $(1,1,\dots,1)$. By definition $$ R_{0it}=[\Delta X\mathcal P]_{it}=\Delta X_{it}-\frac{1}{T}\sum_{\tau=1}^T \Delta X_{i\tau},\quad R_{1it}=[\tilde X\mathcal P]_{it}=\tilde X_{it}-\frac{1}{T}\sum_{\tau=1}^T \tilde X_{i\tau}. $$ Using the fact that $\mathcal P^2=\mathcal P$,
Let us emphasize that our test differs from the original Johansen test in the fact that we use $\tilde X_t$ instead of $X_{t-1}$. Note that $\tilde X_t$ can be viewed as a rank $1$ perturbation of $X_{t-1}$. Hence, our test statistic is a finite rank perturbation of the original Johansen procedure.
In this section we formulate our main asymptotic results, explain and discuss the role of the precise setting we chose, and indicate directions for generalisations. We provide a sketch of the proof of our asymptotic theorem in Section (ref).
Our results use the Airy$_1$ point process. Thus, let us introduce it before we formulate the main theorems. The Airy$_1$ point process is a random infinite sequence of reals $$ \mathfrak a_1>\mathfrak a_2>\mathfrak a_3>\dots $$ which can be defined through the following proposition.
The law of the first coordinate $\mathfrak a_1$ is known as the Tracy-Widom $F_1$ distribution; its distribution function can be written in terms of a solution of the Painleve $II$ differential equation.
We remark that from the computational point of view (ref) gives an efficient way to access the distribution of various functions of $\{\mathfrak a_i\}_{i=1}^{\infty}$.\footnote{An even faster way uses tridiagonal matrix models dumitriu_edelman.} From the theoretical perspective, one would like to have a more structural definition, which can be used for the analysis. Such definitions exist, yet, unfortunately, none of them is particularly simple.\footnote{ There are several equivalent ways to define Airy$_1$ point process: through Pfaffian formulas for the correlation functions Forrest_spectr, Tracy_Widom, through combinatorial formulas for the Laplace transform Sodin, through eigenvalues of the Stochastic Airy Operator RRV.}
Figure (ref) shows densities for the random variables $\sum_{i=1}^r \mathfrak a_i$ for $r=1,2,3$. As $r$ grows, the skewness of $\sum_{i=1}^r \mathfrak a_i$ decreases, the expectations of $\sum_{i=1}^r \mathfrak a_i$ tend to $-\infty$ and the variances tend to $+\infty$.
Quantiles of the distribution of $\sum\limits_{i=1}^r \mathfrak a_i$ serve as critical values for testing the hypothesis $H_0$ of no cointegrations against the alternative of at most $r$ cointegrations. We report the quantiles for $r=1,2,3$ in Table (ref).
Note the shifts by $\frac{2}{N}$ in the definition (ref) of $\mathfrak p$ and $\mathfrak q$. Such finite sample correction is inspired by Johnstone_Jacobi. Theorem (ref) remains valid without these shifts, i.e., for the simplified $\mathfrak p=2$, $\mathfrak q=\frac{T}{N}-1$. However, our simulations show that for $T/N<6$ the shifts improve the speed of convergence and, thus, the finite sample behavior of the test, see Sections (ref) and (ref) for details.
We sketch the proof of Theorem (ref) in Section (ref) and give full details in Appendix. The two main ideas that we use are:
It is natural to ask about an extension of Theorem (ref) for the case of $r$ growing together with $N$. We distinguish two subcases here:
In both cases we expect that (after proper adjustment of $c_1$ and $c_2$) the limit in (ref) becomes Gaussian. Although we are not going to pursue this direction here, for the slow growth case the proof of asymptotic normality can probably be obtained by the same methods as we use in Theorem (ref): we can again use the Jacobi ensemble for the computation. For the Jacobi ensemble individual eigenvalues $\lambda_i$ (in the regime of growing $i$ and $N-i$) become asymptotically Gaussian, hence, also their sums.\footnote{ORourke proves asymptotic Gaussianity for eigenvalues of Gaussian random matrices and we expect the same proof to work for the Jacobi ensemble case, see also BPZ for more details in the complex setting.}
For the linear growth case the situation is more complicated. Our present tools only allow us to prove an asymptotic upper-bound of the form
for any $\epsilon>0$ (for $\epsilon=1$ this also follows from the results of onatski_ecta, onatski_joe),\footnote{The value of the constant $c_3$ can be computed through certain integrals involving the equilibrium measure of the Jacobi ensemble, see (ref) and (ref) for more details.} while we expect the expression to be $O(1)$. Yet, there exist very general statements on asymptotic Gaussianity of linear statistics of functions of random matrices (see, e.g., Guionnet_Novak, Mingo_Popa), and, therefore, there is little doubt in the fact that asymptotic distribution is Gaussian. Note, however, that these methods usually give very limited information about asymptotic variance and it is unclear at this moment how to find a reasonably simple explicit formula for it.
An important feature of our setting and our asymptotic results is that $r$ is kept finite as $N,T\to\infty$. (As far as the authors know, Eq. (ref) in Theorem (ref) is the very first statistic for which the computation of the distributional limit became possible in the $N,T\to\infty$, $T/N\in [2+\gamma_1,\gamma_2]$ regime.) In practice, when the dimensions are given to us in the data, this means that $r$ should be much smaller than $N$ and $T$. Notice that for the purpose of rejecting $H_0$ any statistic with known asymptotic distribution can be useful. It could be tempting to instead choose $r=N$ and test $H_0$ of no co-integrations against the alternative of arbitrary number of cointegrations, yet, we believe that in many situations small $r$ allows to infer much more and is a good choice by the following reasons:
The ultimate goal of the cointegration analysis is to find the cointegrating relationships. Turns out, the feasibility of that task crucially depends on their true number. Imagine for a second that the true number of cointegrating relationships is large, say $N/2$, then they span an $N/2$--dimensional subspace of $\mathbb{R}^N$. Thus, the number of unknown parameters which need to be estimated is $N^2$ up to a multiplicative constant. On the other hand, the total number of observations that we have is $NT$. Hence, in our asymptotic regime $N,T\to\infty$, $T/N\in[2+\gamma_1,\gamma_2]$, we have only finitely many observations per estimated parameter and consistent estimation is unlikely. The conclusion from this dimension-based heuristics is that in the situation when the number of cointegration relationships is large, our data does not contain enough information for the consistent estimation of all cointegrating relationships and, in a sense, the problem is hopeless. On the other hand, the same heuristics shows that we are in a much better scenario, if we have some a priory reasons to expect that either the number of cointegrating relationships is small, or if this number is in a small neighbourhood of $N$ (in the latter case instead of the space of cointegrating relationships we can estimate its low-dimensional orthogonal complement). The former situation precisely corresponds to our choice of small $r$ in Theorem (ref); the first step towards the asymptotic analysis in the latter situation is discussed in Section (ref).
In the small $N$ cointegration analysis described in detail in johansen1988, johansen1991, johansen_book, rejection of the hypothesis of no cointegration is the first step. The next step is the estimation of the true rank $r_0$ of $\Pi$ in (ref) by subsequent rejections of the hypotheses of the rank being smaller than $2,3,\dots, r_0$. Finally one wants to consistently estimate $r_0$ cointegrating relationships, closely related to $\Pi$ itself.
Switching to our limit regime $N,T\to\infty$, $T/N\in [2+\gamma_1,\gamma_2]$, there are two separate cases. If $r_0$ also grows linearly in $N$, then the full success of the above program is unlikely: as outlined in Section (ref), we are trying to consistently estimate too many parameters from too few data points.
On the other hand, if true $r_0$ is small, then the situation is different (another potentially feasible case is that of small $N-r_0$, which we do not discuss here). We believe that our present results --- Theorem (ref) and the new approach underlying it --- open a path to the full realization of the cointegration analysis for growing $N$. Let us outline the reasons.
From the technical point of view a key ingredient which needs to be developed is the asymptotics of squared sample canonical correlations $\lambda_1,\dots,\lambda_N$ and corresponding eigenvectors for the model (ref) with $\Pi$ of a fixed and not growing rank $r_0$. In particular Theorem (ref) should be extended from $H_0:\:r_0=0$ to arbitrary finite $r_0$. A natural approach here is by further developing the link to random matrices established in our proof of Theorem (ref) via the Jacobi ensemble. This leads to the theory of spiked random matrices. A basic question of this theory is to infer the information on a large-dimensional small rank matrix $B$ from observing $A+B$, where $A$ is a matrix of pure noise (e.g., one can use $A=\frac{1}{2}\left(Y_N+Y_N^*\right)$ in the notations of Proposition (ref)). This is a well-developed theory, see e.g., johnstone_spiked, BBP for seminal contributions and spiked_china, johnstone_onatski2020 for the most recent progress. In order to obtain a connection with the VAR-framework (ref), we treat the $rank(\Pi)=r_0$ case as a finite rank deformation of the “pure noise case” $\Pi=0$, which we consider in this text. Hence, by combining our present results with the spiked random matrix theory, we expect to establish the complete asymptotic theory for any finite as $N\to\infty$ value of $r_0$. (This will require further technical and conceptual efforts and, therefore, is left for future work.)
Suppose that the true number of cointegrating relationships is small, say, $r_0=1$, yet we are trying to reject the null of no cointegrations using the LR-statistic with $r=N$, i.e.,
in the notations of Theorem (ref). As we explain at the end of Section (ref), the standard deviation of (ref) is $O(1)$, i.e., after subtracting a constant as in Eq. (ref), statistic (ref) stays finite as $N,T\to\infty$. Hence, in contrast to (ref) there should be no $N$--dependent rescaling in a test based on (ref). This should lead to very different powers of tests based on statistics (ref) and (ref).
Indeed, we expect the main difference in the joint law of $(\lambda_1,\dots,\lambda_N)$ under the $H_0$ and under the alternative of one cointegrating relationship to be in the behavior of $\lambda_1$ (which is in line with $\ln(1-\lambda_1)$ being the likelihood ratio test statistic in such setting). Change of $\lambda_1$ by $\nu\in\mathbb{R}$, i.e., $\lambda_1\to \lambda_1+\nu$, is on the same scale as the standard deviation of (ref), and $\nu$ needs to be large in order for the test statistic to detect this change and to reject $H_0$. On the other hand, (ref) is changing in the same situation by a number of order $N^{2/3} \nu$, which is much larger, and, hence, rejection of $H_0$ is more likely. Therefore, one could anticipate that (ref) leads to a test of higher power than (ref) in this situation.\footnote{We do not analyze (ref) in this text and, thus, we do not provide any rigorous justification of the above. We remark that the existing in the literature similar comparisons, e.g., LST2001, paruolo2001, do not apply in our situation, since they only deal with the small $N$ case.}
The fact that we assume $r$ is finite and only use first $r$ largest eigenvalues allows us to use the remaining ones to verify the plausibility of our model specification. That is, we can check whether the statistical properties of $\lambda_{r+1},\dots,\lambda_N$ agree with the predictions of our model: for any finite $r$ we expect the empirical distribution function $\frac{1}{N-r} \sum_{i=r+1}^N \mathbf 1_{\lambda_i\le x}$ to converge to the CDF of the Wachter distribution (under $H_0$ this can be deduced from our Theorems (ref) and (ref), while more general result under $H_1$ with finite $r$ follows from onatski_ecta). Our empirical example in Section (ref), indeed, shows the match of the S$\&$P$100$ data with the Wachter distribution (see Figure (ref)).
The proof of Theorem (ref) rests on the notion of the Jacobi ensemble. Let us first define it and then provide a sketch of the proof of Theorem (ref).
When $N=1$, (ref) is the Beta distribution. For general $N$, the eigenvectors of random Jacobi-distributed $M$ are uniformly distributed, while $N$ eigenvalues $x_1\ge \dots\ge x_N$ admit an explicit density with respect to the Lebesgue measure given by
where $Z(N;p,q)$ is an (explicitly known) normalization constant.
The Jacobi ensemble is widely used in statistics Muirhead_book, theoretical physics Mehta, and random matrix theory forrest. There are numerous tools for studying it\footnote{Just to mention some: Pfaffian point processes Mehta, combinatorics of moments Bai_Silverstein, variational problems for log-gases BenArous, Kadell integrals Kadell, Schwinger-Dyson/Loop equations Johansson1998, tridiagonal models Killip_Nenciu, Macdonald processes BorG.}, and as a result its asymptotic properties (usually as $N\to\infty$) are known in great detail.
Two particularly important appearances of the Jacobi ensemble in statistics are in the multivariate analysis of variance (MANOVA) and in the canonical-correlation analysis, see Johnstone_Jacobi for more detailed discussions and more examples. Let us explain the latter setting, as it has some resemblance with the Johansen test. We start with two rectangular matrices of data: $X$ of size $N\times T$ and $Y$ of size $K\times T$. One can interpret $X_{it}$ (and $Y_{it}$) as the $i$th variable at time $t$.
If $N=K=1$, then one can think about data sets $X$ and $Y$ as $T$ observations of two variables. Then one can measure the (linear) dependence between these variables by computing the (squared) sample correlation coefficient:
A direct computation shows that if the elements of $X$ and $Y$ are i.i.d. mean zero Gaussians, then the law of (ref) is given by Beta distribution for any $T=1,2,\dots$.
A generalization of (ref) to $N,K\ge 1$ defines the squared sample canonical correlations $\{r_i\}$ as solutions to the equation
where $S_{XY}= X Y^*$, $S_{YY}=Y Y^*$, $S_{YX}= Y X^*$, $S_{XX}= X X^*$ (by $X^*$ we mean transpose of $X$ when $X$ is a real matrix and conjugate-transpose of $X$ when $X$ is complex). Generically the equation (ref) has $\min(N,K)$ non-zero solutions. They have a variational meaning. For instance, the maximal solution of Eq. (ref) is the maximal squared correlation coefficient between $a^* X$ and $b^*Y$, where the maximization goes over all $N$--dimensional vectors $a$ and $K$--dimensional vectors $b$.
Now suppose that the columns of $X$ are i.i.d. $N$--dimensional mean $0$ Gaussians with (non-degenerate) covariance matrix $\Lambda$ and the columns of $Y$ are i.i.d. $K$--dimensional mean $0$ Gaussians with (non-degenerate) covariance matrix $\Lambda'$. Assume for simplicity that $N\le K$ and $N+K \le T$. If $X$ and $Y$ are independent, then it can be shown that the (squared) sample canonical correlations $r_1>r_2>\dots>r_N$ have the law of the eigenvalues of the Jacobi ensemble $\mathbf J(N; \frac{K-N+1}{2},\frac{T-N-K+1}{2})$.
At this point we observe both a similarity and a difference between the above instance of the Jacobi ensemble and the matrix appearing in the Johansen test. On one hand, the latter also deals with sample canonical correlations of two data sets. On the other hand, the data sets are no longer independent, instead one is obtained from another by a deterministic linear transformation. So are the roots of Eq. (ref) related to the Jacobi ensemble? There is evidence in both directions. First, computer simulations quickly reveal that in one-dimensional case the distribution of the single eigenvalue in the Johansen test is not the Beta distribution. Yet, second, recent results of onatski_ecta,onatski_joe show that the Law of Large Numbers for the empirical distribution of the eigenvalues appearing in the Johansen test (with roots of Eq. (ref) being a particular case) matches the one for the Jacobi ensemble (with shifted dimension parameters) in the limit as $N,T\to\infty$. Those articles were asking for an explanation.
To prove Theorem (ref), we first need to establish the following central statement: a small perturbation of the Johansen test matrix, obtained by replacing the deterministic summation matrix in its definition by a random analogue, exactly matches the Jacobi ensemble. We further show that the perturbation vanishes in the limit as $N,T\to\infty$, thus, allowing us to obtain the asymptotics of the variants of the Johansen test from the known asymptotic results for the Jacobi ensemble.
The proof of Theorem (ref) is based on the following idea. Looking at Eq. (ref) when $\Pi=0$, one can notice that matrices entering into the test and given by Eq. (ref) can be expressed in terms of the matrix of data $X$ and the lag operator mapping $X_t\to X_{t-1}$. A computation shows that since we deal with the de-trended and de-meaned data, we can replace the lag operator with its cyclic version, which maps $X_1$ to $X_{T}$ rather than $X_0$. The latter is an orthogonal operator whose eigenvalues are roots of unity of order $T$. Then, the idea is to replace this operator by uniformly random orthogonal operator. From that we proceed in two steps:
The full proof of Theorem (ref) is given in Appendix. In addition to the real case, we also simultaneously prove a similar statement for complex matrices, encountering Jacobi ensemble of Hermitian matrices. Generally complex settings are rare guests in economics. Yet, they play a major role in the spectral analysis of time series data and in other areas, such as high energy physics.
By combining Theorem (ref) with known asymptotic results for the Jacobi ensemble we can obtain the asymptotics of test statistic (ref) in various regimes.
In this section we illustrate the performance of our test via Monte Carlo simulations. We consider both size (rejection rate) and power.
First, we compare the finite sample performance of our approach versus Johansen's LR test and one of its corrected versions for (reasonably small) $T=30$ and $N=5,\dots,10$. The finite sample correction takes the form $\frac{T-N}{T}$ and was suggested in Reinsel_Ahn. (Let us note that there are several more advanced empirical finite sample correction procedures for Johansen's LR test and its variations, we refer to onatski_joe for a recent comparison of some of those for a wide range of values of $N$ and $T$.) Table (ref) summarizes the results (numbers closer to $5$ mean better performance). The $LR_{N,T}$ column is our test and the last two columns are from gonzalo_pitarakis. They correspond to the same null of no cointegration, but use a different from ours alternative hypothesis. Our alternative is at most $r$ cointegrating relationships (in Table (ref) we use $r=1$, which means we look at the largest eigenvalue). The LR and RALR tests consider “at most $N$” cointegrations as $H_1$, which means that they use the sum of all eigenvalues. We can see that in finite samples when both $N$ and $T$ are of a similar magnitude, our approach significantly outperforms the alternatives based on small $N$, large $T$ asymptotics.
To illustrate the performance of our test as the sample size increases, we fix the ratio $T/N\equiv c$ and plot the empirical size as a function of $N$. This is shown in Figure (ref), where the target is $5\%$ rejection rate. The picture suggests that the test has rejection rate close to $5\%$. The green solid line corresponding to $c=4$ achieves $5\%$ rejection rate very fast. Other three curves overshoot $5\%$ by couple percents (e.g., both $c=5$ and $c=6$ are always below $7\%$). Moreover, the larger is $c$, the higher is the corresponding rejection curve. E.g., $c=10$ curve (blue, short dash) is strictly above other curves.
To improve the finite-sample behavior for large $c$ we suggest to ignore the $\frac{2}{N}$ correction in $\mathfrak p$ and $\mathfrak q$ in Theorem (ref), Eq. (ref), when $T/N$ is large. That is, to use simplified formulas instead: $\mathfrak p=2,\,\mathfrak q=T/N-1$. We recalculate the empirical rejection rate under the simplified formulas for $\mathfrak p,\mathfrak q$ in Figure (ref). Under the simplified formulas, the larger is $c$, the smaller is over-rejection and the closer the curve is to $5\%$ line. As can be seen by comparing Figures (ref) and (ref), for small $c$ the formula with the correction leads to better rejection rates, while for large $c$ it is the opposite, and we do not gain from finite sample correction. The conclusion is to use the simplified formula when $c=T/N$ is at least $6$.
An important feature of Figures (ref) and (ref) is that as $N,T\to\infty$ with fixed $T/N$ the results improve towards the perfect match to $5\%$. This is in contrast to various finite sample corrections of Johansen's LR test used earlier. In particular, examination of Table $2$ in onatski_joe reveals that each procedure has its own pairs $(N,T)$ where the results are close to perfect (empirical size of the test is very close to the desired $5\%$) and some others where the results are much worse, yet rejection rates in general do not improve as the sample size increases keeping $T/N$ fixed.
When one analyzes the power of a test, the question is which data generating process (dgp) to use under an alternative $H_1$. In our setting the space of alternative dgps has growing with $N$ dimension. Thus, there is no clear choice of the proxy alternative hypothesis to use in simulations. Hence, instead we proceed with various random alternatives.\footnote{Another approach would be to stick to some ad hoc matrix $\Pi$. Yet, we do not expect that to affect our simulations in a qualitatively significant way.} Therefore, the corresponding power is also random. For illustrational convenience, we only report averages of those powers.
To analyze the power of our test we conduct several experiments. In all of them the errors $\varepsilon_{it}$ are generated as i.i.d. $\mathcal{N}(0,1)$. In the first experiment we randomly sample a matrix $\Pi$ of rank $1$. We do this by generating two uniformly random $N$--dimensional unit vectors, $u,v$, such that the non-zero eigenvalue of $uv^{\ast}$ is negative. Then we set $\Pi=uv^{\ast}$. By construction, $\Pi$ has rank $1$, singular value $1$, and one eigenvalue between $-1$ and $0$ (others are zero),\footnote{As $N$ goes to infinity the non-zero eigenvalue is of order $N^{-1/2}$.} so that $X_t$ is non-stationary, while $\Delta X_t$ is stationary. The average power constructed from such random alternative is shown in Figure (ref). As in the previous subsection we fix the ratio $T/N$ and plot the average power as a function of $N$. Following our recommendation, we use simplified formulas for $\mathfrak p$ and $\mathfrak q$ when $T/N\geq6$. We can see that the average power quickly reaches $100\%$ for all ratios of $T/N$. The larger is that ratio, the faster we reach $100\%$. This is due to the fact that higher ratio means higher time span. This gives the process more chances to accumulate the effects of the presence of cointegration.
In the second experiment, we randomly generate a symmetric matrix $\Pi$ of rank $1$. We do this by generating a uniformly random $N$--dimensional unit vector $v$. Then we set $\Pi=-\lambda vv^{\ast}$, where $\lambda$ goes from $0$ to $2$. The coefficient $\lambda$ equals to the non-zero eigenvalue of $-\Pi$. The fact that it lies between $0$ and $2$ guarantees that $X_t$ is non-stationary, while $\Delta X_t$ is stationary. Figure (ref) shows the power as a function of $\lambda$ for $N=100$ and various values of $T$ corresponding to $T/N=\{4,5,6,10\}$ as in the previous experiments. The larger is $-\lambda$, the larger is power, which eventually reaches $100\%$. When $\lambda=0$, the dgp corresponds to the null $H_0$. Thus, all curves start from $\approx5\%$, which corresponds to the empirical size of the test. We can see that again the larger is $T/N$, the faster we reach $100\%$. The reason is the same as in the previous simulation.
Drawing the intuition from the unit root testing literature (e.g., unitroot_X0_ME, unitroot_X0_HLT), we also analyze the sensitivity of power to the choice of initial condition, $X_0$.\footnote{Note that $X_0$ does not affect rejection rates in Section (ref), because under $H_0$ it cancels out.} In two previous experiments we set $X_0=0$. Figure (ref) shows how the power against random alternative of $1$ cointegrating relationship for symmetric $\Pi$ changes for various magnitudes of $X_0$. To be more specific, we redo the same simulations as in the previous paragraph, but start with $X_{i0}\thicksim$ i.i.d. $\text{std}_0\cdot\mathcal{N}(0,1)$ for each Monte Carlo draw. We consider $T/N=5$ and $T/N=10$. Curves with std$_0=0$ are the same as on Figure (ref) (they are also represented with the same style and color on both pictures). We can see that the larger is the magnitude of $X_0$, as measured by its standard deviation std$_0$, the slower the power reaches $100\%$.
In contrast to the previous paragraph, the power against random alternative of $1$ cointegrating relationship when $\Pi$ is asymmetric (as in Figure (ref)) does not exhibit any substantial changes with respect to the magnitude of $X_0$. Hence, we do not redraw the analogue of Figure (ref) for non-zero $X_0$.
Overall, the simulations suggests the good asymptotic performance of our test procedure both under $H_0$ and $H_1$. The theoretical analysis of the power remains an open and challenging question.
In this section we illustrate our testing strategy by analyzing cointegration in log prices of various stocks. The search for cointegrations is a part of a stock market strategy called pairs trading, see, e.g., pairs_trading and references therein. For us this is a convenient testing ground, as both $N$ and $T$ are large.
We use logarithms of weekly S$\&$P100 prices over ten years: $01.01.2010-01.01.2020$, which gives us $522$ observations across time. The time range is chosen so that we do not need to worry about potential structural breaks due to financial crisis of $2007-2008$ and due to COVID-19. S$\&$P100 includes $101$ (because one of its component companies, Google, has 2 classes of stock) leading U.S. stocks with exchange-listed options. We use 92 of those stocks (those which were available for the whole ten years span and only one of two Google's stocks). More details on those stocks are in Section (ref) in Appendix. Therefore, $N=92$, $T=521$ and $T/N\approx 5.66$.
Figure (ref) shows the histogram of eigenvalues which solve Eq. (ref) for the S$\&$P100 data. The key feature of this histogram is that it closely resembles the Wachter distribution, which density is shown by thick orange line in Figure (ref). This distribution governs the asymptotics of the eigenvalues of the Jacobi ensemble (see Section (ref) in Appendix for the details). In particular, we see a precise match between supports of the histogram and of the Wachter distribution. The resemblance is in line with our Theorem (ref). Indeed, if we believe that the true data-generating process (ref) has $\Pi=0$, then the theorem combined with the asymptotics of the Jacobi ensemble of Theorem (ref) predicts convergence to the Wachter distribution. The conclusion is that Figure (ref) shows a close match between theoretical predictions and real data. Simultaneously, this figure is consistent with $H_0$ of no cointegration.
We compute our test statistic for $r=1,2,3$ using the S$\&$P100 data. In neither of the cases the value of the statistic is large enough for a statistically significant rejection of $H_0$. If $r=1$, then the value is $-0.27$, which corresponds to approximately the $0.78$ quantile of the asymptotic distribution shown in Figure (ref). For $r=2,3$ the values are closer to the right tails of the distribution, but still below the (one-sided) $0.95$ quantiles. Hence, we do not see evidence towards the presence of cointegration in S$\&$P100 stock prices for the last 10 years.
Let us describe possible extensions and modifications of our results. In Subsection (ref) we consider non-Gaussian errors $\varepsilon_t$. Subsection (ref) looks at the model without an intercept $\mu$ (linear trend in $X_t$) and considers the effect of de-trending vs. no de-trending in such setting. Subsection (ref) investigates the performance of our test when the true process follows higher order of autoregression. Finally, in Subsection (ref) we discuss testing the hypothesis $\Pi=-I_N$ using the same approach as for $\Pi=0$.
The result of Theorem (ref) is obtained under the assumption that the errors $\varepsilon_t$ in Eq. (ref) are Gaussian. However, we believe that it should be possible to remove this restriction and it is reasonable to expect that the very same statement should hold for any (independent and identical across time $t$) distribution of $\varepsilon_t$, as long as it has sufficiently many moments. The underlying reason for this belief is the so-called universality phenomenon in the random matrix theory: asymptotic local spectral characteristics of a random matrix are almost independent of the distributions of the matrix elements, see, e.g., ErdosYau, Tao_Vu for general reviews and HanPanZhang_2016, HanPanYang_2018 for the recent work in contexts of multivariate analysis of variance and canonical correlations. In particular, for the Wigner matrices ($Y+Y^*$, where $Y$ is a square matrix with real i.i.d. entries), it is known that the asymptotic behavior of the largest eigenvalues depends only on the first two moments (expectation and variance) of the distribution of an individual matrix element. Note, however, that our asymptotic result in Theorem (ref) holds for any choice of the first two moment of Gaussian noise $\varepsilon_t$: in Eq. (ref) the covariance matrix $\Lambda$ is arbitrary and any shift in expectation can be absorbed into the parameter $\mu$. Hence, we conjecture that the result of Theorem (ref) would hold for any distribution of $\varepsilon_t$, as long as it is sufficiently well-behaved.
In order to test this conjecture we made simulations for the case when elements of $\varepsilon_t$ are non-gaussian, but i.i.d. across both $i$ and $t$ (corresponding to a diagonal covariance matrix $\Lambda$). We ran Monte-Carlo simulations for three different distributions for $\varepsilon_{it}$: uniform on the interval $[0,1]$, uniform on $3$ points $\{1,2,3\}$, and the product of two independent $\mathcal{N}(0,1)$ random variables. In each case for $T=900$, $N=300$, and small values of $r$, we do not see any significant changes in the distribution of $\sum_{i=1}^r \ln(1-\lambda_i)$ from the limit in Theorem (ref). However, things go differently when the distribution has heavy tails. In the forth experiment we chose $\varepsilon_{it}$ to be Cauchy-distributed, and then the distribution of $\sum_{i=1}^r \ln(1-\lambda_i)$ dramatically changed from what we saw in the Gaussian case. Hence, we conclude, that the existence of at least some number of moments of the errors $\{\varepsilon_t\}$ should be necessary for the validity of the conjecture. Rigorous proof of the conjecture remains an important problem for the future research.
One of the important steps in our testing procedure is de-trending the data. Moreover, the exact form of the de-trending that we use (we subtract the slope of the line which connects $X_0$ at $t=0$ and $X_T$ at $t=T$) is an ingredient substantially used in the proof of Theorem (ref). Yet, we expect that this is just a technical artifact and statements similar to Theorem (ref) should also be true with other forms of de-trending or in models where this step is not needed at all. Let us provide some evidence.
Our model allows for any linear trend and works even if the true value of $\mu$ is zero. Yet, when one has such prior knowledge it may seem natural to ignore the de-trending and de-meaning steps (Steps $1$ and $2$ in Section (ref)) and proceed without them. In this section we compare tests with and without de-trending for the model without $\mu$:
If both de-trending and de-meaning steps are omitted, we get (the small $r$ version of) the classical Johansen test statistic for the model (ref). If only de-trending is omitted, we get the Johansen test statistic corresponding to our original model (ref). When both de-trending and de-meaning are implemented, we get our procedure described in Section (ref). To compare the asymptotic properties of those procedures, we perform a Monte Carlo simulation. Results are reported in Figure (ref).
As Figure (ref) suggests, the densities have almost identical shape. Figure (ref) plots all three of them together as well as their versions with subtracted mean. As we can see, after subtracting the mean, all three densities are identical. This suggests robustness of the Airy$_{1}$ point process in our asymptotic results of Theorem (ref) (which correspond to Figure (ref)) and predicts that some modifications, which preserve the limiting distribution, are possible.
Let us emphasize that our limit theorems currently only apply to Figure (ref), but not to the settings of Figures (ref) and (ref). The mismatch of the means in Figure (ref) makes one suspect that some modifications of the constants in the asymptotic theorems are needed as soon as we start slightly adjusting the setting.
In this subsection we discuss the performance of our test when the true data generating process (dgp) is VAR($k$), $k>1$. That is
and the no cointegration situation corresponds to $\Pi\equiv 0$.
Even if $k>1$, we expect to see the Tracy-Widom distribution and marginals of Airy$_1$ process (as in Theorem (ref)) in the asymptotic behavior of the squared sample canonical correlations from Section (ref) under mild restrictions on $\Gamma_i$. The belief is based on the universality intuition of the random matrix theory. However, we do not expect the scalings (such as the coefficients $c_1$ and $c_2$ in Theorem (ref)) to remain the same. The most plain analogy is the dependence of centering and scaling in the classical Central Limit Theorem on the underlying process. Closer to our context is the asymptotic behavior of sample covariance matrices: when the data is i.i.d., the empirical distribution of the eigenvalues of the sample covariance matrix converges to the Marchenko-Pastur law (and the largest eigenvalues concentrate near the right edge of this distribution), while data with general covariance structure leads to much richer limits, see Bai_Silverstein and references therein. Figuring out (even heuristically) any formulas for the scaling coefficients $c_1$ and $c_2$ as functions of $\Gamma_i$ is a challenging open problem for the future research.
There is an important family of cases where one can hope that the formulas for $c_1$ and $c_2$ from Theorem (ref) remain valid (perhaps, with minor modifications). This is when $\Gamma_i$ are small and evolution of $X_t$ given by (ref) can be treated as a small perturbation of the VAR($1$) process. One way to formalize the “smallness' of $\Gamma_i$ is by requiring them to be of small rank (cf.\ discussion of rank in Section (ref)) and of small norm.
In order to investigate the above conjecture, we run Monte Carlo simulations for $k=2$ and $\Gamma_1$ of rank $1$. We consider the null $\Pi\equiv0$ and calculate our test statistic based on $r=1$ under VAR($2$) and VAR($1$) data generating processes and compare their asymptotic distribution. Under VAR($1$) $\Gamma_1\equiv0$, while for VAR($2$) we use $\Gamma_1=0.5E_{11}$ and $\Gamma_1=0.5E_{12}$, where $E_{ij}$ is a matrix which has $1$ on the intersection of row $i$ and column $j$ and $0$ everywhere else; the components of the noise $\varepsilon_{it}$ are i.i.d. $\mathcal{N}(0,1)$. The asymptotic distribution of our test statistic with $r=1$ for various dgps is shown in Figure (ref). We see that the distributions on each panel of Figure (ref) are close to each other. Hence, testing based on Theorem (ref) in such a VAR($2$) setting remains valid. Yet, if we consider more general situation of $\Gamma_1=\theta E_{11}$ or $\Gamma_{1}=\theta E_{12}$, then we observe in simulations (not shown) that the quality of approximations significantly deteriorates as $\theta$ grows to $1$.
The main result of Theorem (ref) is a development of a test for the hypothesis $\Pi=0$ in the VAR($1$) model (ref). One could also try to understand for which other $\Pi$ can the testing be possible. The asymptotic distribution depends on the choice of $\Pi$, and it is not possible to estimate $\Pi$ consistently in our regime of $T$ and $N$ growing to infinity proportionally. Thus, for general $\Pi$ the problem seems infeasible at this point. However, there is another particular choice of $\Pi$ for which an approach very similar to Theorem (ref) still works: $\Pi=-I_N$. Denoting this hypothesis $H_0^{w.n.}$, where $w.n.$ stays for the white noise, the data generating process (Eq. (ref)) becomes:
In other words, we are now testing the hypothesis that the time series $X_t$ is independent across time $t$ against various VAR($1$) alternatives. Here is one setup where such testing can be relevant. Suppose that we want to forecast some variable $Y$ and we chose some model for it. After estimating parameters of the model we obtain residuals $X_t$. If we know that $X_t$ are independent, then they are unforecastable and we cannot further improve our forecasting model. To check the above we can take the residuals $X_t$ and then apply our white noise hypothesis testing procedure.
Let us introduce an adaptation of the Johansen statistic to $H_0^{w.n.}$. As in Eq. (ref), we use the notation $\mathcal P$ for the de-meaning operator projecting on the hyperplane orthogonal to $(1,1,\dots,1)$.
Following johansen1988,johansen1991 we are going to use the de-meaned data $X_t \mathcal P$. For the increments $\Delta X_t$ in addition to the conventional de-meaning we use an extra modification: we deal with cyclic increments $\Delta^c X_t$ defined as: $$ \Delta^c X_t=
$$ While it might seem bizarre to subtract the last observation from the first one, if we recall that our current hypothesis of interest \eqref{eq_wn_hypothesis} ignores the time ordering, then this becomes less controversial. Note also a shift of index by $1$, as compared to the conventional $\Delta X_t$, which is compensated by the lack of shift $t\to t-1$ in $X_t$, as compared to Eq.~\eqref{eq_detrending}. Our choice of definition of $\Delta^c$ is important for the following precise asymptotic results. As in Section (ref), the conventional Johansen statistic should be thought of as a finite rank perturbation of the modified version that we now introduce.
We further define $N$ numbers $\lambda_1\ge \lambda_2\ge \dots\ge \lambda_N$ as $N$ roots to the equation
Equivalently, $\{\lambda_i\}$ are eigenvalues of $S_{10}^{w.n.} (S_{00}^{w.n.})^{-1} S_{01}^{w.n.}(S_{11}^{w.n.})^{-1}$.
The proof of Theorem (ref) follows a similar strategy as Theorem (ref) and we refer to Section (ref) in Appendix for details; in particular, the proofs rely on yet another novel appearance of the Jacobi ensemble.
The remaining straightforward step to obtain the asymptotics of various statistics built on the eigenvalues $\{\lambda_i\}$ is to combine Theorem (ref) with asymptotic results for the Jacobi ensemble presented in Section (ref). This is in the spirit of Theorem (ref).
Note that the hypothesis $H_0^{w.n.}$ implies the maximal amount of cointegrating relationships: each of the $N$ components of $X_t$ is already stationary. A reasonable alternative hypothesis $H_1$ is the presence of $N-r$ cointegrating relationships. For simplicity, let us concentrate on the case $r=1$. Then the alternative can be also interpreted as a presence of a single growing factor. In this situation we expect the smallest eigenvalue $\lambda_N$ to be a good test statistic. We see in numerical simulations that $\lambda_N$ is bounded away from $0$ under $H_0^{w.n.}$. It can also be formally proved by combining Theorem (ref) with asymptotics of the Jacobi ensemble from Section (ref). Thus, if $\lambda_N$ is close to $0$, then we are able to reject $H_0^{w.n.}$. The same simulations indicate that $\lambda_N$ starts to be close to $0$ when there are at most $N-1$ cointegrating relationships. Hence, the test based on $\lambda_N$ should have a good asymptotic power. We leave rigorous justifications of this observation till future research, and for now only mention the following heuristics: the stationary linear combinations of $X_t$ are strongly correlated with the same linear combinations of $\Delta^c X_t$; on the other hand, the growing linear combinations of $X_t$ have very weak correlation with the same linear combinations of $\Delta^c X_t$ (cf. correlations of a one dimensional random walk with its increments). Hence, if the latter are present, the smallest canonical correlations of $X \mathcal P$ and $\Delta^c X\mathcal P$, which coincide with the eigenvalues of the matrix $S_{10}^{w.n.} (S_{00}^{w.n.})^{-1} S_{01}^{w.n.} (S_{11}^{w.n.})^{-1}$, should become small leading to close to $0$ value of $\lambda_N$.
The paper presents a cointegration test which has desirable empirical size when $N$ and $T$ are of the same magnitude. To our knowledge, this is a first paper which constructs and analyzes asymptotic properties of a test that does not suffer from significant distortions (such as over-rejection) for comparable $N$ and $T$. The test is based on the Johansen LR test and incorporates some additional steps. First, our procedure reinforces the importance of de-trending in cointegration testing. It turns out, that de-trending is crucial for deriving desirable asymptotic properties. (E.g., only after de-trending one can rewrite the lagged process as a linear function of its first differences.) Second, our asymptotic results reveal and explain an unexpected connection between the Johansen cointegration test and the Jacobi ensemble --- a classical ensemble of the random matrix theory whose previous appearances in statistics include multivariave analysis of variance (MANOVA) and sample canonical correlations for independent sets of data.
On the theoretical side the next step would be to go from null hypothesis of zero cointegration to analyzing the behaviour of our test under $r$ cointegrations. This will allow us to calculate the power of the test, reinforcing our simulational findings in Section (ref), as well as to perform tests of $r$ versus $r+1$ cointegrations.
On the empirical side it would be interesting to apply our test to other data sets beyond what is presented in Section (ref). Annual cross-country data provides a natural example of our setting where the number of years and countries is comparable. Another example arises if one considers network-type settings which evolve over time (e.g., as in bykh). Data on trade or on foreign direct investment can potentially be non-stationary, especially if we focus on largest and the most active countries. Moreover, although such monthly data is available, for many countries it only covers $~20$ years. Thus, we have $T\approx200$. If we look at directed pairs across $10$ largest countries, this gives us $N=90$ cross-section units, which fits ideally in our setting. Classical cointegration tests are known to perform poorly in the above settings. However, the asymptotic results of our paper open up a possibility of detecting the presence of cointegration in such time-series data.