Rank based cointegration testing for dynamic panels with fixed T

In this paper, we show that the cointegration testing procedure of Binder et al. (Econom Theory 21:795–837, 2005) for Panel Vector Autoregressive model of order 1, PVAR(1) is not valid due to the singularity of the hessian matrix. As an alternative we propose a method of moments based procedure using the rank test of Kleibergen and Paap (J Econom 133:97–126, 2006) for a fixed number of time series observations. The test is shown to be applicable in situations with time-series heteroscedasticity and unbalanced data. The novelty of our approach is that in the construction of the test we exploit the “weakness” of the Anderson and Hsiao (J Econom 18:47–82, 1982) moment conditions. The finite-sample performance of the proposed test statistic is investigated using simulated data. The results indicate that for most scenarios the method has good statistical properties. The proposed test provides little statistical evidence of cointegration in the employment data of Alonso-Borrego and Arellano (J Bus Econ Stat 17:36–49, 1999).


Introduction
In this paper, we consider the cointegration testing problem for Panel VAR model of order 1 with a fixed time dimension. Up to date the only testing approach in this Financial support from the NWO MaGW grant "Likelihood-based inference in dynamic panel data models with endogenous covariates" is gratefully acknowledged. case is the likelihood ratio test based on the transformed maximum likelihood (TML) estimator of Binder et al. (2005) [hereafter BHP]. However, in the univariate setup it is known that for data with autoregressive parameter close to unity, the likelihood approach does not have a gaussian asymptotic limit, see e.g. Kruiniger (2013). We extend that result to multivariate setting and argue that the cointegration testing procedure of Binder et al. (2005) is not valid due to the singularity of the corresponding expected hessian matrix.
To the best of our knowledge, in the fixed T dynamic panel data (DPD) literature no feasible method of moments (or least-squares) alternative to likelihood based cointegration testing procedures is available. The main reason for the absence of method of moments based alternatives is that jacobian matrix of the Anderson and Hsiao (1982) moment conditions is of reduced rank, when process is cointegrated. It is natural to use this information and consider a rank based cointegration test, based on the rank of the jacobian matrix. In this paper, we propose such a test and show that it is applicable in situations with time-series heteroscedasticity and, unlike the likelihood based tests, the new test does not require any numerical optimization. At the same time, this procedure cannot provide inference that is uniform over the parameter space, as the asymptotic distribution of the test depends on the properties of the initial condition.
In the Monte Carlo section of this paper, we investigate the finite sample properties of the proposed procedure. We find that the new testing procedure provides a good size control as well as high power in most of the designs considered. However, in some setups this test lacks power if the data generating process for the initial condition substantially deviates from stationarity.
The paper is structured as follows. In Sect. 2, we briefly present the model, the testing problem at hand and the results for the testing procedure of Binder et al. (2005). Rank-based cointegration testing procedure is formally introduced in Sect. 3. In Sect. 4, we continue with the finite sample performance by means of a Monte Carlo analysis. In Sect. 5, we illustrate the testing procedure using the data of Alonso-Borrego and Arellano (1999). Section 6 concludes.
Here we briefly discuss notation. Bold upper-case letters are used to denote the original parameters, i.e. {Φ, Σ, Ψ }, while the lower-case letters {φ, σ , ψ} denote vec (·) (vech (·) for symmetric matrices) of corresponding parameters, in the univariate setup corresponding parameters are denoted by {φ, σ 2 , ψ 2 }. We use ρ( A) to denote the spectral radius 1 of a matrix A ∈ R n×n . We defineȳ i− ≡ (1/T ) T t=1 y i,t−1 and similarlyȳ i ≡ (1/T ) T t=1 y i,t . We usex to indicate variables after Within Group transformation (for exampleỹ i,t ≡ y i,t −ȳ i ), whileẍ are used for variables after a "quasi-averaging" transformation. 2 For further details, see Abadir and Magnus (2002). Where necessary, we use the 0 subscript to denote the true value of the parameters, e.g. Φ 0 .

Cointegration testing for fixed T panels 2.1 The model
In this paper, we consider the following PVAR(1) specification where y i,t is an [m × 1] vector, Φ is an [m × m] matrix of parameters that are of main interest, η i is an [m × 1] vector of unobserved individual specific covariates, and ε i,t is an [m × 1] vector of innovations independent across i, with zero mean and some finite covariance matrix. If one sets m = 1, the model reduces to the univariate linear dynamic panel data model with autoregressive dynamics. Throughout this paper we assume that η i satisfy the so-called "common dynamics" assumption where Π ≡ Φ − I m . If at least one eigenvalue of Φ is equal to unity this assumption ensures that there is no discontinuity in the data generating process (DGP), for further discussion, see e.g. BHP.
Assuming common dynamics we can rewrite the model in (1) as y i,t = Πu i,t−1 + ε i,t , i = 1, . . . , N , t = 1, . . . , T. (2) Here we define u i,t−1 ≡ y i,t−1 − μ i . We say that series y i,t are cointegrated if the Π matrix is of reduced rank 0 < r < m. 3 In particular, there exist full column rank (of rank r ) matrices α r and β r such that: 4 where r is the rank of Π. Matrices α r and β r are not unique, as for any [r ×r ] invertible matrix U This is the so-called "rotation problem". As a result, it is a usual practice in the literature to introduce identifying restrictions on α r or β r . To construct the test statistic that we formally introduce in Sect. 3.1, we follow common practise and assume that β r = (I r , Δ), where Δ is an [r × m − r ] matrix.
As a natural starting point one can consider the fixed effects (FE) estimator for Φ parameterΦ It is well known that the estimator of this type is inconsistent for fixed T , see e.g. Nickell (1981) and Hahn and Kuersteiner (2002). For the unit root case, i.e. Φ 0 = I m (correspondingly r = 0), one can show that suggesting that for T > 2 we can obtain a consistent estimate of Φ 0 = I m , by considering the bias corrected estimator of the form The bias corrected estimatorΦ BC can be used for unit root testing, see e.g. Harris and Tzavalis (1999) and Kruiniger and Tzavalis (2002) for the univariate case. However, no comparable result is available when series are cointegrated, i.e. 0 < r < m. Alternatively, one can rely on the bias-corrected approaches of Bun and Carree (2005), Ramalho (2005) or Dhaene and Jochmans (2016) to obtain a consistent estimator of Φ. However, these procedures generally fail to guarantee correct asymptotic inference for reasons similar to those discussed in the next section. 5

