Inference in heavy-tailed non-stationary multivariate time series

Inference in heavy-tailed non-stationary multivariate time series

\address{Matteo Barigozzi, University of Bologna, [email removed]} \address{Giuseppe Cavaliere, University of Bologna, [email removed]} \address{Lorenzo Trapani, University of Nottingham, [email removed]} \subjclass{Primary 62M10; Secondary 62G20}

\doublespacing

Introduction

Since the seminal works by englegranger, SW1988 and johansen1991, investigating the presence and number $m$ of common stochastic trends has become an essential step in the analysis of multivariate time series which are non-stationary over time. Inference on $m$ is of great importance on its own, as it has a \textquotedblleft natural\textquotedblright\ interpretation in many applications: for example, it can provide the number of non-stationary factors in Nelson-Siegel type term structure models; it allows to assess the presence of (long-run) integration among financial market (kasa), or, in the context of asset pricing, the number of portfolios that allow hedging of long-run risks (alexander). In terms of theory, available estimators are based either on sequential testing (see e.g. boswijk2015improved), or on information criteria (see aznar and qu2007modified).

With few exceptions, a standard feature of existing techniques is that they strongly rely on the assumption that some moments of the data (usually at least the second, or even the fourth) exist. This assumption, however, often lacks of empirical support. In fact, data exhibiting heavy tails, which do not have finite second (or even first) moment, are often encountered in many areas: macroeconomics (see e.g. ibragimov), finance ( rachev2003handbook), urban studies (gabaix); and also insurance, actuarial science, telecommunication network traffic and meteorology (see e.g. the review by embrechts).

Violation of the moment assumptions may result in (possibly severe) incorrect determination of the number of common trends - see e.g. the simulations in caner and the empirical evidence in falk. Unfortunately, contributions which explicitly deal with inference on common stochastic trends under infinite variance are rare. caner derives the asymptotic distribution of Johansen's trace test under infinite variance and shows that it depends on the (unknown) tail index of the data. ling study the (non-standard) rate of convergence of estimators in non-stationary VAR models and show that the limiting distributions depend on the tail index in a non-trivial fashion; similar results are also found in paulauskas and fasen. In the univariate case, jach2004subsampling and cavaliere2018unit show that suitable (respectively, $M$-out-of-$N$ and wild) bootstrap approaches could be used to test whether data are driven by a stochastic trend; knowledge of the tail index is not needed, but extensions of these bootstrap approaches to multiple time series are not available, and are likely to be extremely difficult to develop and implement. Distribution-free approaches could also help overcome this difficulty. hallin2 (see also hallin1) apply the rank transformation to the VECM residuals, obtaining nuisance-free statistics, although this approach requires the correct specification of the VECM.

Contribution

Our key contribution is the estimator of the number $m$ of common trends of an $N$-variate time series in the possible presence of heavy tails. Crucially, our procedure does not require any a priori knowledge as to whether the variance is finite or not, or as to how many moments exist, thus avoiding having to estimate any nuisance parameters or even pre-testing for (in)finite moments. We also show that, in contrast to most of the literature on time series with heavy tails, our methodology also applies to time series with heterogeneous innovations. Specifically, we allow for changes in the scale of the innovations of a very general form, which covers, e.g., multiple shifts and smooth scale changes. As far as we are aware, this paper is the first one where heterogeneity in the scale is allowed under infinite variance.

Details are spelt out in the remainder of the paper; here, we offer a heuristic preview of how the methodology works. The starting point of our analysis is a novel result concerning the properties of the sample second moment matrix of the data in levels. Specifically, we show that the $m $ largest eigenvalues of the matrix diverge to positive infinity, as the sample size $T $ passes to infinity, faster than the remaining eigenvalues by a factor (almost) equal to $T$. Importantly, this result always holds, irrespective of the variance of the innovations being finite or infinite. Building on this, for each eigenvalue we construct a statistic which diverges to infinity under the null that the eigenvalue is diverging at a \textquotedblleft fast\textquotedblright\ rate, and drifts to zero under the alternative that the relevant eigenvalue diverges at a \textquotedblleft slow\textquotedblright\ rate. Although the limiting distribution of our statistic is bound to depend on nuisance parameters such as the tail index, the relative rate of divergence between the null and the alternative does not depend on any nuisance parameters. Therefore, in order to construct a test, we can rely on rates only, and randomise our statistic using a similar approach to bandi2014. Thence, our estimator of $m$ is based on running the tests sequentially; in this respect, it mimics the well-known sequential procedure advocated in johansen1991 for the determination of the rank of a cointegrated system. Our approach has some commonalities with the literature on factor models too, where the spectrum of the covariance matrix of the data is employed to estimate the number of common factors. zyr18 propose a methodology to determine the rank of cointegration in a possibly fractionally integrated system with finite second moments, based on the eigenvalues of (functions of) the sample autocovariance matrices (see also tu2020error).

Our methodology has at least four desirable features. Firstly, as mentioned above, our technique can be applied to data with infinite variance (and even infinite expectation), with no need to know this a priori. Secondly, our procedure does not require at any stage the estimation of the tail index of a distribution, which is a notoriously delicate issue. Thirdly, our procedure does not require the correct specification of the lag structure of the underlying VECM model, and it is therefore completely robust to misspecification of the dynamics. Finally, our results are based only on first order asymptotics, i.e. rates, which makes our procedure extremely easy to implement in practice.

All these advantages are illustrated via Monte Carlo simulation, where we show that even in samples of moderate size our method performs extremely well, and through four empirical applications, which cover commodity prices, interest rates data, purchasing power parity in the long run, and cryptocurrency markets.

Structure of the paper and notation

The remainder of the paper is organised as follows. Assumptions and asymptotics are in Section (ref). The main results on the number of common trends and the estimation of the common trends (and associated loadings) are presented in Section (ref), and extensions are collected in Section (ref). We provide Monte Carlo evidence in Section (ref), and we validate our methodology through four real data applications in Section (ref). Section (ref) concludes. Further simulations, and all technical lemmas and proofs, are relegated to the Supplement.

Throughout the paper, we make use of the following notation. For a given matrix $A\in \mathbb{R} ^{n\times m}$, we denote its element in position $\left( i,j\right) $ as $ A_{i,j}$; we also let let $\lambda ^{\left( j\right) }\left( A\right) $ denote the $j$-th largest eigenvalue of $A$. We denote with $c_{0},c_{1},...$ positive, finite constants whose value can change from line to line. We denote the $k$-iterated logarithm of $x$ (truncated at zero) as $\ln _{k}x$ - e.g. $\ln _{2}x=\max \left\{ \ln \ln x,0\right\} $. The backshift operator is denoted as $L$.

Theory

Consider an $N$-dimensional vector $y_{t}$ with $MA\left( \infty \right) $ representation

