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.
65,913 characters · 15 sections · 56 citation commands
Largevars: An R Package for Testing Large VARs for the Presence of Cointegration
\sloppy
Vector Autoregressions (VARs) are a fundamental tool in econometrics and time series analysis, providing a framework for modeling the dynamic interrelationships among multiple time series. However, as the number of variables in a VAR increases, the complexity of the model grows significantly, posing challenges for both estimation and inference. One critical aspect of analyzing VAR models is testing for the presence of cointegration, which can inform whether a set of non-stationary series share a long-run equilibrium relationship. That is, whether a set of non-stationary time series has a stationary linear combination, see, e.g., johansen_book.
There are several ways to test for the presence of cointegration (see, e.g., maddala for the detailed description of various methods). Yet traditional tests for the presence of cointegration (e.g., likelihood ratio of johansen1988,johansen1991) are not suitable for analyzing large systems, as they tend to significantly over-reject the null hypothesis, see, for example, ho_sorensen1996, gonzalo_pitarakis. To address this issue, bykhovskaya_gorin_1, bykhovskaya_gorin_k propose an approach based on alternative asymptotics, where both the number of coordinates $N$ and the length $T$ of time series are large, and tailored for high-dimensional time series. {\fontseries{m}\fontseries{b}\selectfont Largevars} package implements this approach.
The test implemented in the {\fontseries{m}\fontseries{b}\selectfont Largevars} package is based on the squared sample canonical correlations between transformed past levels (lags) and changes (first differences) of the data, as outlined in Section (ref). Its asymptotic distribution (derived under $N,T\to\infty$ jointly and proportionally) is given by the partial sums of the Airy$_1$ point process, a random matrix object defined in Section (ref). The quantiles of the sums are necessary in order to implement the test; we have tabulated them and included both in the package and in the tables presented in this article. These tables are of independent interest and, with the exception of Table (ref), have not appeared in the literature before. Table (ref) corresponds to the quantiles of the Tracy-Widom distribution, which were also tabulated in Bejan.
Another \proglang{R} package that provides quantiles of random matrix origin is {\fontseries{m}\fontseries{b}\selectfont RMTstat} rmstatpackage, which offers density, distribution, and quantile functions for the Tracy–Widom distribution with parameters $\beta=1,2,4$, based on precomputed tables. In comparison, our package provides Tracy-Widom quantiles for $\beta=1$ in Table (ref), along with nine additional tables containing quantiles necessary for cointegration testing based on the Airy$_1$ partial sums. The method used in {\fontseries{m}\fontseries{b}\selectfont RMTstat} to generate quantile tables is not applicable in our more general setting. Our alternative nontrivial algorithm, along with its \proglang{MATLAB} implementation, is described in detail in Section (ref).
Various \proglang{R} packages have been developed for cointegration testing. For instance, the {\fontseries{m}\fontseries{b}\selectfont urca} package urcapackage includes functions such as \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{ca.jo} (implementing the procedure of johansen1988, johansen1991), \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{ca.po} (for the test of phillips1990asymptotic), and \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{cajolst} (implementing the procedure of lutkepohl2004testing). The {\fontseries{m}\fontseries{b}\selectfont ARDL} package ardlpackage provides functions \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{bounds_f_test} and \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{bounds_t_test}, which implement the Wald bounds-test and t-bounds test for no cointegration by pesaran2001bounds. Additionally, the {\fontseries{m}\fontseries{b}\selectfont bootCT} package bootctpackage offers \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{boot_ardl} function for the bootstrap version of ARDL cointegration tests, as in bertelli2022bootstrap. Several other implementations of cointegration tests exist in other software languages. Despite these offerings, the majority of existing packages are not tailored for high-dimensional settings, and their performance on time-series with large $N$ remains uncertain (see onatski_ecta, onatski_joe for detailed theoretical discussions on over-rejection in classical tests for large $N$). Some of the packages explicitly prohibit the use of large $N$, e.g.\ \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{ca.jo} in the {\fontseries{m}\fontseries{b}\selectfont urca} package does not output test results for $N>10$. Therefore, our new package complements existing software by providing specialized tools with theoretical assurances for high-dimensional settings.
This section explains the theoretical foundations of the cointegration testing and challenges in their practical implementation.
We consider a $N \times (T+1)$ data set represented by columns $X_t$, $0 \leq t \leq T$. These columns are interpreted as observations of an $N$-dimensional vector at $T+1$ time points or as $N$ scalar time series. $X_t$ is said to be cointegrated if there exists a linear combination with coefficients $\beta$ of these time series such that the scalar time series $\beta^\top X_t$, for $0 \leq t \leq T$, is stationary over time, potentially after detrending. Conversely, if every linear combination is non-stationary, we conclude no cointegration exists. In stationary scenarios, there is no growth over time, and correlations exhibit short-range behavior. Non-stationary scenarios, however, show growth and wider correlations, see brockwell1991time,johansen_book for rigorous treatments. In particular, when $N=2$, cointegration implies a long-term equilibrium where the first and second coordinates of $X_t$ move together.
Cointegration has been extensively studied in econometrics, beginning with seminal works by granger1981,engle_granger1987. Many variables in macroeconomics and finance, such as price levels, consumption, output, trade flows, and interest rates, may exhibit cointegration. A classic example is the relationship between interest rates for 3-month and 1-year US Treasury bills, which are cointegrated (see, e.g., BG_review for an illustration). In portfolio management cointegrated stocks give rise to the strategy known as “pairs trading”. This article and its accompanying software discuss statistical procedures for determining whether a given data set $X_t$ demonstrates cointegration. The main advantage of our approach is its applicability to settings with large $N$, in contrast to many existing packages and articles.
For the mathematical setup we assume that the data set is a realization of an $N$-dimensional vector autoregressive process of order $k$, denoted VAR($k$). The process is driven by a sequence of i.i.d. mean-zero errors ${\varepsilon_t}$ with non-degenerate covariance matrix $\Lambda$. In its error correction form, the model reads:
where $\Delta X_t := X_t - X_{t-1}$. The parameters $\mu \in \mathbb{R}^N$, and $\Gamma_1, \ldots, \Gamma_{k-1}, \Pi \in \mathbb{R}^{N \times N}$ are unknown. The process is initialized with fixed values $X_{1-k}, \ldots, X_0$.
(Some authors use a slightly different form: $\Delta X_t=\mu+\Pi X_{t-1}+\sum\limits_{i=1}^{k-1}\tilde{\Gamma}_i\Delta X_{t-i}+\varepsilon_t$, but the distinction is not crucial for our discussion.)
It is well known (see engle_granger1987,johansen_book) that, under technical conditions, the process $X_t$ is cointegrated if and only if $\Pi\ne 0$. In particular, testing for the absence of cointegration can be recast as testing the hypothesis $\Pi = 0$ in model (ref).
A widely used procedure for testing this hypothesis is the Johansen likelihood ratio test, introduced in johansen1988,johansen1991 and based on Gaussian maximum likelihood (see also Anderson). The method involves computing the eigenvalues $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_N$ of a certain matrix constructed from the observed data ${X_t}$. These eigenvalues lie in the interval $[0,1]$ and can be interpreted as squared sample canonical correlations between transformed lagged levels and first differences of the time series. We describe a slightly modified version of this procedure in Section (ref).
The statistic for testing the hypothesis $$H_0:\:{\rm rank}(\Pi)=0\qquad (\text{i.e.},\, \Pi\equiv0)$$ against the alternative $$ H(r):\quad {\rm rank}(\Pi)\in[1,r] $$ is based on the $r$ largest eigenvalues and takes the form: $$ \sum_{i=1}^r\ln(1-\lambda_i). $$ Under $H_0$ the $\lambda_i$'s tend to be small, making the statistic close to zero. Under the alternative $H(r)$, one expects some $\lambda_i$ to be close to one, which makes the sum more negative. Thus, the test rejects $H_0$ when the statistic is sufficiently negative.
The parameter $r$ is user-specified, and for fixed $N$ this gives rise to $N-1$ different tests corresponding to $r = 1, 2, \dots, N-1$. The asymptotic distribution of the statistic $\sum_{i=1}^r \ln(1 - \lambda_i)$ under $H_0$ as $T \to \infty$, with $N$ fixed was derived in johansen1988, johansen1991. It involves the eigenvalues of a certain matrix of It\^{o} integrals. This result underpins the classical cointegration testing procedure: compute the test statistic from data and compare it to the quantiles of its theoretical asymptotic distribution under $H_0$; if the observed value is smaller than the critical threshold, reject $H_0$ and conclude that the series $X_t$ exhibits cointegration.
In practice the Johansen procedure and its associated software are typically applied only when $N$ is small. The main reason is that the quality of the approximation based on the asymptotic distribution involving It\^{o} integrals deteriorates rapidly as $N$ increases. Specifically, the distribution of the statistic $\sum_{i=1}^r \ln(1 - \lambda_i)$ can deviate significantly from its theoretical large-$T$ limit even for moderate values of $N$ (for example, the case $N = 10$, $T = 100$ already yields poor performance). As a result, the likelihood ratio test tends to over-reject the null hypothesis $H_0$ in finite samples. See onatski_ecta, onatski_joe for a detailed theoretical analysis of this phenomenon.
This breakdown highlights the need to adapt the testing procedure and software for high-dimensional settings. In this work, we address this challenge by developing a new approach that is reliable for large $N$, building on recent theoretical results from bykhovskaya_gorin_1, bykhovskaya_gorin_k.
We now discuss a procedure implemented in {\fontseries{m}\fontseries{b}\selectfont Largevars} package, which is a modification of the construction of eigenvalues $\lambda_i$ used in the Johansen likelihood test. Starting with the data set $X_t$, $0\le t \le T$, we fix a number $r=1,2,\dots$ and perform the following steps.
Step 1 (Detrending). We de-trend and shift the data by defining
Step 2 (Cyclic indexing and regressor construction). We define cyclic indices modulo $T$: for any $a\in\mathbb Z$, set $$ a\mid T= a+ m T,\quad \text{where } m\in\mathbb Z\text{ is such that } a+ m T\in \{1,2,\dots,T\}. $$ Using this notation, we construct the following regressor matrices: $$ \tilde{Z}_{0t}=\Delta X_{t\mid T}\equiv\Delta X_{t},\quad\tilde{Z}_{kt}=\tilde X_{t-k+1\mid T},\quad\tilde{Z}_{1t}=(\Delta X_{t-1\mid T}^\top,\ldots,\Delta X_{t-k+1\mid T}^\top,1)^\top, \quad 1\le t \le T. $$ Here and below $^\top$ denotes matrix transposition. For each fixed $t$, the vector $\tilde{Z}_{1t}$ is a column of dimension $((k-1)N+1)\times 1$. The index $k$ in $\tilde{Z}_{kt}$ is used symbolically to reflect the VAR($k$) structure and does not refer to a specific numerical value --- this convention follows johansen1988, johansen1991.
We emphasize that due to the use of cyclic indices, values of $X_t$ at $t = 0, -1, \ldots$ are replaced by values at $t = T, T-1, \ldots$. However, when $k = 1$ (i.e., for a VAR($1$) model), no negative indices arise for $1 \le t \le T$, and the cyclic indexing becomes irrelevant.
Step 3 (Regression residuals). We compute the residuals from the regressions of $\tilde{Z}_{0t}$ and $\tilde{Z}_{kt}$ on $\tilde{Z}_{1t}$:
Step 4 (Canonical Correlations). Let $\tilde{R}_i$ be the $N \times T$ matrix whose columns are $\tilde{R}_{it}$ for $1 \le t \le T$, for $i = 0, k$. Define the cross-product matrices
and form the matrix
The eigenvalues $\tilde{\lambda}_1 \ge \ldots \ge \tilde{\lambda}_N$ of $\tilde{\mathcal{C}}$ represent the squared sample canonical correlations between $\tilde{R}_k$ and $\tilde{R}_0$. Equivalently, they solve the eigenvalue problem
Step 5 (Test Statistic). We construct the modified likelihood ratio statistic:
The subscript $(N,T)$ indicates that this version of the Johansen LR test is tailored for the high-dimensional regime where both $N$ and $T$ are large. After centering and scaling, the statistic $LR_{N,T}(r)$ is compared with suitable critical values to decide whether to reject the null hypothesis $H_0$. Heuristically, rejection corresponds to the case where the largest eigenvalues $\tilde{\lambda}_i$ are significantly large and well-separated from the rest.
We now describe the asymptotic distribution theory underlying the critical values used in our procedure. These formulas are based on the limiting behavior of the test statistic $LR_{N,T}(r)$ as both $N, T \to \infty$. The relevant asymptotics were developed in bykhovskaya_gorin_1, bykhovskaya_gorin_k, and we briefly recall the key results.
The limiting distribution involves a stochastic object known as the Airy$_1$ point process, denoted by $\{\mathfrak a_i\}_{i=1}^{\infty}$. This is a random, strictly decreasing sequence of real numbers: $\mathfrak a_1 > \mathfrak a_2 > \mathfrak a_3 > \dots$. We discuss this process in more detail in the next section.
Let $T$, $N$, and $k$ be such that $\frac{T}{N} > k+1$. Define the following constants:
$$ c_1\left(N,T\right)=\ln\left(1-\lambda_+\right), \qquad c_2\left(N,T\right)=-\frac{2^{2/3} \lambda_+^{2/3}}{(1-\lambda_+)^{1/3} (\lambda_+-\lambda_-)^{1/3}} \left(\mathfrak p+\mathfrak q\right)^{-2/3} <0. $$ Then, under appropriate assumptions, it follows from bykhovskaya_gorin_1 and bykhovskaya_gorin_k that
The practical applicability of this result depends on whether the theoretical assumptions hold for a given data set $X_t$. In Section (ref) we discuss model diagnostics that users can perform to assess the validity of the asymptotic approximation in real data.
To carry out the test in practice, one starts with the statistic $LR_{N,T}(r)$ defined in (ref). We recommend choosing small values of $r$ (e.g., $r = 1, 2,$ or $3$). The approximation in (ref) assumes that $r$ is fixed as $N, T \to \infty$ --- the rationale for this and its implications are discussed in detail in bykhovskaya_gorin_1.
The testing procedure then proceeds by computing the rescaled statistic: $$ \frac{LR_{N,T}(r)- r \cdot c_1(N,T)}{ N^{-2/3} c_2(N,T)} $$ and comparing it to the quantiles of the distribution $\sum_{i=1}^r \mathfrak a_i$. If the rescaled value exceeds the $\alpha$-quantile, we reject the null hypothesis of no cointegration at the $(1 - \alpha)$ significance level. The function \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{largevar()} in the {\fontseries{m}\fontseries{b}\selectfont Largevars} package implements this procedure.
The asymptotic formula (ref) shows that implementing our cointegration testing procedure requires knowledge of the distribution of the random variables $\sum_{i=1}^r \mathfrak a_i$. In this section we discuss how this distribution can be computed.
The Airy$_1$ point process is a random infinite sequence of real numbers $\mathfrak a_1 > \mathfrak a_2 > \mathfrak a_3 > \dots$ that can be defined via the following proposition:
The distribution of $\mathfrak a_1$ is known as the Tracy–Widom distribution $F_1$. tw1996 showed that the cumulative distribution function of $F_1$ can be expressed as the solution to a Painlevé differential equation. Numerical solutions to this equation were used by Bejan to compute highly accurate tables of quantiles of $F_1$. See also dieng2005distribution, bornemann2009numerical, trogdon2024computing for further numerical advances.
Several software packages incorporate precomputed tables of the Tracy–Widom distribution for practical use. For example, the {\fontseries{m}\fontseries{b}\selectfont RMTstat} package in \proglang{R} provides such functionality; see rmstatpackage.
When $r > 1$, much less is known about the distribution of $\sum_{i=1}^r \mathfrak a_i$, and it is unclear whether any of the approaches described in the previous paragraph remain applicable. Therefore, in the {\fontseries{m}\fontseries{b}\selectfont Largevars} package, we embed precomputed quantile tables for $\sum_{i=1}^r \mathfrak a_i$ obtained through direct simulation based on the definition in (ref).
The convergence rate in (ref) is of order $N^{-1/3}$, which is relatively slow. For instance, even with a large matrix of size $N=1000$, the approximation error remains around $N^{-1/3} = 0.1$. johnstone2012fast proposed computational techniques that accelerate convergence to $N^{-2/3}$ for $\mathfrak a_1$, but it is unknown whether such techniques yield similar improvements for $\sum_{i=1}^r \mathfrak a_i$ with $r > 1$. Moreover, testing this is nontrivial: while johnstone2012fast could compare their numerics against the known distribution of $\mathfrak a_1$, no such benchmark exists for the sum of the top $r$ points. Consequently, to ensure accurate quantile estimation via (ref), we opted for a very large matrix size: $N = 10^8$.
Using (ref) with $N = 10^8$ presents a computational challenge: no modern system can compute the eigenvalues of a $10^8 \times 10^8$ dense matrix. To circumvent this, we employ a numerical technique based on the tridiagonalization of the symmetric matrix $\frac{1}{2}(Y_N + Y_N^\top)$, following the approach of dumitriu_edelman. This reduces the problem to computing the eigenvalues of a real symmetric tridiagonal matrix of size $N \times N$:
where all entries on or above the diagonal are independent. Here, $\mathcal{N}(0,2)$ denotes a normal random variable with mean $0$ and variance $2$, and $\chi_\ell$ denotes the square root of a chi-squared random variable with $\ell$ degrees of freedom.
Directly computing the eigenvalues of the full matrix in (ref) for $N = 10^8$ remains computationally infeasible. However, one can instead consider the eigenvalues of its top-left $\sqrt{N} \times \sqrt{N}$ submatrix. Owing to the specific structure of (ref), the largest eigenvalues of the full $N \times N$ matrix and those of the $\sqrt{N} \times \sqrt{N}$ submatrix have the same asymptotic distribution; see edelman2005numerical and johnstone2021spin for theoretical justification.
We leverage this result in our simulations by performing $10^7$ Monte Carlo runs on symmetric tridiagonal random matrices of size $10^4 \times 10^4$, corresponding to the top-left corner of the tridiagonal matrix in (ref) with $N = 10^8$. After computing their eigenvalues, we rescale them according to the transformation in (ref) with $N = 10^8$. This yields approximations for the distribution of the individual $\mathfrak a_i$ values and, consequently, for the sums $\sum_{i=1}^r \mathfrak a_i$. We carried out this procedure for $r = 1, 2, \dots, 10$.
To assess the quality of our approximation, we do not run all $10^7$ simulations in a single batch. Instead, we perform $10^6$ simulations ten times using different initial random seeds. The average of the resulting sample quantiles provides our estimates for the quantiles of $\sum_{i=1}^r \mathfrak a_i$, while the standard deviation across the ten runs (not shown here) offers a measure of their reliability. The resulting quantile tables are embedded within the package and also presented separately in Section (ref). The standard deviations suggest that the error is at most $\pm 1$ in the third significant digit, i.e., $\pm 0.01$ for $r=1$ and $\pm 0.1$ for $r=10$. For $r=1$, our results closely match those of Bejan.
The simulation scripts were written in \proglang{MATLAB} and executed on a computing cluster provided by the Department of Economics at Duke University. Due to the substantial runtime, these simulations cannot be executed “on the fly” within the \proglang{R} package and also it limits our ability to further increase the precision of the quantile tables. For reproducibility, we include both the full script and a reduced version (with smaller $N$ and a fixed random seed), along with its output --- a low-precision version of the tables presented in Section (ref).
Finding higher-precision quantiles for $\sum_{i=1}^r \mathfrak a_i$ remains an open problem, whether by theoretical or numerical means. We hope that the tables in Section (ref) will help stimulate further interest in this question.
The procedure for cointegration testing relies on the validity of the approximation (ref) under the null hypothesis $H_0$ of no cointegration. This approximation, in turn, depends on theoretical assumptions made in its derivation. A natural concern for users is whether these assumptions are reasonable for a given dataset $X_t$, $0 \le t \le T$.
bykhovskaya_gorin_k provides a mathematical justification for (ref) under the setting where both $T$ and $N$ are large, with $T/N \in (k+1, \infty)$ bounded away from the endpoints. The proof assumes Gaussian errors $\varepsilon_t$ in (ref) and imposes the restriction $\Gamma_1 = \Gamma_2 = \dots = \Gamma_{k-1} = 0$. The discussion following the theorem in the same article argues that the approximation should remain valid when the $\Gamma_i$ matrices are of low rank, and bykhovskaya_gorin_1 further argues that the Gaussianity assumption on $\varepsilon_t$ is likely not essential.
Moreover, onatski_ecta and bykhovskaya_gorin_k demonstrate that under these relaxed assumptions an additional result holds --- independent of whether cointegration is present (i.e., whether $\Pi = 0$ in (ref) or not), as long as $\Pi$ remains of small rank. Specifically, they show that the histogram of all eigenvalues $\tilde{\lambda}_1, \dots, \tilde{\lambda}_N$ from (ref) converges to a deterministic limiting distribution.
In more detail, the Wachter distribution is a probability distribution on the interval $[0,1]$ that depends on two parameters $\mathfrak p > 1$ and $\mathfrak q > 1$, and is defined by the density
where the support $[\lambda_-, \lambda_+] \subset (0,1)$ is given by
onatski_ecta and bykhovskaya_gorin_k show that if the parameters are chosen according to (ref), then the empirical distribution of the eigenvalues $\tilde{\lambda}_i$ from (ref) converges to the Wachter distribution:
This convergence provides a practical check for the applicability of our procedure. If, for a given data set $X_t$, the histogram of eigenvalues $\tilde{\lambda}_i$ resembles the shape of the Wachter distribution (possibly excluding a few outlying values, which may correspond to cointegration), then it is reasonable to trust the assumptions behind our cointegration test. If not, the modeling assumptions are likely violated, and the test results should not be used.
To facilitate this diagnostic step, the function \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{largevar()} in our package includes an option to plot the empirical histogram of eigenvalues along with the Wachter density.
The latest version of the {\fontseries{m}\fontseries{b}\selectfont Largevars} package can always be found on Github and installed using the \verb|devtools| R package:
The latest stable version from CRAN can be installed by \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{install.packages("Largevars")}. Help for using the functions in the package can be called by running ?? \verb|function name| . The empirical example in Section (ref) of this paper can provide further guidance.
\bgroup\@makeother_\@makeother\ \@makeother\$\@codex{largevar()} is the main function in the package that implements the cointegration test for high-dimensional VARs.
The function \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{largevar()} operates according to the steps laid in out in Section (ref). The test statistic is formed based on the $r$ largest eigenvalues. Any value of $r$ can be used to reject the hypothesis $H_0$ of no cointegration, and the user can try different options. We recommend small values such as $r=1,2,3$, see bykhovskaya_gorin_1 for the detailed discussion.
\bgroup\@makeother_\@makeother\ \@makeother\$\@codex{largevar()} returns a list object that contains the test statistic, a statistical table with a subset of theoretical quantiles ($q=0.90,\,0.95,\,0.97,\,0.99$) presented for $r=1$ to $r=10$, the decision about $H_{0}$ at the significance level specified by the user, and the p-value. These can be accessed by \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{list\$statistic} (numeric value) \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{list\$significance_test\$significance_table} (numeric matrix), \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{list\$significance_test\$boolean_decision} (numeric value of $0$ or $1$, where $1$ means “reject”), and \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{list\$significance_test\$p_value} (numeric value), respectively.
The simulations for the quantiles of the limiting distribution were conducted for $r=1$ to $r=10$ values. For this reason, $p$ values are accessible at inputs $r=1$ to $r=10$ only. For larger $r$ inputs, the function returns the test statistic but not the $p$ value and not the decision about $H_{0}$ at the significance level specified by the user.
To access the test quantile tables for the partial sums of the Airy$_1$ point process, discussed in Section (ref) and in Section (ref), the user can call the \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{quantile_tables()} function. Quantile tables are available for $r=1$ to $r=10$. The function returns a numeric matrix, where the $0.ab$ quantile corresponds to the row $0.a$ and the column $b$. For example:
\bgroup\@makeother_\@makeother\ \@makeother\$\@codex{sim_function()} is an auxiliary function that allows the user to calculate an empirical $p$ value based on a simulation of the data generating process $\widehat H_{0}$ stated in of bykhovskaya_gorin_k; equivalently this is the model (ref) based on mean $0$ Gaussian errors $\varepsilon_t$ and $\Gamma_1=\dots=\Gamma_{k-1}=0$. This function should be used only for quick approximate assessments, as precise computation of the distribution of the test statistic requires a very large number of simulations, as discussed in Section (ref).
The function \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{sim_function()} runs the cointegration test (following the steps of Section (ref)) on simulated data generated under the evolution (ref) based on mean $0$ Gaussian errors $\varepsilon_t$ and $\Gamma_1=\dots=\Gamma_{k-1}=0$ and calculates the empirical $p$ value based on the test statistic (ref) corresponding to \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{r} specified by the user. The empirical $p$ value is defined as the fraction of realizations larger than the specified \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{stat_value}. For comparison purposes, it is advised to specify the same parameters \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{k} and \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{r} as one expects to use for the run of \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{largevar()} for the desired data set.
\bgroup\@makeother_\@makeother\ \@makeother\$\@codex{sim_function()} returns a list object that contains the simulation values, the empirical $p$ value and a histogram of the distribution of simulated test statistic values (which is an approximation of the probability distribution of the test statistic).
This section provides two examples of the usage of the package. Section (ref) replicates the S$\&$P$100$ example from bykhovskaya_gorin_1, bykhovskaya_gorin_k, while Section (ref) uses simulated data. Both examples include the code, which can be copied into \proglang{R}.
We use logarithms of weekly adjusted closing prices of assets in the S&P100 over ten years (01.01.2010--01.01.2020), which gives us $\tau=522$ observations across time. The S&P100 includes 101 stocks, with Google having two classes of stocks. We use 92 of those stocks, those for which data were available for our chosen time period. Only one of Google’s two listed stocks is kept in the sample. Therefore, $N = 92$, $T = 521$ and $T/N \approx 5.66$. We obtained the raw data from Yahoo! Finance and made the sample available in the “data” folder of the package for convenient data loading:
We first make necessary transformations, then convert to a numeric matrix to match function requirements:
There is documentation available for the following function which can be called using
The following code conducts the cointegration test and displays its results:
Since we set \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{plot_output=TRUE}, we obtain a histogram of eigenvalues solving (ref), as shown in Figure (ref). The resemblance of the histogram with the theoretical curve is very good and we expect that our cointegration test should be applicable to this data set. The remaining output of \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{largevar()} is displayed in the console as:
The decision $0$ means that we do not reject $H_0$, i.e., it is likely that the data set has no cointegration. If we want to individually access certain values from the output list, we can do it by referencing the elements of the list:
We can further compare the exact $p$ value (which was outputted by \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{largevar()}) with a sample value obtained trough $1000$ simulations, by running \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{sim_function()}:
As we see, the empirical $p$ value $0.247$ is close to the previous output $0.23$, but there is a small mismatch, as expected.
We also present an example based on simulated data that users can replicate. The code below generates VAR($2$) with $N = 100$, $T = 1500$, and
where $\Delta X_{it}:=X_{it}- X_{it-1}$. The process is initialized by vectors $X_0,\,X_{-1}$ with independent standard normal coordinates. The data generating process (ref) corresponds to a matrix $\Pi$ of rank $2$: $\Pi$ has two nonzero and linearly independent rows. To be more precise, the coefficient matrices in Eq. (ref) correspond to $N-2$ unit root and $2$ stationary components. We create a data set based on evolution (ref) and Gaussian errors, and then run the cointegration test on it, as follows.
The following code constructs matrices $\Pi$ and $\Gamma$ that by construction create two separate cointegrating systems, the first and second, and the fourth and fifth time series:
The initialization of the time series by setting all the below values:
The development of the system up to $T$ is calculated, starting with changes $dX$ for each $t$:
Finally, we conduct the cointegration test and display the results:
Since the decision is $1$, we reject $H_0$ and conclude that the data set has cointegration. The rejection is in line with the largest eigenvalue being significantly to the right from $\lambda_+$. The separation of eigenvalues is clearly visible in the histogram output, see Figure (ref).
If we want to take a look at how the significance of our test statistics vary across different choices of $r$, we can call the table below. The $p$ values for our test statistics stay below $0.01$.
We include the tables discussed in Section (ref), both inside the package and in this section. These tables are used inside \bgroup\@makeother_\@makeother\ \@makeother\$\@codex{largevar()} to obtain the quantiles. The tables below present our simulation results. The $0.ab$ quantile in each table corresponds to the row $0.a$ and the column $b$. The standard deviations of our results suggest that the error is at most $\pm 1$ in the third digit of the elements of the Airy$_1$ sequence, meaning that the error is $\pm 0.01$ for $r=1$ and $\pm 0.1$ for $r=10$.