The likelihood ratio test of BHP
In this section we introduce likelihood based testing and estimation approach advocated by BHP. First, we list the assumptions (that we abbreviate as, TML) used to derive asymptotic distribution of the transformed maximum likelihood estimator: (TML 1) The error terms ε i,t are i.i.d. across i and uncorrelated over time (TML 3) The following moment restrictions are satisfied: Assumption (TML 1) is the no serial correlation assumption. Note how Assumption (TML 2) places moment restrictions only on one linear combination of y i,0 and μ i , rather than separately on y i,0 and μ i . Assumption (TML 3) imposes zero covariance restriction on initial deviation and the error terms, which is a standard assumption in the literature. For Π = O, E[u i,0 ε i,t ] is unrestricted, as y i,1 is not a function of u i,0 . Assumption (TML 5) is the main exception as compared to the setup in e.g. Juodis (2016), as in this paper we allow the maximum eigenvalue of the autoregressive matrix Φ 0 to be 1. The exact components and dimension k of the κ vector are related to a particular parametrization of the parameter space used for estimation. Assumptions TML are almost identical to the corresponding assumption FE in Binder et al. (2005), the only difference is that some of their assumptions are "high-level" (e.g. bounds on covariances), while we consider "low-level" assumptions by imposing restrictions directly on ε i,t and u i,0 .
The quasi log-likelihood function for Y i = vec ( y i,1 , . . . , y i,T ) is then defined as follows (up to a constant): where κ = (φ , σ , ψ ) , hence k = m 2 + m(m + 1) and Ψ is the variance-covariance matrix of the initial observation y i,1 . The Σ τ matrix has a block tri-diagonal structure, with −Σ on first lower and upper off-diagonal blocks, and 2Σ on all but the first (1, 1) diagonal blocks. The first (1,1) block is set to Ψ which takes into account the fact that we do not restrict y i,1 to be covariance stationary. 6 The [mT × mT ] R matrix has I m matrices on the diagonal blocks, and −Φ on the first lower off-diagonal blocks.
In Juodis (2016) it is shown that the log-likelihood function of BHP can be substantially simplified to If the matrix Φ− I m = α r β r is of a reduced rank r (cointegration), 7 this information can be taken into account in estimation and used for testing. To avoid rotational indeterminacy, one can use the same parametrization as BHP and set β r = δ r H r + b r , where H r and b r are known and δ r is an [m −r ×r ] (given 0 < r < m) matrix of parameters. The parameter set in this case is defined as κ r = ((vec α r ) , (vec δ r ) , σ , θ ) . Binder et al. (2005) suggest to use the likelihood function (8) in constructing the likelihood ratio test statistic to consistently estimate rank r . In particular, the likelihood ratio test statistics for the null hypothesis H 0 : r 0 = r against the alternative H A : r A = r + 1 is of the form L R(r, r + 1) = −2( (κ r ) − (κ r +1 )).
Similarly to the classical maximum eigenvalue test of Johansen (1991) for N = 1, the overall testing procedure can be performed sequentially, i.e. by first considering r 0 = 0 (Π 0 = O m ), and if rejected, proceeding with r 0 = 1, and so on. BHP argue that this procedure under H 0 : r 0 = r has a χ 2 (·) asymptotic limit. In particular, in Remark 4.1. they note that: "Unlike in time-series models, first differencing in panels with T fixed still allows identification and estimation of the long-run (level) relations that are of economic interest, irrespective of the unit root and cointegrating properties of the y i,t process.". As we show next, this conclusion is not completely correct. One of the standard regularity conditions for extremum estimators, is that the asymptotic (or expected) hessian matrix, H ≡ E[H N (κ 0 )] is positive definite. In Bond et al. (2005), authors showed that for the TML estimator of Hsiao et al. (2002) (which is a special case of Binder et al. 2005 for m = 1) this regularity condition is violated. In the next theorem we show that the same conclusion extends to a more general case with m ≥ 1.