where $C\left( L\right) =\sum_{j=0}^{\infty }C_{j}L^{j}$. We assume that $y_{0}=0$ and $\varepsilon _{t}=0$ for $t\leq 0$, for simplicity and with no loss of generality; inspecting the proofs, all the results derived here can be extended to more general assumptions concerning e.g. the initial value $y_{0}$. Standard arguments allow to represent ((ref)) as

where $C=\sum_{j=0}^{\infty }C_{j}$, $C^{\ast }\left( L\right) =\sum_{j=0}^{\infty }C_{j}^{\ast }L^{j}$ and $C_{j}^{\ast }=\sum_{k=j+1}^{\infty }C_{k}$. Equation ((ref)) is derived from the multivariate Beveridge-Nelson decomposition of the filter $C\left( L\right)$ (see watson1994vector).

We assume that the $N\times N$ matrix $C$ can have reduced rank, say $m$. This corresponds to assuming that the long-run behavior of the N-dimensional vector $y_{t}$ is driven by m non-stationary common factors.

Part (ii) of the assumption requires that the $MA$ coefficients $ C_{j}$ decline geometrically. This assumption is similar to Assumption 1 in caner, where the $C_{j}$s are assumed to decline at a rate which increases as the tail index of the innovations $\varepsilon _{t}$ decreases. The case $m=0$ correspond to $y_{t}$ being (asymptotically) strictly stationary. Assumption (ref) is also implied by Assumption 2.1 in ling, where a finite-order VAR model under the classic $I\left( 1\right) $ conditions stated e.g. in reisnel is considered.

On account of the possible rank reduction of $C$, a different formulation of ((ref)) may be helpful. It is always possible to write $C=PQ$, where $P$ and $Q$ are full rank matrices of dimensions $N\times m$ and $m\times N$ respectively. Defining the $m$-dimensional process $x_{t}=Q\sum_{s=1}^{t} \varepsilon _{s}$, and using the short-hand notation $u_{t}=C^{\ast }\left( L\right) \varepsilon _{t}$, then ((ref)) can be given the following common trend representation

We now make some assumptions on the error term $\varepsilon _{t}$.

As is generally the case in the analysis of time series with possibly infinite variance, we assume that $\varepsilon _{t}$ is i.i.d. Part (ii) of the assumption implicitly states that the vector $ \varepsilon _{t}$ has a multivariate distribution which belongs to the domain of attraction of a strictly stable, multivariate law (see Theorem 2.1.5(a) in taqqu) with common tail index $\eta $, so that linear combinations of them can be constructed. The assumption implies that, when $ \eta <2$, $E\left\vert \varepsilon _{i,t}\right\vert ^{p}<\infty $ for all $ 0\leq p<\eta $, whereas $E\left\vert \varepsilon _{i,t}\right\vert ^{\eta }=\infty $ (see petrov1974). Also, by Property 1.2.15 in taqqu, it holds that, as $x\rightarrow \infty $

where $L\left( x\right) $ is a slowly varying function in the sense of Karamata (see seneta for a review), and $c_{l,1},c_{l,2}\geq 0$, $ c_{l,1}+c_{l,2}>0$. The condition that $G$ is strictly stable entails that $c_{l,1}=c_{l,2}$ when $\eta =1$, thus ruling out asymmetry in that case (see Property 1.2.8 in taqqu).

Asymptotics

Define

We report a set of novel results for the eigenvalues of $S_{11}$ and $S_{00}$ , which we require for the construction of the test statistics.

Proposition (ref) states that the first $m$\ eigenvalues of $ S_{11}$ diverge at a faster rate than the other ones (faster by an order of \textquotedblleft almost\textquotedblright\ $T$), thus entailing that the spectrum of $S_{11}$ has a \textquotedblleft spiked\textquotedblright\ structure.

In order to study $S_{00}$ and its eigenvalues, we need the following assumption, which complements Assumption (ref)(ii).

The integral Lipschitz condition in Assumption (ref) is a technical requirement needed for $\Delta y_{t}$ to be strong mixing with geometrically declining mixing numbers, and it is a standard requirement in this literature (see e.g. withers and phamtran).

Similarly to Proposition (ref), Proposition (ref) provides bounds for the spectrum of $S_{00}$. Part ((ref)) has been shown in trapaniLIL, and it is a Chover-type Law of the Iterated Logarithm (chover). As shown in trapaniLIL, the bound in ((ref)) is almost sharp. The lower bound implied in ( (ref)) is also almost sharp. \newline

The spectrum of $S_{11}$ (and, in particular, the different rates of divergence of its eigenvalues) can -- in principle -- be employed in order to determine $m$. However, $S_{11}$ has two main problems which make it unsuitable for direct usage. First, by Proposition (ref), the spectrum of $S_{11}$ depends on the nuisance parameter $\eta $. Also, the matrix $S_{11}$ and, consequently, its spectrum depend on the unit of measurement of the data, and thus they are not scale-free. In order to construct scale-free and nuisance-free statistics, we propose to rescale $ S_{11}$ by $S_{00}$. Proposition (ref) ensures that this is possible: by equation ((ref)), the inverse of $S_{00}$ cannot diverge too fast, and therefore the spectrum of the matrix $ S_{00}^{-1}S_{11} $ should still have $m$ eigenvalues that diverge at a faster rate than the others. This is shown in the next theorem.

Theorem (ref) states that the spectrum of $S_{00}^{-1}S_{11}$ has a similar structure to the spectrum of $S_{11}$: the first $m$ eigenvalues are spiked and their rate of divergence is faster than that of the remaining eigenvalues by a factor of almost $T$. More importantly, by normalising $ S_{11}$ by $S_{00}$, dependence on $\eta $ is relegated to the slowly-varying (logarithmic) terms. In essence, apart from the slowly varying sequences, equations ((ref)) and ((ref)) imply that the rates of divergence of the eigenvalues of $S_{00}^{-1}S_{11}$ are of order (arbitrarily close to) $O\left( T\right) $ for the spiked eigenvalues, and (arbitrarily close to) $O\left( 1\right) $ for the other ones. This is the key property of $\lambda ^{\left( j\right) }\left( S_{00}^{-1}S_{11}\right) $: dividing $S_{11}$ by $S_{00}$\ washes out the impact of the tail index $\eta $, which essentially does not play any role in determining the divergence or not of $\lambda ^{\left( j\right) }\left( S_{00}^{-1}S_{11}\right) $. Although this result pertains to first order asymptotics (i.e., rates), it is possible to find an analogy between the result in Theorem (ref) and approaches based on eliminating nuisance parameters using self-normalisation (see e.g. shao).

Inference on the common trends

In this section we collect our main results about estimation and inference on the common trends in the possible presence of heavy tails. In Section (ref), we report a novel one-shot test about the (minimum) number of common trends. Then, in Section (ref) we introduce a sequential procedure for the determination of the number of common trends. Estimation of the common trends and associated factor loading is presented in Section (ref).

Testing hypotheses on the number of common trends

We consider testing hypotheses on the number of common trends $m$. The tests we propose herein will form the basis of our sequential procedure for the determination of the number of common trends -- see Section (ref). Specifically, we consider the null hypothesis