Theorem 1 (Singularity) Let Assumptions TML be satisfied. Then at
Proof In the Appendix.
As the TML estimator can be seen as a non-linear MM estimator with the score vector defining the moment conditions, singularity of the H matrix can be seen as a "weak instrument" problem (using the GMM notation). The singularity result in Theorem 1 is of special interest when the inference regarding the rank of I m − Φ 0 is concerned. It is important to note that despite singularity of H , the TML estimatorκ T M L E remains consistent, hence the identification part of Remark 4.1. in BHP is correct. However, as a result of singularity the limiting distribution for this estimator is nonstandard. Using the approach of Roznitzky et al. (2000), Ahn and Thomas (2006) showed that in the univariate model (i.e. m = 1), the TML estimator of φ converges at the N 1/4 rate to a non-standard distribution. 8 Additionally, they show that LR test statistic for H 0 : φ 0 = 1 has a mixture distribution, of a χ 2 (1) random variable and a degenerate random variable that takes value 0 with probability 1, with equal mixing weights of 0.5. In this paper, we do not attempt to study the distributional consequences of the singularity for the LR test and leave it for future research. 9 Based on results in Dovonon and Renault (2009) (for GMM), it is known that for general rank deficiencies the maximal rate of convergence is N 1/4 . However, no results regarding the behavior of the estimator (see discussion in Dovonon and Hall 2016) and the LR ratio test in cases like ours are available. As a result, it is not obvious that using the critical values from the χ 2 distribution with (m − r ) 2 degrees of freedom results in a conservative test.
Although the unit root model is not of prime importance for the main topic of this paper, Theorem 1 provides a natural starting point for intuition of the next result. For the unit root case (i.e. Φ 0 = I m ) the expression for H simplifies dramatically as Σ 0 = Θ 0 . That allowed us in Theorem 1 to show that |H | = 0 for any value of Σ 0 and T . Unfortunately, no result of this type is available when Π is of reduced rank r > 0. However, some special results can be derived for T = 2.
Proposition 1 Let Φ 0 be such that rk Π 0 = r and T = 2 then Proof In the Appendix.
This quantity is smaller than m 2 for all m ≤ 4 (note that the bivariate PVAR model is analyzed in most empirical studies with limited number of time-series observations). It follows that for cases of most empirical value the expected hessian matrix is singular and the corresponding estimator does not have a normal limiting distribution. Although in this paper we do not prove more general results for T > 2, we performed numerous numerical evaluations of H for larger values of T and different combinations of population matrices in the bivariate setup. 10 For all setups we found that the expected hessian matrix is singular for r < m and of full rank otherwise. Given these results the unit root and cointegration testing procedure of BHP that is based on asymptotic χ 2 (·) critical values is not asymptotically valid.
Remark 1 Alternatively, instead of considering likelihood function for observations in first differences one can consider a correlated random effects likelihood function (conditional on y i,0 ) as in Arellano (2016) and Kruiniger (2013). Although we do not formally consider a possible singularity of the hessian matrix for that estimator, we conjecture that the main conclusions of this paper are also applicable to that approach (Ahn and Thomas 2006;Kruiniger 2013 proved this for m = 1).
Remark 2 Note that the results of this section are derived under assumption that Ψ is estimated without any restrictions, i.e. as suggested by Binder et al. (2005). If one instead imposes some restrictions on this parameter matrix, e.g. covariance stationary, it is possible that the expected hessian matrix has full rank. For example, Kruiniger (2008) considers univariate results, where he shows that for φ 0 = 1, the TML estimator retains standard asymptotic properties if the stationarity assumption is used in estimation.