where $j\in \{1,...,N\}$ is a (user-chosen) lower bound on the number of common trends. For instance, testing whether there is at least one common trend would entail setting $j=N-1$; a test of non-stationarity against the alternative of strict stationarity corresponds to $j=1$.

Based on Theorem (ref), we propose to use

where $\kappa \in \left( 0,1\right) $; criteria for the choice of $\kappa $ in applications are discussed in Section (ref). Importantly, the $ T^{-\kappa }$ term in ((ref)) is used in order to exploit the discrepancy in the rates of divergence of the $\lambda ^{\left( j\right) }\left( S_{00}^{-1}S_{11}\right) $ under $H_{0}$ and under $H_{A}$. In particular, it ensures that $T^{-\kappa }\lambda ^{\left( j\right) }\left( S_{00}^{-1}S_{11}\right) $ drifts to zero under $H_{A}$ (i.e., whenever $j>m$ ), whereas it still passes to infinity under $H_{0}$ (i.e., when $j\leq m$). According to ((ref)), this only requires a very small value of $ \kappa $, which would also allow $\lambda ^{\left( j\right) }\left( S_{00}^{-1}S_{11}\right) $ to diverge at a rate close to $T$ under $H_{0}$. On account of Theorem (ref), it holds that $P(\omega :\lim_{T\rightarrow \infty }\phi _{T}^{\left( j\right) }=\infty )=1$ for $ 0\leq j\leq m$; hence, we can assume that under the null that $m\geq j$, it holds that $\lim_{T\rightarrow \infty }\phi _{T}^{\left( j\right) }=\infty $ . Conversely, we have $P(\omega :\lim_{T\rightarrow \infty }\phi _{T}^{\left( j\right) }=0)=1$ for $j>m$, which entails that under the alternative that $j>m$, we have $\lim_{T\rightarrow \infty }\phi _{T}^{\left( j\right) }=0$. In essence, $\phi _{T}^{\left( j\right) }$ diverges to positive infinity, or converges (to zero), according to whether $ \lambda ^{\left( j\right) }\left( S_{00}^{-1}S_{11}\right) $ is \textquotedblleft large\textquotedblright\ or \textquotedblleft small\textquotedblright .

Since the limiting law of $\phi _{T}^{\left( j\right) }$ under the null is unknown, we propose a randomised version of it. The construction of the test statistic is based on the following three step algorithm, which requires a user-chosen weight function $F\left( \cdot \right) $ with support $ U\subseteq \mathbb{R}$.

Step 1 Generate an artificial sample $\{\xi _{i}^{\left( j\right) },1\leq i\leq M\}$, with $\xi _{i}^{\left( j\right) }\sim $ i.i.d.$N\left( 0,1\right) $, independent of the original data.\newline Step 2 For each $u\in U$, define the Bernoulli sequence $ \zeta _{i}^{\left( j\right) }\left( u\right) =I(\phi _{T}^{\left( j\right) }\xi _{i}^{\left( j\right) }\leq u)$, and let

Step 3 Compute

where $F\left( \cdot \right) $ is the user-chosen weight function.

In Step 2, the binary variable $\zeta _{i}^{\left( j\right) }\left( u\right) $ is created for several values of $u \in U$, and in Step 3, the resulting statistics $\theta _{T,M}^{\left( j\right) }\left( u\right) $ are averaged across $u$, through the weight function $F\left( \cdot \right) $, thus attenuating the dependence of the test statistic on an arbitrary value $u$. The following assumption characterizes $F\left( \cdot \right) $.

A possible choice for $F\left( \cdot \right) $ could be a distribution function with finite second moment; note that part (ii) of the assumption is trivially satisfied if $U$ is bounded. Possible examples include a Rademacher distribution, based on choosing $U=\left\{ -c,c\right\} $ with $F\left( c\right) =F\left( -c\right) =1/2$, or the standard normal distribution function.

Let $P^{\ast }$ denote the conditional probability with respect to the original sample; we use the notation \textquotedblleft $\overset{D^{\ast }}{ \rightarrow }$\textquotedblright\ and \textquotedblleft $\overset{P^{\ast }}{ \rightarrow }$\textquotedblright\ to define, respectively, conditional convergence in distribution and in probability according to $P^{\ast }$.

Theorem (ref) provides the limiting behaviour of $\Theta _{T,M}^{\left( j\right) }$, also illustrating the impact of $M$ on the size and power trade-off. According to ((ref)), the larger $M$ the higher the power. Conversely, upon inspecting the proof, it emerges that $\theta _{T,M}^{\left( j\right) }\left( u\right) $ contains a non-centrality parameter of order $O\left( M^{1/2}\exp \left( -T^{1-\kappa -\epsilon }\right) \right) $, whence the upper bound in ((ref)).

The one-shot test developed in this section has at least three advantages compared to existing methods. First, our approach can also be implement to check (asymptotic) strict stationarity. Indeed, running the test for $j=1$ corresponds to the null hypothesis that the data are driven by at least one common trend; rejection supports the alternative of stationarity. Second, running the test with $j=N$ corresponds to the null hypothesis that the $N$ variables do not cointegrate, thus offering a test for the null of no cointegration against the alternative of (at least) one cointegrating relation. Finally, we point out a further advantage over the well-known method of johansen1991. Johansen's likelihood ratio test allows to test the null of rank $R$ (i.e., of $m=N-R$ common trends), where $R$ is user-chosen, versus the alternative of rank greater than $R$ (i.e., less than $N-R$ common trends). However, while the limiting distribution under the null is well-known, if the true rank is lower than $R$, then the limiting distribution is different (see bernstein2019). Hence, Johansen's test should be use only if the practitioner knows that the rank cannot be lower than $R$. In contrast, since the null hypothesis is formulated as a minimum bound on the number of common trends, our test does not have this drawback.

A final remark on the test is in order. Letting $0<\alpha <1$ denote the nominal level of the test, and defining $c_{\alpha }$ such that $P\left( \chi _{1}^{2}>c_{\alpha }\right) =\alpha $, an immediate consequence of the theorem is that under $H_{A}$ it holds that $\lim_{\min \left( T,M\right) \rightarrow \infty }P^{\ast }(\Theta _{T,M}^{\left( j\right) }>c_{\alpha })=1 $ for almost all realisations of $\left\{ \varepsilon _{t},0<t<\infty \right\} $ - i.e. the test is consistent under the alternative. Conversely, under $H_{0}$ we have, for almost all realisations of $\left\{ \varepsilon _{t},0<t<\infty \right\} $