Regularity conditions
In this section we propose an alternative approach to cointegration testing. To explain the intuition of our approach, consider the following Anderson and Hsiao (1982) moment conditions for Panel VAR(1) model (see e.g. BHP) where only the most recent lag, y i,t−2 is used as an instrument (other choices are discussed later in the paper). The (minus) jacobian of these moment conditions is given by From the properties of the Kronecker product it follows that the rank of this matrix is determined by the rank of the matrix inside the brackets. 11 That term can be expanded as follows (upon redefining t → t + 1, as the previous expression is well defined for Under the no serial correlation assumption, e.g. (TML 1), E[ε i,t y i,t−1 ] = O m , while the first term is the product of rank r and rank r uy ≤ m matrices. As a result rk (E[ y i,t y i,t−1 ]) = min (r, r uy ) < m and this leads to a violation of the "relevance" condition for the Instrumental Variable (IV) estimator, thus the Anderson and Hsiao (1982) moment conditions cannot be used to consistently estimate Φ. 12 However, we can use the jacobian matrix directly to test for cointegration, avoiding the estimation step. Next, we list assumptions that are sufficient to derive the asymptotic properties of the testing approach that we introduce in this section. For the purpose of this section we deviate from (TML) assumptions and restrict moments of μ i and y i,0 separately, rather than their linear combination.
(A.1) The error terms ε i,t are i.i.d. across i and uncorrelated over time, Note that we allow ε i,t to be heteroscedastic over time. However, cross-sectional heteroscedasticity is in general not allowed, as in this case E[ y i,t y i,t−1 ] is individual specific, and we cannot consistently estimate both the mean and the variance of y i,t y i,t−1 . In particular, for this reason we assume that μ i are iid. As we consider fixed T panels, it is important to explicitly specify DGP for the initial conditions. For this, we assume that y i,0 are of the following form Here Υ is an [m × m] matrix that controls the degree of non-stationarity in the initial condition, i.e. if Υ = I m , initial condition is effect non-stationarity. 13 Note that if Υ = I m the assumption (A.2) can be somewhat relaxed allowing heterogenous μ i . What is left is to specify assumptions on ε i,0 . Below we list a few DGPs for ε i,0 that can be used for our purpose.
Here M is assumed to be finite.
Here ξ i is an [m × 1] vector of the (independent) individual-specific initialization effects. 14 In what follows, we assume that all random variables in (DGP.1)-(DGP.3) satisfy assumptions (A.1)-(A.2). Furthermore, in (DGP.3) for simplicity all ε * i,−l are homoscedastic over the time-series dimension. This restriction can be relaxed by assuming that all Σ * −l are appropriately summable as l → ∞. (DGP.3) initialization was used in the Monte Carlo studies of BHP and is motivated by the Granger Representation Theorem, see e.g. Theorem 4.2 in Johansen (1995). The (DGP.2), was e.g. used by Hayakawa (2016). It is important to emphasize that all three DGP are well defined for all values of r . 15

Rank test
In this paper, we use the generalized rank test of Kleibergen and Paap (2006) as a basis for our testing procedure. Here we briefly introduce their testing procedure and later apply it to our problem. In construction of the rank test Kleibergen and Paap (2006) use the property that any [k × f ] matrix D can be decomposed as: ] matrix and all ⊥ matrices are defined in the usual way. For Λ q = O the rank of D is determined by the rank of A q B q . The procedure in Kleibergen and Paap (2006) with matrices A q , B q , Λ q obtained using the singular value decomposition (SVD). In our case, we consider singular value decomposition of jacobian matrix, i.e. D = E[ y i,t y i,t−1 ].
13 Also referred as "mean non-stationarity", see e.g. Bun and Sarafidis (2015). Next, we define the sample analogue of D: Applying the standard Lindeberg-Levý CLT, it follows that: Here the full rank matrix V can be consistently estimated using its finite sample counterpart: Consequently, the estimator y i,t y i,t−1 satisfies sufficient conditions in Kleibergen and Paap (2006). 16 As a result one can apply Theorem 1 of Kleibergen and Paap (2006) to the problem at hand: Matrices A and B in Theorem 2 are obtained from the SVD of y i,t y i,t−1 . An operational version of the rk(r ) test statistic is obtained by replacing the (unknown) matrix Ω r with some consistent estimator. An obvious choice forΩ r is given by: 16 y i,t y i,t−1 satisfies Assumption 1 (asymptotic normal distribution), while V satisfies Assumption 2(full rank). Note that if V is of reduced rank the result in Theorem 2 can be modified as in Andrews (1987).
Note that this test statistic uses the unrestricted estimate of E[ y i,t y i,t−1 ], hence we do not explicitly specify the alternative hypothesis (as it is done for LR test of e.g. Johansen 1991or Binder et al. 2005. However, as we discuss next, this testing approach has power only towards alternative with r 0 > r , i.e. test rejects if the true rank is larger than the hypothesized one. The result in Theorem 2 suggests that one can use a sequential testing procedure in order to determine the rank of E[ y i,t y i,t−1 ]. In particular, one can begin with H 0 : rk E[ y i,t y i,t−1 ] = 0, and if rejected, consider H 0 : rk E[ y i,t y i,t−1 ] = 1, and so on, until the first non-rejection. The construction of such sequential procedure does not differ from the one suggested in Johansen (1991) or Binder et al. (2005). However, E[ y i,t y i,t−1 ] is not of the prime interest for us, as we are interested in testing the rank of Π. This begs the question:

Under which conditions one can interpret rejection/non-rejection of the rk(r ) test as an evidence regarding the rank of Π?
If one rejects the null hypothesis has a reduced rank if and only if y i,t are cointegrated (the "if" part was established above). Hence, it still remains to be investigated under which conditions the E[ y i,t y i,t−1 ] term is of reduced rank if and only if Π is of reduced rank. Let us investigate this issue more closely by expanding the In the effect-stationary setup (Υ = I m ) all terms involving Υ are equal to O m . Furthermore, the third term is a p.d. matrix as all Σ s matrices are positive definite. Moreover, is also at least positive semi-definite (p.s.d.) matrix. Thus we conclude that the "only if" part is also valid for Υ = I m .
Unfortunately, there is a lot of evidence in the DPD literature suggesting that in general this assumption can be too restrictive, see e.g. Arellano (2003) and Roodman (2009). If Υ = I m , the first, third and fourth terms are p.s.d. matrices, while it is not immediately clear what happens with the second term. For the procedure to consistently estimate the rank of Π (and not only underestimate it), one has to place additional restriction Note that restriction is "high-level" as it imposes restrictions not directly on parameters, but instead on their non-linear function. Intuitively, assumption (IDN) is satisfied if Φ t−1 (Υ − I m )Σ μ is "large", with eigenvalues sufficiently bounded away from zero. However, it is not a trivial task to identify the parameter space of {Φ, Υ , Σ μ } for the aforementioned condition to be satisfied. One special case is obtained for Υ = I m (effect stationarity) with other matrices being unrestricted (at least finite). If we can ensure that has reduced rank r if and only if y i,t−1 are cointegrated. 17 In the Monte Carlo section of this paper, we check the adequacy of the proposed procedure by considering different values of Υ that are mentioned in the literature.
The test statistic in Theorem 2 is based only on one time series observation (in a sense that if T > 2, then we can construct a test statistic for every value of t, but t = 1). However, it is not the most efficient way of using the time series information provided. Instead, all time series observations can be pooled into one test statistic to test the rank of: 18 For any fixed value of T , the y i,t y i,t−1 T term satisfies the sufficient conditions for the CLT, so that the results of Theorem 2 can be extended, with V N for this case given by: In the next section we use "rk-J" to denote the jacobian based cointegration test for y i,t y i,t−1 T . Until now we considered only the jacobian of Anderson and Hsiao (1982) moment conditions, however, for T > 2 further lags y i,t− j , can be used. Nevertheless, it is not clear that the use of lags j larger than j > 1 still ensures that, even in the effect stationary case, E[ y i,t y i,t− j T ] has reduced rank r if and only if rk Π = r . Moreover, the power of the test might be substantially affected by the choice of lags, as with any alternative close to the unit circle we encounter the weak instruments problem for any distanced lags. On the other hand, we can expect a better test power to the alternatives with substantially lower ρ(Φ).
Remark 3 If the model contains time effects λ t , the test statistic needs to be modified using variables in deviations from their cross-sectional averagesy Remark 4 One important advantage of the proposed test statistic is the additional flexibility while dealing with unbalanced panels. As long as for every individual i at least one y i,t y i,t−1 (t > 1) term is available, the test statistic can be computed. The only difference as compared to the balanced case is that an individual contribution to y i,t y i,t−1 T is no longer a simple averages with T − 1 terms, but has an individual specific number of observations T i − 1.

Remark 5
The testing procedure remains valid if, as suggested by Kleibergen and Paap (2006) One interesting special case is obtained when we set G N = I m and F −1 , as in this case we are testing the rank of the pooled OLS estimatorΠ. Even though the estimator itself is inconsistent (due to the presence of the unobserved heterogeneity), as we show in this paper, it can be used for estimation of rk Π 0 .

Discussions
In this section we summarize some of the underlying assumptions, and related problems, for the rk-J test.
Effect non-stationarity Recall that results in Theorem 2 are written in terms of the rank of y i,t y i,t−1 T rather than Π. This suggests that if one uses this rank test to perform the sequential procedure in testing the rank of Π the procedure is "conservative", i.e. as for some values of Υ , the jacobian can be of a reduced rank, even if Π is of full rank. However, the rank of y i,t y i,t−1 T can never be larger than the rank of Π. In such situations, the rk-J procedure controls the size of the test for Π, i.e. it rejects in at most α% cases, and it never gets larger than the nominal level, thus this test controls the size uniformly over (Σ, Υ ). However, for some combinations of the nuisance parameters (Σ, Υ ), the power of this test does not converge to 1 as N → ∞ even if rk(Π) > r 0 , and as a result such testing procedure lacks power. Thus, it is difficult to draw general conclusions about the properties of the rk-J test when one does not reject the null hypothesis and it is likely that y i,t process is not effect stationary.
Common dynamics assumption Throughout the paper we maintain the common dynamics assumption for η i . In the univariate case, it is known that if this assumption is satisfied the moment conditions are not relevant at unity, see a more detailed discussion in Bun and Kleibergen (2016). On the other hand, if the common dynamics assumption is violated, it is possible to have a full rank jacobian even for Φ 0 = I m , see the aforementioned paper and the discussion in Hayakawa and Nagata (2016). Hence, even if Π matrix is of reduced rank r < m the rank of the jacobian matrix can be of full rank m, when the common dynamics assumption is violated. In this case the rejection of the null hypothesis of the rk-J test, is not informative about the underlying rank of the Π. 19 Initialization For initialization we assumed that var ε i,0 is well defined irrespective of the time-series properties of the data. In the univariate setting, it is known that if e.g. Anderson and Hsiao (1982) moment conditions have a full rank jacobian matrix. Nevertheless, as discussed in Bun and Kleibergen 2016) this does not imply that φ parameter is identified. 20 Note that the initializations of this type would imply that the cross-sectional average of y i,t is not well defined for any t ≥ 0 which is a rather unrealistic assumption to make.
These issues cannot be underestimated in empirical work. However, at the same time we acknowledge that in order to obtain testing procedures that controls size uniformly 21 one would have to rely on procedures that are numerically challenging, i.e. subset inference using the continuously updated GMM estimator. We should also emphasize that most of the testing procedures for dynamic panel data (especially for persistent data) fail to guarantee uniform inference over the parameter space of autoregressive parameter and/or initialization of the initial condition. 22

Monte Carlo simulations
To the best of our knowledge only the BHP study provides results on cointegration analysis for panels with fixed T . 23 Hence, for the main building blocks of the finitesample studies performed in this paper we take the setups from BHP, but we provide an extended set of scenarios. Only bivariate panels are considered, thus the only null hypothesis we are testing is: For simplicity we use (DGP.2) for initialization: 19 However, in order to accommodate η i that does not satisfy the common dynamics assumption the data generating process of the initial condition needs to be modified. 20 For example if var ε i,0 = σ 2 /(1 − φ 2 ). 21 As e.g. in Andrews and Cheng (2012). 22 For example in the panel unit root testing literature this issue is prominent, see discussion in Moon et al. (2007) and also Westerlund (2016). 23 Other studies, like Mutl (2009) adapted setups of BHP. ε i,t for all i are generates as follows We assume that the error terms are normally distributed i.i.d. both across individuals and time with zero mean and variance-covariance matrix Σ (to be specified later). We set M = 50 24 and the number of replications B = 10000. We generate the individual heterogeneity (μ i ) using the exactly same procedure as in BHP: where we set vech Σ = (.05, .03, .05) . Following BHP, the variance-covariance matrix ofη i coincides with the corresponding variance-covariance matrix used in generating ε i,t . Before summarizing the design parameters for this Monte Carlo study, recall that Π can be rewritten as (for m = 2): We set λ = 0 to study the size of the test, while non-zero values of λ are used to investigate power. In particular, the following values of λ are considered: In order to reduce the dimensionality of the parameter space we assume that vectors α and β are of the following structure: As we discussed in Sect. 3.2 in the effect non-stationary case the particular choice of {Υ , Σ} and τ might substantially influence the performance of the test statistic. To address these concerns the following five choices of Υ are considered: 25 Υ = 0.5I 2 ; I 2 ; 1.5I 2 ; I 2 − Φ 10 ; . 85 .15 .00 .85 . 24 Results for M = 5 are qualitatively and quantitatively similar to the ones presented in this paper. 25 Later we use notation Υ (q) with q indicating the particular element of this set.
The choice of Υ (4) is motivated by the finite start-up assumption, so that the individual specific effects are accumulated only over 9 periods. The particular choice of S = 10 was rather arbitrary and is not empirically or theoretically motivated. 26 Comparing our setup to BHP, we can see that design 3 of BHP is achieved when α = −0.5 and λ = 0.0 (as they consider size). In order to match our designs with the empirical application, we also considered N = 750, however the results are qualitatively and quantitatively similar to N = 500, thus omitted. Other design parameters are also chosen to match some of the properties of the empirical application, as Υ (5) is based on the estimates in Arellano (2016) obtained from the bivariate panel of Spanish firm data.
In terms of the test power, we suspect that it should be decreasing with |λ|, with almost no power against alternatives with λ ≈ 0. However, it is very likely that for general Υ matrices the power curve might not be monotonic because λ not only controls the rank of Π but as well (indirectly) the eigenvalues of the E[u i,t−1 y i,t−1 ] matrix. Hence, for some specific choices of Υ we can observe the weak instruments problem of Anderson and Hsiao (1982) moment conditions that is not caused by the reduced rank of Π matrix.