As noted also in HT16, our test is constructed using a randomisation which does not vanish asymptotically, and therefore the asymptotics of $ \Theta _{T,M}^{\left( j\right) }$ is driven by the added randomness. Thus, different researchers using the same data will obtain different values of $ \Theta _{T,M}^{\left( j\right) }$ and, consequently, different $p$-values; indeed, if an infinite number of researchers were to carry out the test, the $p$-values would follow a uniform distribution on $\left[ 0,1\right] $. This is a well-known feature of randomised tests. In order to ameliorate this, HT16 suggest to compute $\Theta _{T,M}^{\left( j\right) }$ $S$ times, at each time $s$ using an independent sequence $\left\{ \xi _{i,s}^{\left( j\right) }\right\} $ for $1\leq j\leq M$ and $1\leq s\leq S$, thence defining

Based on standard arguments, under $H_{0}$ the LIL\ yields

Hence, a \textquotedblleft strong rule\textquotedblright\ to decide in favour of $H_{0}$ is

Decisions made on the grounds of ((ref)) have vanishing probabilities of Type I and Type II errors, and are the same for all researchers: having $ S\rightarrow \infty $ washes out the added randomness. No restrictions are required on $S$ as it passes to infinity, so it can be arbitrarily large.

Determining $m$

In order to determine the number of common trends $m$, we propose to cast the individual one-shot tests discussed above in a sequential procedure, where different values $j=1,2,...$ for $m$ are tested sequentially (note that the individual tests must be based on artificial random samples independent across $j$, see below). \newline The estimator of $m$ (say, $\widehat{m}$) is the output of the following algorithm:

Algorithm 1.

Step 1 Run the test for $H_{0}:m\geq 1$ based on $\Theta _{T,M}^{\left( 1\right) }$. If the null is rejected, set $\widehat{m}=0$\ and stop, otherwise go to the next step.

Step 2 Starting from $j=2$, run the test for $H_{0}:m\geq j$ based on $\Theta _{T,M}^{\left( j\right) }$, constructed using an artificial sample $\{\xi _{i}^{\left( j\right) }\}_{i=1}^{M}$ generated independently of $\{\xi _{i}^{\left( 1\right) }\}_{i=1}^{M},$ $...,$ $\{\xi _{i}^{\left( j-1\right) }\}_{i=1}^{M}$. If the null is rejected, set $ \widehat{m}=j-1$\ and stop; otherwise, if $j=N$, set $\widehat{m}=N$; otherwise, increase $j$ and repeat Step 2.

Consistency of the proposed procedure is presented in the next theorem.