Results
The results for all designs are summarized at the top part of Tables 6, 7, and 8 (θ = 0). All rejection frequencies are rounded up to two digits. Empty entries indicate maximal power of 1, 00.
General patterns First of all, we can observe that rejection frequencies are monotonically decreasing in |λ| for the vast majority of designs without spatial dependence. As we discussed in Sect. 3.2 this property should not be taken as granted for the rk-J test (as dependence on Φ is non-linear). For lower values of N the test tends to be undersized for T = 3 and oversized for T = 7. 27 In the effect stationary case τ does not play substantial role and only affects the V matrix, but we can still observe that higher value of τ is associated with slightly lower power. For N = 500, the rk-J test has notable power even when λ is very close to 0. For instance, all rejection frequencies in the effect stationary designs at λ = 0.005 are above 30% and 25% for α = −0.5 and α = −0.1 respectively. In the vast majority of cases with size distortions being of similar magnitude, the test power for α = −0.5 tends to be higher than for α = −0.1.
Effect non-stationarity and non-monotonic power curves First, we consider rejection frequencies for Υ = 0.5 × I m as this case is most exceptional in terms of observed patterns. In this case we observe power curves that are not monotonic for α = −0.1 (especially for N = 250) and sharply decreasing for α = −0.5 if τ = 5 and T = 3. It can be intuitively explained as in this case the effect non-stationarity term in E[ y i,t y i,t−1 ] is negative, driving the whole expression towards the zero matrix (recall the analysis in Hayakawa 2009 for the univariate case). Thus, we have a weak instrument problem under the alternative hypothesis that is not induced by cointegration. 28 By varying λ parameter we directly vary the relative contributions of time invariant and time varying parts of the variance components in var y i,t . For larger values of |λ| the time invariant part is more pronounced, resulting in substantial effects of the "negative" effect stationarity. On the other hand, for |λ| ≈ 0 the idiosyncratic part is dominant and there is no substantial effects of the "negative" effect non-stationary initialization.
Remark 6 This non-monotonicity is further illustrated in Tables 6, 7, where we show how the minimum eigenvalue of the jacobian matrix changes for different nuisance parameters (for very larger N ). Those patters resemble power curves of the rk-J test as presented in Fig. 1.
As it can be expected, the results for Υ = 1.5 × I m are more straightforward. In this case the power curves are monotonic, and rejection frequencies are uniformly dominating the ones from effect stationary case irrespective of other design parameters. Results for Υ (4) seem to combine the properties of both Υ (3) and Υ (1) . 29 Finally, the results of Υ (5) are somewhat in between those of Υ (1) and Υ (2) , but are slightly closer to Υ (2) . It serves as an indication that the off-diagonal element in Υ (5) is not of any great importance (given the choice of other design parameters).
Remark 7 In this paper, we do not provide extensive results for the TML estimator of Binder et al. (2005). The main reason for this (besides theoretical problems discussed in Sect. 2.2) is possibly bimodal log-likelihood function (see e.g. Calzolari and Magazzini 2012;Bun et al. 2016;Juodis 2016). For model with stable dynamics, Juodis (2016) presents several alternatives how one can choose the maximizer of the log-likelihood function from the set of local minimizers. Unfortunately, no results are available for non-stationary dynamics analyzed in this paper. Thus, in order to avoid the situation in which unintentionally test procedure based on the TML estimator performs suboptimally, we present only some limited results, see Table 5. Results suggest that for alternatives close to the null hypothesis LR test has low power, as the critical value from the χ 2 (1) distribution is too large. On the other hand, for some very distant alternative (where the rk-J test struggles to reject the null hypothesis), LR test has sizeable power.
Remark 8 As a robustness check in 6, we also consider model with spatial dependence in the error terms. Evidence of the uniform upward shift in the size can be observed when designs with spatial dependence are considered. 28 Some preliminary MC results, not presented in this paper suggest that effect of τ in this setup is notmonotonic. In the sense that higher values of τ lead to increase of power rather than further decrease. At least for this particular design it seems that τ = 5 represents the close to worst possible scenario as minimum is reached for τ ≈ 6.2. 29 From Υ (1) some non-monotonicities are inherited. Apart from that, the superior test power properties (as compared to the effect stationary case) of Υ (3) are dominant. This combined behavior is due to the fact that Υ (4) is changing with λ. In designs with λ substantially lower than 0 we have Υ (4) ≈ I m , consecutively the weak instrument problem under alternative is less pronounced. In this section, we analyze the Spanish firm panel dataset covering 1983-1990 of 738 manufacturing companies from Alonso-Borrego and Arellano (1999). This datasets constitutes a balanced panel of manufacturing companies recorded in the database of We construct a bivariate PVAR(1) model with logarithms of employment and wages as dependent variables. Table 1 contains year specific descriptive statistics for these two variables. 30 Given that cross-sectional means for both variables differ substantially (especially that of wages) between the beginning and the end points, we follow other papers that considered this dataset (e.g. Arellano 2016) and include the time effects in the model, i.e. we consider variables in their deviations from the corresponding year specific cross-sectional averages. The sensitivity to the cross-sectional demeaning is discussed later in this section.