Theorem (ref) states that $\widehat{m}$ is consistent, as long as the nominal level $\alpha $\ of the individual tests is chosen so as to drift to zero: no specific rates are required. This can be better understood upon inspecting the proof of the theorem: letting $\alpha $ denote the level of each \ individual test, in ((ref)), we show that, as $\min \left( T,M\right) \rightarrow \infty $, $P^{\ast }(\widehat{m}=m)=\left( 1-\alpha \right) ^{N-m}$, whence the requirement $c_{\alpha }\rightarrow \infty $, which entails$\ \alpha \rightarrow 0$. This can also be read in conjunction with Johansen's procedure (johansen1991) and its bootstrap implementations (cavaliereECTA), whose outcome is an estimate of $m$, say $\widetilde{m}$, such that, asymptotically, $P( \widetilde{m}=m)=1-\alpha $ for a given nominal value $\alpha $ for the individual tests. On account of ((ref)), in our case choosing a fixed nominal level $\alpha $ would yield, as mentioned above, $P^{\ast }( \widehat{m}=m)=\left( 1-\alpha \right) ^{N-m}$, which depends on the unknown $m$ and is, for $m>1$, worse than Johansen's procedure. On the other hand, in our procedure the individual tests are independent (conditional on the sample). Thus, in order to match the $1-\alpha $ probability of selecting the true $m$, one can use a Bonferroni correction with $\alpha /N$ as nominal level for each test, rather than $\alpha $. In this case, the same calculations as in the proof of Theorem (ref) (and Bernoulli's inequality) yield

Finally, as an alternative to Bonferroni correction, one could consider the following top-down algorithm.

Algorithm 2.

Step 1 Run the test for $H_{0}:m=N$ based on $\Theta _{T,M}^{\left( N\right) }$. If the null is not rejected, set $\widetilde{m} =N $\ and stop, otherwise go to the next step.

Step 2 Starting from $j=N-1$, run the test for $ H_{0}:m\geq j$ based on $\Theta _{T,M}^{\left( j\right) }$, constructed using an artificial sample $\{\xi _{i}^{\left( j\right) }\}_{i=1}^{M}$ generated independently of $\{\xi _{i}^{\left( 1\right) }\}_{i=1}^{M},$ $ ..., $ $\{\xi _{i}^{\left( j-1\right) }\}_{i=1}^{M}$. If the null is rejected, set $\widetilde{m}$\ and stop; otherwise, if $j=1$, set $ \widetilde{m}$; otherwise, decrease $j$ and repeat Step 2.

The proof of the consistency of $\widetilde{m}$, and some Monte Carlo evidence (showing that the performance of $\widetilde{m}$ is virtually indistinguishable from that of $\widehat{m}$), are in the Supplement.

Estimation of the common trends

Recall the common trend representation provided in ((ref)),

After determining $m$, it is possible to estimate the non-stationary common stochastic trends $x_{t}$ by using Principal Components (PC), in a similar fashion to penaponcela and zyr18.

In particular, let $\widehat{\upsilon }_{j}$ denote the eigenvector corresponding to the $j$-th largest eigenvalue of $S_{11}$ under the orthonormalisation restrictions $\left\Vert \widehat{\upsilon } _{j}\right\Vert =1$ and $\widehat{\upsilon }_{i}^{\prime }\widehat{\upsilon } _{j}=0$ for all $i\neq j$. Then, defining $\widehat{P}=(\widehat{\upsilon } _{1},...,\widehat{\upsilon }_{m})$, the estimator of the common trends $ x_{t} $ is

The next theorem provides the consistency (up to a rotation) of the estimators of $P$ and $x_{t}$. Interestingly, the convergence rate of $ \widehat{P}$ is not affected by the tail index.

Theorem (ref) states that both $\widehat{P}$ and $\widehat{x}_{t}$ are consistent estimators of $P$ and $x_{t}$ -- up to an invertible linear transformation, since it is only possible to provide a consistent estimate of the eigenspace, as opposed to the individual eigenvectors. Equation ((ref)) states that the estimated loadings\ $\widehat{P}$ provide a superconsistent estimator of (a linear combination of the columns of) $P$. This result, which is the same as in the case of finite variance, is a consequence of the fact that $x_{t}$ is an \textquotedblleft integrated\textquotedblright\ process, and it is related to the eigen-gap found in Proposition (ref). Equation ((ref)) could also be read in conjunction with the literature on large factor models, where -- contrary to our case -- it is required that $N\rightarrow \infty $. In that context, bai04 obtains the same result as in ((ref)) albeit for the case of finite variance: thus, in the presence of integrated processes, the PC\ estimator is always superconsistent, irrespective of $N$ passing to infinity or not.

According to ((ref)), $\widehat{x}_{t}$ also is a consistent estimator of the space spanned by $x_{t}$. The \textquotedblleft noise\textquotedblright\ component does not drift to zero and, when $\eta <1$ , it may even diverge; however, the \textquotedblleft signal\textquotedblright\ $x_{t}$ is of order $O_{P}\left( t^{1/\eta }\right) $, and therefore it dominates the estimation error (in fact, when $ \eta <1$, the estimation error is smaller by a factor $T$). This result can be compared to the estimator proposed by gg95, which is studied under finite second moment and requires a full specification of the VECM. Furthermore, ((ref)) can also be compared with the findings in the large factor models literature (see, in particular, Lemma 2 in bai04). As far as uniform rates in $t$ are concerned, in the proof of the theorem we also show the strongest result

This result is unlikely to be sharp, since it is based on the (rather crude) fact that the maximum of any $T$-dimensional sequence with finite $p$-th moment is bounded by $O_{P}\left( T^{1/p}\right) $.

Extensions

The framework developed in the previous section does not allow for deterministic terms in the data, and requires $\varepsilon _{t}$ to be identically distributed. We now discuss possible extensions of our set-up, to accommodate for heterogeneous innovations and deterministics. The main result of this section is that our procedure can be implemented even in these cases, with no modifications required.

Heterogeneous innovations

We consider a novel framework where we allow for innovation heterogeneity of a very general form. Specifically, we assume that

where $u_{t}$ satisfies Assumption (ref) and $h\left( \cdot \right) $ is a deterministic function. The representation in ((ref)) has also been employed in order to consider the presence of heteroskedasticity in the case of data with finite variance (see e.g. cavaliere2009).

The only requirement on the scale function $h\left( \cdot \right) $ is that it has bounded variation within the interval $\left[ 0,1\right] $. Functions of bounded variations and their properties are very well studied, and we refer to e.g. hewitt for details. The design in ((ref)) includes several cases which can be of potential interest: $h\left( r\right) $ can be piecewise linear, i.e. $h\left( r\right) =\sum_{i=1}^{n+1}h_{i}I\left( c_{i-1}\leq r<c_{i}\right) $, with $c_{0}=0$ and $c_{n+1}=1$, thus considering the possible presence of jumps/regimes in the heterogeneity of $ \varepsilon _{t}$; or it could be a polynomial function.

Repeating the proof of Theorem (ref), the results in Corollary (ref) entail that, using the Algorithm 1 in Section (ref), $P^{\ast }(\widehat{m}=m)\rightarrow 1$ with probability tending to 1: $ \widehat{m}$ is still a consistent estimator of $m$. Similarly, it can be shown that the asymptotic properties of our estimator of the common trends and associated loadings are unaffected by heterogeneity in the scale, as in ( (ref)).

Deterministics

In this section we show that our procedure can be also applied to data which contain a deterministic term. Specifically, we consider the representation

where $C$ and $C^{\ast }\left( L\right) $ are defined as before. Equation ( (ref)) is derived from the multivariate Beveridge-Nelson decomposition of the filter $C\left( L\right) $. As is known, ((ref)) can also be obtained from a VECM representation (see ling and yapreisnel)

under the constraint $\mu =\alpha \rho $ with $\rho$ an $m \times 1$ vector.

In this case, our procedure still yields the same results as without the deterministic term.

Monte Carlo evidence

In this section, we discuss the implementation/specification of our procedure and illustrate its finite sample properties through a small scale Monte Carlo exercise. To save space, we report only a limited number of results; further results are in the Supplement.

Design

Following ling, we simulate the $N$-variate $VAR\left( 1\right) $ model

initialized at $y_{0}=0$ and where $A=I_{N}-PP^{\prime }$, $P$ being an $ N\times \left( N-m\right) $ matrix with orthonormal columns (i.e., $ P^{\prime }P=I_{N-m}$).\footnote{ Specifically, in order to create $P$, we have used $P=D\left( D^{\prime }D\right) ^{-1/2}$, where $\left( M\right) ^{-1/2}$ denotes the Choleski factor of a matrix $M$, and have set $D\sim \mathbf{1}_{N\times (N-m)}+d_{N\times (N-m)}$, where $\mathbf{1}_{N\times (N-m)}$ is an $N\times (N-m)$ matrix of ones and $d_{N\times (N-m)}$ is an $N\times (N-m)$ matrix such that $vec\left( d_{N\times (N-m)}\right) \sim N\left( 0,\mathbf{1} _{N(N-m)}\right) $. In particular, $d_{N\times (N-m)}$ is kept fixed across Monte Carlo iterations.}.The innovations $\varepsilon _{t}$ in ((ref)) are i.i.d. and coordinate-wise independent, from a power law distribution with tail index $\eta \in \left\{ 0.5,1,1.5,2\right\} $. We follow the procedure proposed by clauset and generate $\varepsilon _{i,t}$ as

where $v_{i,t}$ is $i.i.d.U\left[ 0,1\right] $; $\varepsilon _{i,t}$ is subsequently centered.\footnote{ In unreported experiments, we have also considered the case $\varepsilon _{t}\sim i.i.d.N\left( 0,I_{N}\right) $, but results are essentially the same as in the case $\eta=2$.}

First, we note from unreported experiments that our procedure for the determination of the number of common trends is not particularly sensitive to the choice of the various specifications. In our experiments, we have used $M=100$ to speed up the computational time, but we note that results do not change when setting e.g. $M=T$, $M=T/2$ or $M=T/4$. In ((ref)), we have used $\kappa =10^{-4}$. This is a conservative choice, whose rationale follows from the fact that, in ((ref)), dividing by $ T^{\kappa }$ serves the purpose of making the non-spiked eigenvalues drift to zero. The upper bound provided in ((ref)) for such non-spiked eigenvalues is given by slowly varying functions, which suggests that even a very small value of $\kappa $ should suffice. Indeed, altering the value of $ \kappa $ has virtually no consequence. In order to compute the integral in ( (ref)), we use the Gauss-Hermite quadrature.\footnote{ In our case, we have used

where the $z_{s}$s, $1\leq s\leq n_{S}$, are the zeros of the Hermite polynomial $H_{n_{S}}\left( z\right) $ and the weights $w_{s}$ are defined as

Thus, when computing $\theta _{T,M}^{\left( j\right) }\left( u\right) $ in Step 2 of the algorithm, we construct $n_{S}$ of these statistics, each using $u=\pm \sqrt{2}z_{s}$. The values of the roots $z_{s}$, and of the corresponding weights $w_{s}$, are tabulated e.g. in salzer. In our case, we have used $n_{S}=2$, which corresponds to $u=\pm 1$ with equal weight $\frac{1}{2}$; we note that in unreported experiments we tried $ n_{S}=4$ with the corresponding weights, but there were no changes up to the $4$-th decimal in the empirical rejection frequencies.} Finally, as far as the family-wise detection procedure is concerned, the level of the individual tests is $\alpha \left( T\right) =0.05/T$, as also used in bt1; this corresponds to having a critical value $c_{\alpha }$ which grows logarithmically with $T$. All routines are based on $1,000$ iterations and are written using GAUSS 21.

Results

Results are reported in Tables (ref)-(ref), where we analyse the properties of our estimator of $m$ with $N\in \left\{ 3,4,5\right\} $. The reported frequencies of the estimates of $m$ show that the finite sample properties are largely satisfactory. Our procedure seems to be scarcely affected by the value of $m$, although, especially for the smaller sample sizes, it appears to be marginally better when $m=0$ as opposed to the case $m=N$. This difference, however, vanishes as $T$ increases. The impact of $N$ is also very clear: as the $VAR$ dimension increases, the performance of $\widehat{m}$ tends to deteriorate, as expected.\ Inference improves for larger values of $T$. Indeed, whilst results for $N=3$ are good even when $\eta =0.5$ and $T=100$, when $N=5$ the estimator $\widehat{m}$ requires at least $T=200$ in order to have a frequency of correctly picking the true value of $m$ higher than $90\%$. This is, as noted above, more pronounced when $m=N$, and less so when $m=0$. As it can also be expected, our procedure improves as $\eta $ increases; results are anyway very good even in the (very extreme) case $\eta =0.5$, and the impact of $\eta $ is less and less important as $T$ increases. Finally, although Tables (ref)-(ref) focus only on the i.i.d. case, unreported experiments showed that results are essentially the same when allowing for serial dependence.

In the Supplement, we report a broader set of results which, in addition to serial dependence in the errors $\varepsilon _{i,t}$, also compare the proposed method with classic information criteria. Results are in Tables (ref)-(ref). Broadly speaking, our procedure is very good on average at estimating $m$ - and better than the best performing information criterion, BIC -- for all values of $N$ and $T$ (and $\eta $). This is true across all values of $m$, including the stationary case ($m=0$) and the no-cointegration case ($m=N$). Information criteria seem to perform marginally better when $m=0$, but this is more than offset when considering that they tend to overestimate $m$ in general, especially so when $m=N$ and $ m=N-1$. When errors are serially correlated (see Tables (ref)- (ref)), results are affected, albeit marginally, but the relative performance of the various methods remains as described above. Finally, in order to evaluate our procedure for medium-large systems, in Tables (ref)-(ref) in the Supplement we report results from a further experiment where we set $N\in \left\{ 10,15,20\right\} $. We find that our estimator delivers a superior performance with respect to information criteria for $m=N $ and $m=N-1$, especially for small values of $ \eta $ (when BIC\ grossly understates $m$, even for large values of $T$); BIC is (albeit marginally) superior when $m=0$.

Real data examples

We illustrate our methodology through four empirical applications to: commodity prices (Section (ref)),\ U.S. interest rate data (Section (ref)), long run purchasing power parity (Section (ref)) and cryptocurrency markets (Section (ref)). In these applications, $N$ ranges from $3$ to $7$.

Comovements among commodity prices

We consider a set of $N=7$ commodity prices: three oil prices (WTI, Brent crude, and Dubai crude) and the prices of four metals (copper, gold, nickel, and cobalt). The presence of common trends can be anticipated due to global demand factors (e.g. growth in emerging Asian countries and especially in China; or changes of preferences towards greener energy sources, which increase demand for copper and decrease demand for oil), and also due to global supply factors (e.g. related to the effect that oil prices have on transportation costs of other commodities; or driven by technological innovations which often require the use of cobalt -- see e.g. alquist). Moreover, the three oil prices should exhibit strong comovements, and similarly should the prices of metals, which are often used in combination in industry (e.g. copper and nickel). In order to study the presence of such common trends, we use a dataset consisting of monthly data from January 1990 to March 2021, corresponding to a sample of $T=373$ monthly observations.\footnote{ Data have been downloaded from https://www.imf.org/en/Research/commodity-prices} We use the logs of the data, which are subsequently demeaned and detrended.

We have applied our methodology using the same specifications as described in Section (ref), i.e. $\kappa =10^{-4}$, $M=100$ and $n_{S}=2$ in ((ref)). In order to assess robustness to these specifications, we have also considered other values of $M$ (including $M=T$ ) and $n_{S}=4$.

We report the results for the $7$-dimensional series in Table (ref). Initially, we report the (Hill's) estimates of the tail indices for the seven series; the associated confidence sets are quite large, but the test by trapani16 supports the hypothesis that all series have infinite variance.

Estimation of $m$ based on Johansen's sequential procedure for the determination of the cointegration rank (using either the trace tests or the maximum eigenvalue tests) provides ambiguous results, with the estimate of $ m $ ranging between $5$ and 7 (which corresponds to no cointegration). These results, in the presence of heavy tails, are not reliable: the simulations in caner show that Johansen's procedure does not control size under infinite variance, thus making the estimate of $m$ inaccurate (in particular, it tends to be over-sized, thus leading to overestimation of $m$ ). In contrast, through our test we find strong evidence of $m=4$ common stochastic trends. As shown in the table, our results are broadly robust to different values of $M$, $\kappa $. In (much) fewer cases, we find $m=5$, which might suggest the presence of a slowly mean reverting component in the data.

In order to shed more light on these findings, we split the series into two sub-groups: one of dimension $N=3$ (comprising the three crude prices -- Brent, Dubai and WTI crude), and one of dimension $N=4$ (containing the four metal prices). Results for the $3$-dimensional series of crude prices are in Table (ref). On the one hand, Johansen's tests in this case identifies (at $5\%$ level, and only using the trace test) two common trends ($m=2$). On the other hand, our methodology provides evidence of a single ($ m=1$) common stochastic trend (and, in some, more rare, cases, of $m=2$). Results concerning the $N=4$ metals are in Table (ref); in this case, evidence of $m=3$ common trends emerges from all the procedures considered.

Overall, most of the evidence points towards $m=4$ common stochastic trends; there is much less evidence in support of $m=5$. The $m=4$ common trends can be estimated as explained in Section (ref). In order to identify the trends, based on the results above we propose to order the series as follows: WTI, gold, cobalt, copper, Brent crude, Dubai crude, nickel. Then, we constrain the upper $m\times m$ block of the estimated loadings matrix $ \widehat{P}$ to be the identity matrix. In this way the trends are identified with the first four series, i.e. WTI, gold, cobalt, and copper. In Figure (ref), we plot the estimated common trends $\widehat{x} _{t}$ and in Table (ref) we report the associated loadings $ \widehat{P}$.

By construction, the first and second trends ($\widehat x_{1,t}$ and $ \widehat x_{2,t}$) are associated with oil prices and cobalt repsectively. The third one ($\widehat x_{3,t}$) is associated with gold (by construction), and nickel, with a negative loading; finally, the fourth trend ($\widehat x_{4,t}$) is associate with copper by construction, and with nickel (with a positive loading). The trends driving metals are also common to oil prices, albeit with smaller loadings.

The term structure of US interest rates

In this example, we illustrate our methodology through a small scale application based on the same dataset as in ling, consisting of $N=3 $ monthly, seasonally unadjusted series: the 3-month Treasury Bill rate, the 1-year Treasury Bill rate (both from sales on the secondary market) and the Effective Federal Fund rate, spanning the period between February 1974 and February 1999 (thus corresponding to $T=301$ monthly observations).\footnote{ The data have been downloaded from the Federal Reserve Economic Data website, https://fred.stlouisfed.org}

As in ling, we consider the logs of the original data, say $ Y_{t}=\left( y_{1,t},y_{2,t},y_{3,t}\right) ^{\prime }$ where $y_{1,t}$ is the log of the 3-month Treasury Bill rate, $y_{2,t}$ is the log of the 1-year Treasury Bill rate, and $y_{3,t}$ represents the log of the Fed Fund rate. ling carry out their analysis using a VAR(2)\ in error correction format:

with $\mu =\alpha \rho$, such that the data do not have a linear trend but may have a non-zero mean in the stationary directions. In the light of Corollary (ref), in this case our set-up can be applied with no modifications. Moreover, we do not need to make specific assumptions on the lag structure for $\Delta Y_{t}$.

As demonstrated above, the main feature of our contribution is that we do not require any prior knowledge of the index $\eta $. This advantage can be further understood upon observing the estimates of $\eta $ obtained in \citet *{ling} -- see in particular, their Figure 6, where $\eta $ ranges between $ 1 $ and $1.5$ (thereby implying that the data have infinite variance and, possibly, infinite mean). On account of this information, ling report the critical values for the Likelihood Ratio tests in caner for $\eta =1 $ and $\eta =1.5$, using both values to build a decision rule.

In particular, ling conclude that there is some evidence in favour of $m=1$ (corresponding to two cointegration relationships). Upon inspecting their results, however, it seems that they are driven by the choice of the $ 5\%$ nominal level: evidence in favour of $m=1$ is much weaker when considering a $1\%$ nominal level (in particular, in the latter case, should one use the critical values computed for $\eta =1$, there would be no evidence of cointegration at all, whereas using the critical values computed for the case $\eta =1.5$, one would conclude that $m=1$). Interestingly, assuming finite variance, Johansen's sequential procedure rejects $m=2$ and $ m=3$ in favour of $m=1$, even at the $1\%$ nominal level.

In contrast to these mixed results, we find $\widehat{m}=2$ common stochastic trends, corresponding to one cointegration relationship. This result is robust to different values of $M$ and $\kappa $.

Testing for the weak form of the PPP

In this application, we check the validity of the weak form of the PPP, using the same dataset as in falk. In particular, for a set of 12 countries, we individually verify the presence of cointegration in the vectors $\left( p_{c,t},p_{US,t},FX_{c,t}\right) ^{\prime }$, where $p_{c,t}$ is the (log of the) CPI\ of country $c$ at time $t$, $p_{US,t}$ is the (log of the) US CPI index, and $FX_{c,t}$ is the log of the exchange rate of country $c$ currency vis-a-vis the dollar. The series are monthly, spanning from January 1973 to December 1999, corresponding to $T=324$ for each country.

A minimum requirement for the weak form of the PPP to hold is that there is at least one cointegrating relation between the three series -- i.e., at most $m=2$ common stochastic trends. falk find evidence of heavy tails in this dataset (see their Table III), with series typically having finite mean but infinite variance. Thence, the authors test for cointegration using critical values based on allowing for heavy tails (notice that, even in this case, the authors require prior knowledge of the tail index $\eta $). We compare results from falk and Johansen's procedure with our proposed tests. As in the previous sections we use $M=T$, $\kappa =10^{-2}$; the level for the individual tests is set to $0.05/\ln T$ . Tests are based on the deviations from the initial values. Results are reported in Table (ref).

First, we notice that Johansen's method at the $10\%$ ($5\%$) nominal level under the assumptions of finite variance indicates that the necessary condition for the weak form of the PPP holds for 10 (8) out of 12 countries. Second, if the critical values are adjusted using the estimated tail index, as in falk, then the support for the weak from of the PPP decreases. Specifically, they found that the number of countries with $\widehat{m}\leq 2 $ ranges between $3$ (when using a maximum eigenvalue test at the $5\% $ nominal level) and $8$ (when using a trace test at the $10\%$ nominal level). Our test provides support of PPP in $5$ countries out of 12, which is closely in line with falk. Our empirical evidence reinforces the view that the assumption of finite variance can lead to spurious support to the PPP.

Statistical arbitrage in cryptocurrency markets

In the context of asset pricing, since the seminal contribution by diebold1994cointegration, common trends and cointegration have played a prominent role. Indeed, the presence of cointegration entails the possibility of constructing portfolios which are mean reverting, so that a profit can be made when such portfolios depart from their mean -- the so-called \textquotedblleft statistical arbitrage\textquotedblright . Although the literature usually focuses on portfolios of two assets (\textquotedblleft pairs trading\textquotedblright ), it is possible to generalise the notion of statistical arbitrage to the multi-asset case (see alexander). In such a case, having $m$ common stochastic trends allows to construct $N-m$ mean reverting portfolios, each having weights given by a user-chosen linear combination of the cointegrating vectors. Cryptocurrency data have been paid increasing attention (makarov); however, despite the evidence that returns on cryptocurrencies exhibit heavy tails (see pele), this is not usually accounted for in applications.

We consider a portfolio of $N=6$ cryptocurrencies, chosen among those with the largest market capitalisation within our sample period: Cardano (ADA), Bitcoin Cash (BCH), BitCoin (BTC), EOS, Ethereum (ETH), and XRP.\footnote{ Data have been downloaded from https://www.coingecko.com/} We use daily data from October 18, 2017 until September 18, 2020, which is equivalent to a sample of $T=1,066$ observations. As in the previous applications, we consider deviations from the initial value. Table (ref) shows that all cryptocurrencies have heavy tails, with the estimated tail indices around $1$; despite some apparent heterogeneity in the tail indices, confidence intervals for the Hill's estimator show that the common tail index assumption which characterises stable distributions is satisfied.

The results in Table (ref), based on our test using different configurations of the tuning parameters $M$ and $\kappa $, indicate $m=5$ common trends for all cases considered, i.e. the presence of one cointegrating relation between the $6$ cryptocurrencies considered. Note that, in contrast to our findings, Johansen's approach based on a VAR($p$) with $p=10$ (as suggested by BIC) broadly indicates $m=4$ ($m=3$ at the $ 10\% $ nominal level or for larger values of $p$). However, the results by caner indicate a tendency to understate $m$ in the presence of heavy tails; therefore, the additional cointegrating relations found by Johansen's method is likely to be spuriously driven by the heavy-tailedness of the cryptocurrency returns.

Conclusions

In this paper, we propose a methodology for inference on the common trends in multivariate time series with heavy tailed, heterogeneous innovations. We develop: (i) tests for hypotheses on the number of common trends; (ii) a sequential procedure to consistently detect the number of common trends; (iii) an estimator of the common trends and of the associated loadings. A key feature of our approach is that estimation of the tail index of the innovations is not needed, and no prior knowledge as to whether the data have finite variance or not is required. Indeed, the procedure can be applied even in the case of finite second moments. \newline Our method is based on the eigenvalues of the sample second moment matrix of the data in levels, the largest $m$ ($m$ being the unknown number of common trends) of which are shown to diverge at a higher rate, as $T$ increases, than the remaining ones. Based on such rates, we propose a randomised statistic for testing hypotheses on $m$; its limiting distribution is chi-squared under the null, and diverges under the alternative. Combining these individual tests into a sequence of tests, we prove consistency of the estimator of $m$ by simply letting the nominal level to shrink to zero at a proper rate. We also show that, once $m$ is determined, estimation of the common trends and their loadings can be done using PCA. Our simulations show that our method has good properties, even in samples of small and moderate size.

{

}

\setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{theorem}{0}

Further Monte Carlo evidence

Top-down estimation of $m$

As mentioned in Section (ref) in the main paper, in order to determine the number of common trends $m$, it is also possible to consider a top-down algorithm as an alternative to the \textquotedblleft bottom-up\textquotedblright\ Algorithm 1, i.e. Algorithm 2. Its asymptotic validity is demonstrated next.

In Tables (ref)-(ref), we report the frequencies of the estimates of $m$, obtained from the same design as in Section (ref) in the main paper. As can be seen, results are as good as with $ \widehat{m}$.

Further evidence on $\widehat{m}$

Here, we use the same DGP as in Section (ref), viz.

with $A=I_{N}-PP^{\prime }$ and $P$ an orthonormal $N\times \left( N-m\right) $ matrix; however, as well as considering i.i.d. errors, we also allow for serial correlation in $\varepsilon _{t}$. Specifically, the innovations $\varepsilon _{t}$ in ((ref)) have again been generated as coordinate-wise independent, and we consider serial dependence through an $AR\left( 1\right) $ specification

for $1\leq i\leq N$; we set $\theta \in \left\{ 0,0.5,-0.5\right\} $. The innovation $e_{i,t}$ is generated as i.i.d. across $t$, with a power law distribution of tail index $\eta $, using $\eta \in \left\{ 0.5,1,1.5,2\right\} $ and again, as a benchmark, we consider the case $ \varepsilon _{t}\sim i.i.d.N\left( 0,I_{N}\right) $. In order to simulate the power law distribution, we generate

for $1\leq i\leq N$, where $v_{i,t}$ is $i.i.d.U\left[ 0,1\right] $; $ e_{i,t} $ is subsequently centered. The first $1,000$ datapoints are discarded in order to avoid dependence on initial conditions.

All other specifications are the same as in the main paper.

We use information criteria as a comparison. In particular, we consider the three most popular ones (i.e., AIC, BIC and HQ), but in the tables we only report the results corresponding to BIC, which is consistently the best out of the three for all experiments\footnote{ We note that AIC and HQ tend to grossly overstate $R$ in all experiments, with this tendency worsening for lower values of $\eta $. Results are available upon request.}. We point out that the use of information criteria in the context of cointegration analysis is well-studied in the finite variance case (see e.g. \citealp*{aznar} and \citealp*{ cavaliere2018determining}); conversely, the case of infinite variance has not been studied in the case of cointegrated systems (see however the contribution by \citealp*{knight}). In our experiments, we use an \textquotedblleft infeasible\textquotedblright\ version of each information criterion, assuming that the lag structure in ((ref)) is known, so that the only quantity to be determined is the number of common trends $m$.

For all rank estimates, we report three measures of performance. Denoting the estimate of $m$ at iteration $1\leq mc\leq 1,000$ as $\widehat{m}_{mc}$, we define

Results are reported in Tables (ref) and (ref), where we analyse the properties of our estimator of $m$ for a small-scale $ VAR$ with $N\in \left\{ 3,4,5\right\} $, and for the simple case of i.i.d. errors. Broadly speaking, as far as the $PCW$ and $ME$ indicators are concerned, the tables show that our procedure is very good on average at estimating $m$ - and better than the best information criterion, BIC - for all values of $N$ and $T$, even for the (rather extreme) case $\eta =0.5$. This is true across all values of $m$, including the case of stationarity ($ m=0$) and lack of cointegration ($m=N$). Such evidence is further corroborated by considering the $STD$ indicator, which shows that our procedure is systematically better than BIC in this respect. Indeed, the BIC performs (very) marginally better when $m=0$, but this is more than offset when considering that BIC tends to overstate $m$ in general (as the values of $ME$ indicate), and especially when $m=N$ and $m=N-1$; in particular, BIC seems to perform puzzlingly badly when errors are Gaussian. This is even more remarkable as information criteria are based on the (infeasible) assumption that the lag structure is known. In line with the theory, results are better for larger $T$; also, results do improve as $\eta $ increases, although this becomes less evident as $T$ increases. The impact of $N$ is less clear, although generally speaking larger values of $N$ worsen the performance of both BIC and our procedure when $m=N$, whereas it has virtually no impact for other values of $m$. \newline

In Tables (ref)-(ref) we consider serially correlated innovations. When $\theta =0.5$ in ((ref)), our procedure fares better than in the i.i.d. case when $m=N$, worsening instead when $m=0$. In the latter case, as $T$ increases (and as $ \eta $ becomes closer to $2$), the impact of serial correlation is mitigated -\ see Tables (ref) and (ref). Conversely, the case of negative autocorrelation (Tables (ref) and (ref) ), results worsen in a more pronounced way, especially for small values of $ T $ and/or as $N$ increases. Even in this case, however, our procedure is better than BIC for the case $m=N$ (and in general for large $m$), with an almost indistinguishable performance when $m=0$.

Finally, and as also mentioned in the main paper, in Tables (ref)- (ref), we report results using $N\in \left\{ 10,15,20\right\} $. As overview of the results is also in Section (ref) in the main paper.

\setcounter{equation}{0} \setcounter{lemma}{0}

Technical lemmas

Henceforth, we denote the $L_{p}$-norm of a scalar random variable $X$ as $ \left\vert X\right\vert _{p}=\left( E\left\vert X\right\vert ^{p}\right) ^{1/p}$; also, the subscript \textquotedblleft $_{\perp }$ \textquotedblright\ denotes the orthogonal complement of a matrix. Further, we often use the notation $R=N-m$.

\setcounter{equation}{0} \setcounter{lemma}{0}

Proofs

Henceforth, we let $E^{\ast }$ and $V^{\ast }$\ denote the expected value and the variance with respect to $P^{\ast }$ respectively.

\textwidth 20.5cm \oddsidemargin 2cm \evensidemargin 2cm \textheight 20.3cm \topmargin 1cm