Results
In contrast to the previous studies that used this data, we investigate the time-series properties of this dataset in a greater detail. In particular, previous studies assumed that GMM and ML estimators are well-behaved, i.e. unit roots and cointegration were excluded a priori. However, estimation results in Arellano (2016) indicate that some estimated parameter values can be close to unit circle, thus non-stationary behaviour cannot be excluded beforehand. In order to elaborate on those observations, we consider a slightly larger set of estimators to obtain point estimates for Φ that are valid under different sets of assumption. The results in Table 2 are in line with those in Arellano (2016), with many close to unity point estimates of φ 11 . The estimates for φ 22 , on the other hand, are further away from unity, suggesting that both variables can be potentially cointegrated. In order to investigate for possible cointegration in this 30 For a more detailed description of the data, please refer to Alonso-Borrego and Arellano (1999). Here "AB(·)" and "Sys(·)" are the estimators of Arellano and Bond (1991) and Blundell and Bond (1998), respectively. The numbers in brackets indicate, whether these are "two-step" or "one-step" estimates. "FE" denotes the fixed effects estimator. "FEBC" (from top to the bottom) are the bias-correcting fixed effects estimators of Hahn and Kuersteiner (2002), Kiviet (1995) (using "Sys(2)" as the plug-in estimator), Dhaene and Jochmans (2015), Bun and Carree (2005), Juodis (2013). "TMLE" are the Transformed Maximum likelihood estimators with and without rank restrictions imposed on Φ dataset, we make use of the rk-J procedure that was introduced earlier in this paper. 31 Specifically, we test if there is a single cointegrating relationship between firms level employment and wages, i.e. H 0 : r 0 = 1. 32 First, we apply the rk test of Kleibergen and Paap (2006) directly to GMM esti-matesΠ. We restrict the set of GMM estimators to two step estimators that are also presented in BHP: "AB-GMM" stands for the estimator of Arellano and Bond (1991), while "Sys-GMM" is the estimator of Blundell and Bond (1998) that incorporates moment conditions in levels. Second, the LR tests based on the transformed maximum likelihood function of BHP (LR-TMLE) and maximum likelihood function of Arellano (2016) (as mentioned in Remark 1), (LR-RMLE), are considered. Finally, the "rk-J" test of Sect. 3.2 is considered. Under H 0 : rk Π 0 = 1, if no singularities in the corresponding asymptotic distributions are present, all tests have a χ 2 (1) limit. Note that we present results for "AB-GMM" for informal comparison only, as under H 0 : r 0 = 1 this estimator is not consistent. Results are summarized in Table 3.
From Table 3 we can see that only the rk-J test based on the Anderson and Hsiao (1982) moment conditions rejects H 0 . Results for system GMM estimator are mixed, as based on Windmeijer (2005) corrected standard errors the null hypothesis is not rejected, while it is rejected when using the conventional two-step standard errors. Numerous reasons might account for differences in conclusions. First of all, we suspect 31 Other testing procedures are described below. 32 We focus only on testing r = 1 vs. r = 2, as e.g. using LR test based on the TML estimator the H 0 : r 0 = 0 is rejected against the alternative H A : r 0 = 1 with the value of the test statistic equal to 117.561. This value is substantial even if the χ 2 (1) does not provide a correct asymptotic approximation.  Table 3 might be the low power of cointegration test used directly on the estimate of Π. Now we turn our attention to the likelihood ratio tests. Based on analytical results in this paper for T = 2 we can suspect that the likelihood procedures under H 0 of cointegration lack power for close alternatives (recall limited MC results in Table 5), as χ 2 (1) is a poor approximation of the finite sample distribution. Furthermore, we know that both likelihood methods are robust to violations of mean stationarity, but are not so to time-series heteroscedasticity. Thus, we can not rule out the possibility that it can be one of the reasons for divergence in conclusions. 34

Sub-sample analysis
In the previous section we investigated the relationships between the firm level employment and wages in the model estimated using the full length of the dataset. Using the rk-J procedure, no significant statistical evidence was found favouring cointegration between the two variables. In this section, we investigate the sensitivity of this conclusion to smaller values of T by means of the analysis over sub-samples of the original data. We also check the sensitivity towards inclusion/non-inclusion of the time effects in the model. The results of this section are summarized in Table 4, where in total 6 different sub-samples are considered. First of all, we observe that the non-inclusion of time effects leads to an increase of the test statistic in all sub-samples. The difference is especially pronounced for all subsamples with 1990 as the final year. Overall, irrespective of the time span considered the rk − J statistic rejects the null hypothesis if the time-effects are not included, i.e. data is not cross-sectionally demeaned before estimation. 33 However, this testing procedure cannot be used if series are cointegrated. 34 Arellano (2003) presents some evidence of time-series heteroscedasticity in this dataset. The same cannot be generally said when data is cross-sectionally demeaned. Note how the value of test statistic increases as T increases for sub-sample ending in 1990. In particular, for T = {4; 5} the null hypothesis is not rejected at any conventional significance level. This behavior emphasizes the value of additional time-series observations and possible lack of power for small values of T . As it can be seen from Table 4 the observations for 1983 are especially informative about the properties of the bivariate system, as for all sub-samples starting in 1983 the test statistic always rejects the null hypothesis.
Overall, omission of time-effects from the model does not affect the conclusions from Sect. 5.2. However, a moderate amount of time variation in the magnitude of test statistics suggests that this conclusion is sensitive to different estimation horizons.

Conclusions
In this paper, we study the properties of the standard Anderson and Hsiao (1982) moment conditions in a PVAR(1) for cointegrated processes. Under the assumptions similar to Binder et al. (2005) we show that these moment conditions are of reduced rank if the process is cointegrated. Based on this observation we propose a rank based test for the null hypothesis of cointegration. We prove that testing procedure in Binder et al. (2005) is invalid due to the singularity of the hessian matrix for persistent data. Monte Carlo results suggest that for most designs, the new test is reasonably sized and has good power properties but might exhibit non-monotonic power curves for models with substantial effect non-stationarity. We apply our testing procedure to the Spanish manufacturing data of Alonso-Borrego and Arellano (1999) and, unlike the test of BHP, we find no evidence of cointegration. and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix Proofs: Transformed maximum likelihood Notation
The commutation matrix K a,b is defined such that for any First, we define a set of new auxiliary variables, which are handy during the derivations of differentials: Results from Binder et al. (2005) and Juodis (2016) Define then the score vector associated with the log-likelihood function (8) is given by: For 0 < r < m we have the following result.

Corollary 1 Let Assumptions TML be satisfied. Then the restricted score vector associated with the log-likelihood function (8) under cointegrating restrictions is given by:
Proof of this corollary is a a trivial extension of the corresponding result in Juodis (2016) for the unrestricted model. In the special case where m = 2 and r = 1, δ 1 is a scalar while α is a [2 × 1] vector implying that the corresponding entries of the ∇ r (κ r ) vector do not have a vec (·) operator in them.

Identification with unit roots
Proof of Theorem 1 In order to evaluate the expected hessian matrix, we need to separately calculate the expected value for every term presented in the formula for the second differential at the DGP values κ 0 . First of all, we note that E hence the contribution to the hessian matrix of the first four terms in the expression for d 2 (κ) is zero. Following, e.g. Juodis (2013) it can be shown: Due to the fact that exact expressions for other terms are rather messy in the general case, they are derived only for the particular case where Φ 0 = I m . Under this assumption and assumptions of Binder et al. (2005), it then follows that Θ 0 = Σ 0 . The DGP in this case under restriction of common dynamics, simplifies to: At first we consider expectations of the N N (κ 0 ) term which we can evaluate based on the general result derived in Lemma A.2 Juodis (2016): It then similarly follows from the general result in Lemma A.1 Juodis (2016) that: What is left is to evaluate the term involving R N : In Juodis (2016) it is shown that the hessian matrix is of the following form (if we neglect terms that evaluated at the true values κ 0 have expectation O): Plugging in these expressions into the formula for the hessian and evaluating it at the true value κ 0 : Using the formula for the determinant of partitioned matrix, that: where ∝ is the proportionality sign. For future reference, define Thus, H in this case given by: Then by means of the matrix determinant Lemma: where in the first line we used the fact that K m = K −1 m and the second line follows from the fact that | K m | = (−1) 0.5m(m−1) , while |I m 2 − K m | = 0 follows trivially from the fact that rk (I m 2 − K m ) = 0.5m(m − 1). Hence we have proved that the H matrix is not invertible.

Cointegration and T = 2
Proof of Proposition 1 Using derivations of previous section it follows that: Corresponding H matrix: We can express Θ 0 = Σ 0 + aa with a being of rank r . Then: Thus using basic inequalities for the rank of matrix sum we can deduce that rk H is no greater than 0.5m(m − 1) + r 2 , due to the full rank of (2Σ 0 ⊗ Σ 0 + Σ 0 ⊗ aa ) matrix.

Robustness check: spatial dependence
To allow for possible spatial dependence, the ε i,t for all t are generated with Spatial MA process: The θ parameter controls the degree of cross-sectional dependence between units. For θ = 0 we have i.i.d. dataset, while for θ = 0 the cross-sectional units are weakly correlated. To illustrate the effect of spatial dependence we consider one value for θ , namely θ = 0.5. Spatial correlation matrix W N with a typical (i, j) element given by ω i, j is assumed to be 1 ahead -1 behind circular, so that every individual i is directly linked only with individuals i − 1 and i + 1. 35 The particular choice of the spatial matrix W N is motivated by the study in Baltagi et al. (2007), where in the context of the panel unit root testing it is shown that the tests are mostly distorted for this choice of spatial matrix. 36 The results are presented at the bottom of Tables 6, 7, and 8.
Evidence of the uniform upward shift in the size can be observed when designs with spatial dependence (θ = 0.5) are considered. This upward movement does not come as a surprise because similar patterns have been documented in the panel unit root testing literature. However, the same conclusion can not be reached regarding the test power, as for most scenarios it changes marginally and does not show any clear patterns in terms of magnitude and direction. More importantly, major size distortions do not disappear for N = 500, thus as can be expected the fact that we use the variancecovariance matrix that ignores the presence of spatial dependence has a pronounced result. On the other hand, the fact that by ignoring the spatial dependence the rank based test statistic is only mildly oversized, suggests that the procedure developed in this paper is relatively robust to deviations from the i.i.d. assumption.

Tables
See Tables 5, 6 , 7, 8, 9, 10, and 11. 35 The circle is closed by connecting i = 1 with i = N . 36 For a graphical illustration see Figure 2 of the aforementioned paper. is used for initialization and N = 750. Empty entries correspond to 100% empirical rejection frequencies.
As starting values for TML estimators we use the bias corrected Fixed Effects estimators as define in Table 2   .71 .08 See Table 8  Table 10 Monte  .79 .09 See Table 8