Panel Data Cointegration Testing with Structural Instabilities

Abstract

Spurious regression analysis in panel data when the time series are cross-section dependent is analyzed in the article. The set-up includes (possibly unknown) multiple structural breaks that can affect both the deterministic and the common factor components. We show that consistent estimation of the long-run average parameter is possible once cross-section dependence is controlled using cross-section averages in the spirit of Pesaran’s common correlated effects approach. This result is used to design individual and panel cointegration test statistics that accommodate the presence of structural breaks that can induce parameter instabilities in the deterministic component, the cointegration vector and the common factor loadings.

KEYWORDS:

• Common factors
• Cross-section dependence
• Spurious regression
• Structural breaks

1 Introduction

The literature on nonstationary panel data has experienced important developments in recent times, analyzing the properties of parameter estimation under both spurious regression and cointegration. One interesting feature of working in a nonstationary panel data framework is that it is possible to obtain consistent estimates of statistical relationships regardless of whether there is cointegration or spurious regression—see Phillips and Moon (1999). The idea behind the consistent estimation of (the long-run average) parameters in a spurious panel relies on the fact that panel units add information that lead to detect a stronger overall signal than that of the pure time series case—see Phillips and Moon (2000). This article adds to this body of the panel data literature by defining a model specification under spurious regression that allows for the presence of multiple structural breaks affecting the level, trend and/or the slope parameters of the model, as well as in the factor structure generating the cross-section dependence. Panel data units are assumed to be cross-section dependent which is captured by an approximate common factor model.

There are a number of contributions in the literature that are relevant for our proposal. Gengenbach, Urbain, and Westerlund (2016) and Banerjee and Carrion-i-Silvestre (2017) design panel data cointegration tests with cross-section dependence, but without structural breaks. Although the two papers differ in the degree of heterogeneity that is allowed for in the potential cointegration vector—the former paper considers a heterogeneous cointegration vector across the units of the panel, whereas the latter assumes homogeneity—both follow the spirit of the common correlated effects (CCE) estimator advocated in Pesaran (2006) and show that the common factors that drive the cross-section dependence can be proxied through the use of cross-section averages of the variables that appear in the model. CCE-methodology offers the advantage of avoiding the need to estimate the number of common factors $(m^0)$, although m⁰ is required to be bounded by the number of observable stochastic processes $(1+k)$ that appear in the model specification—the so-called rank condition, $m^0\le 1+k$.

Westerlund and Edgerton (2008) and Banerjee and Carrion-i-Silvestre (2015) define panel data cointegration statistics accounting for multiple structural breaks—that can affect the deterministic component and/or the slope parameters of the model—and cross-section dependence—modeled with a common factor model. The framework of these papers allows for complete heterogeneity in the sense that all parameters are panel unit specific—that is, the break dates, the deterministic component and the slope parameters are allowed to be heterogeneous across units. The common factors are estimated using principal components as in Bai and Ng (2004) so that it is possible to estimate the number of common factors using information criteria. However, it is worth noting that Westerlund and Edgerton (2008) assume that all common factors are I(0) stationary, whereas Banerjee and Carrion-i-Silvestre (2015) allow for a combination of I(0)/I(1) common factors.

The contribution of the present article can be summarized as follows. First, the article shows that it is possible to obtain a consistent pooled estimate of the long-run average coefficient, which captures a statistical relationship among non-individually cointegrated variables—see Phillips and Moon (1999). The model specification allows for the presence of multiple structural breaks that can affect the long-run average coefficient. This result is of special relevance from an empirical point of view, since it implies that practitioners can obtain consistent estimates of the slope parameters using a panel data pooled estimator regardless of whether the variables define a cointegration relationship. This is in direct contrast to the findings in the papers mentioned in the previous paragraph that assume heterogeneous potential cointegration vectors and estimate the parameters unit by unit, which delivers inconsistent estimators of the parameters under spurious regression. Second, the fact that the panel cointegration statistic is based on the pooled estimation of the potential cointegration vector, removes the presence of the stochastic regressors in the limiting distribution of the cointegration statistic—that is, we move from a panel data cointegration analysis framework to a panel data unit root one. This produces a spin-off contribution, since the statistic that is proposed here extends the Pesaran (2007) panel data unit root test by allowing for structural breaks in the same way as Perron (1990) and Kim and Perron (2009) extended the Dickey-Fuller framework in time series analysis. Third, our setup also deals with the case where the factor loadings can be affected by structural breaks, a situation that has not been previously investigated in panel data cointegration testing analysis. This adds more flexibility and increases the interest of the model from an empirical point of view. Fourth, the common factors are estimated using the CCE approach which avoids the need to estimate the number of common factors as long as the rank condition is satisfied. While the CCE approach does not pose any problem for the computation of the pooled estimator of the slope coefficients, we show there are potentially important issues at the testing stage. This leads us to suggest a conservative testing strategy that mitigates the need to know both the number and order of integration of the common factors. Finally, the article proves that in some cases it is possible to obtain a consistent estimate of the break fraction parameters under spurious regression when the break dates are unknown. This is an interesting result per se, adding to the results obtained in Phillips and Moon (1999) concerning the consistent estimation of the long-run average coefficient under spurious regression.

Taken together, our article provides a panel cointegration testing framework that explicitly models the effects of structural breaks. Westerlund (2018) had argued that an approximate common factor model—estimated using the CCE approach—can result in a simple device that might capture unattended structural breaks. In principle this is an interesting feature that could act as a model-safeguarding device, avoiding the need to define a particular deterministic specification, at the risk of violating the rank condition. We argue however that if the presence of breaks is a reasonable feature of the data generation process, it would be better to undertake explicit modeling of these breaks and free the common factors to capture other potential misspecification errors and dependences that might affect the model.

The article is organized as follows. Section 2 describes the model, Section 3 defines the pooled CCE estimator of the slope parameters and deals with the estimation of the structural break dates. The consistency of the pooled CCE estimator provides the ground upon which individual and panel data cointegration statistics are proposed in Section 4. Section 5 conducts an extensive Monte Carlo simulation experiment to analyse the finite sample performance of the statistics that have been proposed in the article. An empirical illustration that focuses on house prices and per capita disposable income for the U.S. states is conducted in Section 6. Finally, Section 7 concludes. Supplementary material is collected in the appendices with the proofs, tables of critical values, and tables and figures that summarize the finite sample performance of the statistics. A GAUSS code program is available to implement the proposal that is designed in the article.

2 The model

Let $Y_{i,t}={(y_{i,t},x_{i,t}^{\prime })}^{\prime }$ be a $(1+k)$ -vector of $I(1)$ stochastic processes with data generating process (DGP): (1) $$Y_{i,t}={\Theta }_i{D}_t+{\pi }_{i,t}{F}_t+U_{i,t}$$(1) (2) $$\left(I-L\right)F_t=v_t$$(2) (3) $$\left(I-L\right)U_{i,t}=e_{i,t},$$(3) where Θ_i is a $((1+k)\times (J^0+1))$ matrix of coefficients with rows defined as ${\theta }_{i,l}=({\theta }_{i,l,1,1},\dots ,{\theta }_{i,l,J^0+1,1},{\theta }_{i,l,1,2},\dots ,{\theta }_{i,l,J^0+1,2}),l=1,\dots ,(1+k),D_t={(1,D{U}_t^{\prime },t,D{T}_t^{\prime })}^{\prime },D{U}_t={(D{U}_{1,t},\dots ,D{U}_{J^0,t})}^{\prime },D{T}_t={(D{T}_{1,t},\dots ,D{T}_{J^0,t})}^{\prime }$ with $D{U}_{j,t}=1(t>T_j^0),D{T}_{j,t}=(t-T_j^0)1(t>T_j^0),1(·)$ the indicator function, $T_j^0=\lfloor {\lambda }_j^{0}T\rfloor$ the jth break date, ${\lambda }_j^0=T_j^0/T$ the break fraction and $\lfloor ·\rfloor$ the integer part, $j=1,\dots ,J^0$, with the convention that $T_0^0=0$ and $T_{J^0+1}^0=T$—that is, ${\lambda }_0^0=0$ and ${\lambda }_{J^0+1}^0=1$. Throughout the article, the superscript “0” indicates the true value of the corresponding figure. The deterministic component covers two variants: (a) a constant with level shifts—henceforth, constant case—and (b) a time trend with both level and slope shifts—in what follows, time trend case. The set of generic break fraction parameters vector is defined as $\Lambda =\{\lambda =({\lambda }_1,\dots ,{\lambda }_{J^0})|ϵ<{\lambda }_1<\cdots <{\lambda }_{J^0}<1-ϵ$ and $|{\lambda }_j-{\lambda }_{j+1}|>ϵ\}$, with ϵ denoting the amount of trimming that is specified—popular choices are $ϵ\in \{0.1,0.15,0.2,0.25\}$. At this stage, the model set-up assumes that the number and position of the common structural breaks are known, but the procedure that can be implemented to estimate them is discussed below.

To begin with, cross-section dependence in the panel is defined by F_t which is a m⁰-vector of observable common factors and ${\pi }_{i,t}={\pi }_{i,j},T_{j-1}^0<t\le T_j^0$, the $((1+k)\times m^0)$ matrix of factor loadings. Denote by $K=\sigma (F_0,\dots ,F_t,\dots )$ the sigma field generated by the sequence ${\left\{F_t\right\}}_{t=0}^{\infty }$ so that, conditionally on K, $U_{i,t}={(U_{y_i,t},U_{x_i,t}^{\prime })}^{\prime }$ are independent across i—see Urbain and Westerlund (2011). Let M<∞ be a generic positive number, not depending on T and N, and define the Euclidean norm of a generic matrix A as $\left||A\right||=\mathit{\text{trace}}{(A^{\prime }A)}^{1/2}$. The vector of stochastic processes $V_{i,t}={(v_t^{\prime },e_{i,t}^{\prime })}^{\prime },i=1,\dots ,N,t=1,\dots ,T$, is assumed to satisfy the following assumptions—see Bai and Ng (2004) and Banerjee and Carrion-i-Silvestre (2015).

Assumption 1.

(i) $v_t=C_j^F(L)w_t,w_t\sim \text{iid}(0,{\Sigma }_{w_j}),E\left||w_t\right|{|}^4\le M,T_{j-1}<t\le T_j$, and (ii) $\text{var}(\Delta F_t)={\sum }_{l=0}^{\infty }{C}_{j,l}^F{\Sigma }_{w_j}{C}_{j,l}^{F\prime }>0,T_{j-1}<t\le T_j$, (iii) ${\sum }_{l=0}^{\infty }l\left||C_{j,l}^F\right||<M$; and (iv) $C_j^F(1)$ has rank $m_1^0\forall j,0\le m_1^0\le m^0,j=1,\dots ,J^0+1$.

Assumption 2.

(i) For each i, $e_{i,t}=C_{i,j}^U(L){\epsilon }_{i,t},{\epsilon }_{i,t}\sim \text{iid}(0,{\sigma }_{{\epsilon }_{i,j}}^2),E\left||{\epsilon }_{i,t}\right|{|}^8\le M,{\sum }_{l=0}^{\infty }l\left||C_{i,j,l}^U\right||<M,{\omega }_{i,j}^2=C_{i,j}^U{(1)}^2{\sigma }_{{\epsilon }_{i,j}}^2>0$; (ii) $E({\epsilon }_{i,t}{\epsilon }_{l,t}^{\prime })={\tau }_{i,l}$ with ${\sum }_{i=1}^N\left||{\tau }_{i,l}\right||\le M$ for all l, $T_{j-1}<t\le T_j$; (iii) $E\left||\frac{1}{\sqrt{N}}{\sum }_{i=1}^N\left[{\epsilon }_{i,s}{\epsilon }_{i,t}^{\prime }-E({\epsilon }_{i,s}{\epsilon }_{i,t}^{\prime })\right]\right|{|}^4\le M$, for every (t, s); and (iv) $C_{i,j}^U(1)$ is positive definite almost surely for all i and j, $j=1,\dots ,J^0+1$.

Assumption 3.

The errors ${\epsilon }_{i,t}$, w_t and loadings π_i are three mutually independent groups across i, t and $(1+k)$ dimensions.

Assumption 4.

$E\left||F_0\right||\le M$, and for every $i=1,\dots ,N,E\left||U_{i,0}\right||\le M$.

Assumption 5.

As $T_j^0\to \infty$ and $T\to \infty$, $T_j^0/T\to {\lambda }_j^0,j=1,\dots ,J^0$, with $0<{\lambda }_1^0<\cdots <{\lambda }_{J^0}^0<1$.

Assumptions 1–5 ensure that the central limit theorem (CLT) holds for $V_{i,t}$, so that, for the most general specification that is considered in this article, we have: $$T^{-1/2}\sum _{t=1}^{\left[\mathit{\text{rT}}\right]}{V}_{i,t}\Rightarrow \sum _{k=1}^{j-1}\left[C_{i,k}\left(1\right)-C_{i,k+1}\left(1\right)\right]W_i\left({\lambda }_k^0\right)$$ $$\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}+C_{i,j}\left(1\right)W_i\left(r\right);r\in \left[{\lambda }_{j-1}^0,{\lambda }_j^0\right],$$

as $T\to \infty$ for all i, where $W_i(r)$ is a $(m^0+1+k)$-vector of standard Brownian motion on $r\in \left[{\lambda }_{j-1}^0,{\lambda }_j^0\right]$, and $${\Omega }_{i,j}=C_{i,j}\left(1\right)C_{i,j}{\left(1\right)}^{\prime }=\left[\begin{array}{cc}{C}_j^F\left(1\right)C_j^F{\left(1\right)}^{\prime }& C_j^F\left(1\right)C_{i,j}^U{\left(1\right)}^{\prime }\\ C_{i,j}^U\left(1\right)C_j^F{\left(1\right)}^{\prime }& C_{i,j}^U\left(1\right)C_{i,j}^U{\left(1\right)}^{\prime }\end{array}\right],$$is the long-run conditional covariance matrix with expected value $E({\Omega }_{i,j})={\Omega }_j,j=1,\dots ,J^0+1$, the so-called long-run conditional average covariance matrix—see Lemma 1 in the appendix. The matrix ${\Omega }_{i,j}$ can be partitioned to define ${\Omega }_{U_{i,j}{U}_{i,j}}=C_{i,j}^U(1)C_{i,j}^U{(1)}^{\prime },{\Omega }_{U_{x_{i,j}}{U}_{x_{i,j}}}=C_{i,j}^{U_x}(1)C_{i,j}^{U_x}{(1)}^{\prime }$ and ${\Omega }_{U_{x_{i,j}}{U}_{y_{i,j}}}=C_{i,j}^{U_x}(1)C_{i,j}^{U_y}{(1)}^{\prime },j=1,\dots ,J^0+1$.

The model considers the case where $x_{i,t}$ are assumed to be either cross-section independent—imposing all, but the first, rows of ${\pi }_{i,t}$ to be equal to zero—or cross-section dependent with dependence driven by F_t. Furthermore, it is possible to assume that the set of observable common factors affecting the endogenous variable $y_{i,t}$ is different from those affecting $x_{i,t}$, a situation that is covered if we define ${\pi }_{i,t}$ to be a block-diagonal matrix.

Despite the presence of the operator $(I-L)$ in (2), F_t can be I(0), I(1), or a combination of both, depending on the rank of $C_j^F(1),j=1,\dots ,J^0+1$. Let $m_0^0$ and $m_1^0$ be the number of I(0) and I(1) common factors, respectively, so that $m^0=m_0^0+m_1^0$. If $C_j^F(1)=0$, then F_t is I(0) and $m_0^0=m^0$. If $C_j^F(1)$ is of full rank, then each component of F_t is I(1) and $m_1^0=m^0$. If $C_j^F(1)\ne 0$, but not full rank, then some components of F_t are I(1) and some are I(0). Following Bai, Kao and Ng (2009) and Banerjee and Carrion-i-Silvestre (2015), the inclusion of common factors in the potential long-run relationship implies a change in the standard concept of cointegration. The usual definition of cointegration among $Y_{i,t}$ requires F_t to be I(0), so that $Y_{i,t}$ captures all the common stochastic trends. However, allowing F_t to be I(1) is also relevant from an empirical point of view since F_t might be accounting for effects that are not captured by $Y_{i,t}$ alone. Another interesting feature is that (1) specifies time dependent factor loadings so that the structural breaks can also affect the way in which the common factors affect panel units. Finally, the DGP also allows for the possibility that the structural breaks do not affect all elements in $Y_{i,t}$, and also covers intermediate situations where, for instance, the deterministic component and/or the loadings do not change across regimes for some of the elements in $Y_{i,t}$.

The most general model specification that is admitted in our framework is (4) $$y_{i,t}={\mu }_{i,j,1}^y+{\mu }_{i,j,2}^y(t-T_{j-1}^0)+x_{i,t}^{\prime }{\beta }_j+{\eta }_{i,j}^y{F}_t+{\xi }_{i,t}^y$$(4) (5) $$x_{i,t}={\mu }_{i,j,1}^x+{\mu }_{i,j,2}^x(t-T_{j-1}^0)+{\eta }_{i,j}^x{F}_t+{\xi }_{i,t}^x,$$(5)

${\mu }_{i,j,k}={({\mu }_{i,j,k}^y,{\mu }_{i,j,k}^{x\prime })}^{\prime }$, $k\in \left\{1,2\right\},{\eta }_{i,j}={({\eta }_{i,j}^{y\prime },{\eta }_{i,j}^x)}^{\prime },T_{j-1}^0<t\le T_j^0$, $j=1,\dots ,J^0+1$. Combining (4) and (5) we obtain: $$\begin{array}{c}\left(\begin{array}{c}{y}_{i,t}\\ x_{i,t}\end{array}\right)=\left(\begin{array}{cc}{\beta }_j^{\prime }{\mu }_{i,j,1}^x+{\mu }_{i,j,1}^y& {\beta }_j^{\prime }{\mu }_{i,j,2}^x+{\mu }_{i,j,2}^y\\ {\mu }_{i,j,1}^x& {\mu }_{i,j,2}^x\end{array}\right)d_t\\ \hspace{1em}+\left(\begin{array}{c}{\beta }_j^{\prime }{\eta }_{i,j}^x+{\eta }_{i,j}^y\\ {\eta }_{i,j}^x\end{array}\right)F_t+\left(\begin{array}{c}{\beta }_j^{\prime }{\xi }_{i,t}^x+{\xi }_{i,t}^y\\ {\xi }_{i,t}^x\end{array}\right),\end{array}$$that is, (6) $$\left(\begin{array}{c}{y}_{i,t}\\ x_{i,t}\end{array}\right)=\left(\begin{array}{cc}{\mu }_{i,j,1}^{y\ast }& {\mu }_{i,j,2}^{y\ast }\\ {\mu }_{i,j,1}^x& {\mu }_{i,j,2}^x\end{array}\right)d_t+\left(\begin{array}{c}{\pi }_{i,t}^{y\ast }\\ {\pi }_{i,t}^x\end{array}\right)F_t+\left(\begin{array}{c}{\xi }_{i,t}^{y\ast }\\ {\xi }_{i,t}^x\end{array}\right),$$(6)

$d_t={(1,(t-T_{j-1}^0))}^{\prime }$, which corresponds to (1) given that ${\theta }_{i,1,j,k}={\theta }_{i,j,k}^y={\sum }_{l=1}^j{\mu }_{i,l,k}^{y\ast },{\mu }_{i,l,k}^{y\ast }={\beta }_l^{\prime }{\mu }_{i,l,k}^x+{\mu }_{i,l,k}^y,{\theta }_{i,j,k}^x={\sum }_{l=1}^j{\mu }_{i,l,k}^x$, $k\in \left\{1,2\right\},{\pi }_{i,t}={({\pi }_{i,t}^{y\ast \prime },{\pi }_{i,t}^{x\prime })}^{\prime },{\pi }_{i,t}^{y\ast }={\beta }_j^{\prime }{\eta }_{i,j}^x+{\eta }_{i,l}^{y\prime },{\pi }_{i,t}^x={\eta }_{i,j}^x$, and $U_{i,t}={({\xi }_{i,t}^{y\ast \prime },{\xi }_{i,t}^{x\prime })}^{\prime },{\xi }_{i,t}^{y\ast }={\sum }_{l=1}^j{\beta }_l^{\prime }{\xi }_{i,t}^x+{\xi }_{i,t}^y,T_{j-1}^0<t\le T_j^0,j=1,\dots ,J^0+1$. Although the model has been written as a pure structural change model in which all parameters change, we can impose parameter constraints in (4) to obtain partial structural change specifications. Our article distinguishes six different model specifications depending on whether the structural breaks affect the deterministic component, the slope parameters and/or the factor loadings. Table 1 summarizes all the cases covered in the article. Model A assumes that the structural breaks affect only the deterministic component, Model B considers that the structural breaks affect both the deterministic component and the slope parameters and, finally, Model C permits the structural breaks to affect the deterministic component, the slope parameters and the factor loadings. For each model specification we consider the two variants of the deterministic component—designated by a 1 for the constant case and by a 2 for the time trend case. This defines Models A1, A2, B1, B2, C1, and C2.

Table 1 Model specifications.

Model	Deterministic component		Slope	Factor
	Constant	Trend	parameters	loadings
A1	${\mu }_{i,j,1}\ne {\mu }_{i,j+1,1}\forall j$	${\mu }_{i,j,2}=0\forall j$	${\beta }_j={\beta }_0\forall j$	${\eta }_{i,j}={\eta }_{i,\mathrm{ }0}\forall j$
A2	${\mu }_{i,j,1}\ne {\mu }_{i,j+1,1}\forall j$	${\mu }_{i,j,2}\ne {\mu }_{i,j+1,2}\forall j$	${\beta }_j={\beta }_0\forall j$	${\eta }_{i,j}={\eta }_{i,\mathrm{ }0}\forall j$
B1	${\mu }_{i,j,1}\ne {\mu }_{i,j+1,1}\forall j$	${\mu }_{i,j,2}=0\forall j$	${\beta }_j\ne {\beta }_{j+1}\forall j$	${\eta }_{i,j}={\eta }_{i,0}\forall j$
B2	${\mu }_{i,j,1}\ne {\mu }_{i,j+1,1}\forall j$	${\mu }_{i,j,2}\ne {\mu }_{i,j+1,2}\forall j$	${\beta }_j\ne {\beta }_{j+1}\forall j$	${\eta }_{i,j}={\eta }_{i,\mathrm{ }0}\forall j$
C1	${\mu }_{i,j,1}\ne {\mu }_{i,j+1,1}\forall j$	${\mu }_{i,j,2}=0\forall j$	${\beta }_j\ne {\beta }_{j+1}\forall j$	${\eta }_{i,j}\ne {\eta }_{i,j+1}\forall j$
C2	${\mu }_{i,j,1}\ne {\mu }_{i,j+1,1}\forall j$	${\mu }_{i,j,2}\ne {\mu }_{i,j+1,2}\forall j$	${\beta }_j\ne {\beta }_{j+1}\forall j$	${\eta }_{i,j}\ne {\eta }_{i,j+1}\forall j$

The block-orthogonal model specification given by (4) for each panel unit is: (7) $$y_i=D_{\iota }({\lambda }^0){\mu }_{i,1}^y+D_{\tau }({\lambda }^0){\mu }_{i,2}^y+x_i^l({\lambda }^0)\beta +F^l({\lambda }^0){\eta }_i^{y\prime }+{\xi }_i^y,$$(7) with $l\in \left\{A,B,C\right\}$ denoting the model and y_i and ξ_i being T-vectors of the dependent variable and disturbance term, respectively, for the ith panel unit. $D_{\iota }({\lambda }^0)={[D_{\iota ,1}\cdots D_{\iota ,J^0+1}]}_{T\times (J^0+1)}=\mathit{\text{diag}}({\iota }_1,\dots ,{\iota }_{J^0+1})$ and $D_{\tau }({\lambda }^0)={[D_{\tau ,1}\cdots D_{\tau ,J^0+1}]}_{T\times (J^0+1)}=\mathit{\text{diag}}({\tau }_1,\dots ,{\tau }_{J^0+1})$ are diagonal matrices defined with the elements between parentheses, ι_j is a $(T_j^0-T_{j-1}^0)$-vector of ones, ${\tau }_j={(1,2,\dots ,T_j^0-T_{j-1}^0)}^{\prime }$ is the time trend of the jth regime, and ${\mu }_{i,k}^y={({\mu }_{i,1,k}^y,\dots ,{\mu }_{i,J^0+1,k}^y)}^{\prime },k\in \left\{1,2\right\}$. For Model A, $x_i^A({\lambda }^0)=x_i={[x_{i,1}\dots x_{i,k}]}_{(T\times k)}$ and $\beta ={\beta }_0$ is a k-vector of parameters—subscript “0” denotes that the parameters are not affected by structural breaks. For Models B and C: (8) $$\begin{array}{c}{x}_i^l({\lambda }^0)={\left[x_i\left(1\right)x_i\left(2\right)\dots x_i\left(J^0+1\right)\right]}_{\left(T\times k\left(J^0+1\right)\right)};\\ \hspace{1em}l\in \left\{B,C\right\},\end{array}$$(8) with $x_i(j)=[(D_{\iota ,j}\odot x_{i,1})\dots (D_{\iota ,j}\odot x_{i,k})]$, $j=1,\dots ,J^0+1$, where $\odot$ is the element-wise (Hadamard) product and $\beta ={({\beta }_1^{\prime },\dots ,{\beta }_{J^0+1}^{\prime })}^{\prime }$. Similarly, $F^l({\lambda }^0)=F={[F_1\dots F_{m^0}]}_{(T\times m^0)},l\in \left\{A,B\right\}$, and ${\eta }_i^y={\eta }_{i,0}^y$ is a m⁰-row vector of loadings, whereas for Model C: (9) $$F^C({\lambda }^0)={\left[F\left(1\right)F\left(2\right)\dots F\left(J^0+1\right)\right]}_{\left(T\times m^0\left(J^0+1\right)\right)},$$(9) with $F(j)=[(D_{\iota ,j}\odot F_1)\dots (D_{\iota ,j}\odot F_{m^0})],j=1,\dots ,J^0+1$, and ${\eta }_i^y=({\eta }_{i,1}^y,\dots ,{\eta }_{i,J^0+1}^y)$ is the $m^0(J^0+1)$-row vector of factor loadings. The estimation of β in (7) can be performed using the panel data pooled estimator, although it would require both the common factors and the structural break dates to be known. Certainly, there are some cases where the common factors can be thought to be observable as discussed in Pesaran (2015), but this situation is rarely found in practice. The same applies to the structural break dates, so that the next section designs an estimation procedure that will allow the empirical implementation of our proposal when both the common factors and the structural breaks are unknown—the known common factors and/or structural breaks situation can be seen as particular cases.

3 Unknown Common Factors and Structural Breaks

The discussion that is conducted in this section focuses on Model C2, since the other specifications are obtained as particular cases. To ease the derivations, at this stage it is assumed that the number of structural breaks J⁰ (but not their position) is known, an assumption that will be relaxed below. We follow Pesaran (2006, 2007) and use cross-section averages to capture the unobserved common factors, which can be written in matrix notation from the block-orthogonal model specification as $\overline{Y}=D(\lambda ){\overline{\mu }}^{\prime }+F(\lambda ){\overline{\pi }}^{\prime }+\overline{U}$, where $D(\lambda )=\mathit{\text{diag}}(\left[{\iota }_1,{\tau }_1\right],\dots ,\left[{\iota }_{J^0+1},{\tau }_{J^0+1}\right]),F(\lambda )=F^C(\lambda ),\overline{\pi }={({\overline{\pi }}_1,\dots ,{\overline{\pi }}_{J^0+1})}^{\prime }$, ${\overline{Y}}_t=N^{-1}{\sum }_{i=1}^N{Y}_{i,t},\overline{\mu }=N^{-1}{\sum }_{i=1}^N{\mu }_i$, ${\overline{\pi }}_j=N^{-1}{\sum }_{i=1}^N{\pi }_{i,j}$ and ${\overline{U}}_t=N^{-1}{\sum }_{i=1}^N{U}_{i,t}$.

Assumption 6.

$({\overline{\pi }}_j)=m^0\le (1+k)$ for all N as $N\to \infty$ $\forall j,j=1,\dots ,J^0+1$.

If the rank condition established in Assumption 6 is met, we have: (10) $$F=\left(\overline{Y}-D(\lambda ){\overline{\pi }}^{\prime }-\overline{U}\right)\overline{\pi }{\left({\overline{\pi }}^{\prime }\overline{\pi }\right)}^{-1},$$(10) with ${\overline{U}}_t=O_p(T^{1/2}{N}^{-1/2})$ so that ${\overline{U}}_t\stackrel{p}{\to }0$ for (i) fixed T and $N\to \infty$, and (ii) $T\to \infty$ and $N\to \infty$ with $T/N\to 0$. Then, the observable cross-section averages ${\overline{h}}_t={(D_t,D_{\iota ,1,t}{\overline{Y}}_t^{\prime },\dots ,D_{\iota ,J^0+1,t}{\overline{Y}}_t^{\prime })}^{\prime }$ can be used to proxy the unobserved factors. Adding and subtracting (10) in (7) we have: $$y_i=D_{\iota }(\lambda ){\mu }_{i,1}^y+D_{\tau }(\lambda ){\mu }_{i,2}^y+x_i(\lambda )\beta +F(\lambda ){\pi }_i^{y\prime }$$ $$\pm \left(\overline{Y}-D(\lambda ){\overline{\pi }}^{\prime }-\overline{U}\right)\overline{\pi }{\left({\overline{\pi }}^{\prime }\overline{\pi }\right)}^{-1}{\pi }_i^{y\prime }+{\xi }_i^y$$ (11) $$=D_{\iota }(\lambda ){\mu }_{i,1}^y+D_{\tau }(\lambda ){\mu }_{i,2}^y+x_i(\lambda )\beta +\overline{Y}{\delta }_i^y+{\upsilon }_i,$$(11) with ${\delta }_i^y=\overline{\pi }{({\overline{\pi }}^{\prime }\overline{\pi })}^{-1}{\pi }_i^{y\prime }$ and ${\upsilon }_i={\xi }_i^y+(F(\lambda )-\overline{Y}\overline{\pi }{({\overline{\pi }}^{\prime }\overline{\pi })}^{-1}){\pi }_i^{y\prime }$. This defines an extension of the so-called cross-section augmented regression model in Holly, Pesaran, and Yamagata (2010) and Banerjee and Carrion-i-Silvestre (2017): (12) $$y_{i,t}={\mu }_{i,j,1}^y+{\mu }_{i,j,2}^y\left(t-T_{j-1}\right)+x_{i,t}^{\prime }{\beta }_j+{\overline{Y}}_t^{\prime }{\delta }_{i,j}^y+{\upsilon }_{i,t},$$(12)

$T_{j-1}<t\le T_j,j=1,\dots ,J^0+1$, where $T_B={(T_1,\dots ,T_{J^0})}^{\prime }$ is a vector of generic break dates. Note that it is implicitly assumed that $m^0=k+1$ since all k + 1 cross-section averages are included in (12). Let us define the projection matrix $\overline{M}(\lambda )=I-\overline{H}(\lambda ){(\overline{H}{(\lambda )}^{\prime }\overline{H}(\lambda ))}^{-1}\overline{H}{(\lambda )}^{\prime }$, with $\overline{H}(\lambda )=[D(\lambda )\overline{z}(\lambda )]$, $\overline{z}(\lambda )={[\overline{z}(1)\dots \overline{z}(J^0+1)]}_{(T\times (k+1)(J^0+1))}$, and $\overline{z}(j)=[(D_{\iota ,j}\odot {\overline{x}}_1)\dots (D_{\iota ,j}\odot {\overline{x}}_k)(D_{\iota ,j}\odot \overline{y})]$, $j=1,\dots ,J^0+1$. Then, the estimation of β in (11) can be done using the pooled CCE estimator (PCCE): (13) $${\widehat{\beta }}_{\mathit{\text{PCCE}}}(\lambda )={\left(\sum _{i=1}^N{x}_i{(\lambda )}^{\prime }\overline{M}(\lambda )x_i(\lambda )\right)}^{-1}\left(\sum _{i=1}^N{x}_i{(\lambda )}^{\prime }\overline{M}(\lambda )y_i\right).$$(13)

Following Bai (2010), Kim (2011, 2014), and Baltagi, Feng, and Kao (2016, 2019), the estimation of the break dates is carried out minimizing the global sum of squared residuals (SSR), so that: (14) $${\widehat{T}}_B=\mathrm{\text{arg}}\underset{\lambda \in \Lambda }{\mathrm{\text{min}}}\mathit{\text{SSR}}\left(\lambda \right),$$(14) where ${\widehat{T}}_B={({\widehat{T}}_1,\dots ,{\widehat{T}}_{J^0})}^{\prime }$, $\widehat{\lambda }={\widehat{T}}_B/T,\mathit{\text{SSR}}(\lambda )={\sum }_{i=1}^N{\mathit{\text{SSR}}}_i(\lambda ),{\mathit{\text{SSR}}}_i(\lambda )$ being the SSR for the ith panel unit associated with the estimation of (12). Thus, the estimation of the break dates can be obtained carrying out a grid search over all potential combinations of J⁰ structural breaks.¹ The following theorem analyzes the consistency of ${\widehat{\beta }}_{\mathit{\text{PCCE}}}(\widehat{\lambda })$ and $\widehat{\lambda }$. For completeness, we also include the known breaks case.

Theorem 1.

Let $Y_{i,t}$ be a vector of $(1+k)$ stochastic processes with DGP given by (1)–(3) and satisfying Assumptions 1–6. Then, as $(T,N)\to \infty$:

1. Known break dates case. For Models A, B, and C, ${\widehat{\beta }}_{\mathit{\text{PCCE}}}({\lambda }^0)\stackrel{p}{\to }\beta ={\Omega }_{U_x{U}_x}^{-1}{\Omega }_{U_x{U}_y}.$

2. Unknown break dates case:

   1. Model A1: $(\widehat{\lambda }-{\lambda }^0)$ does not converge to zero in probability.

   2. Models A2, B1, B2, C1, and C2: (i) $(\widehat{\lambda }-{\lambda }^0)\stackrel{p}{\to }0$, with $(\widehat{\lambda }-{\lambda }^0)=O_p(T^{-1/2})$ for Models A2, B2 and C2, $(\widehat{\lambda }-{\lambda }^0)=o_p(N^{-1/2})$ for Models B1 and C1, and (ii) ${\widehat{\beta }}_{\mathit{\text{PCCE}}}(\widehat{\lambda })\stackrel{p}{\to }\beta ={\Omega }_{U_x{U}_x}^{-1}{\Omega }_{U_x{U}_y}$.

See the companion appendix for the proof. Theorem 1 establishes that $\widehat{\lambda }$ is not consistent for Model A1 since the effects of the structural changes are dominated, in the limit, by the stochastic trend components,² whereas it is consistent for the other models. In addition, it is shown that the PCCE estimator converges toward the long-run average coefficient β in the known breaks case for all models, and in the unknown breaks case for the models for which a consistent break fraction vector estimate can be obtained. This establishes the basis for the derivation of our panel data cointegration test. Finally, it is worth noting that Theorem 1 establishes the consistency of the break fraction estimates for Models A2, B1, B2, C1, and C2. As for the break dates estimates, the results in Theorem 1 show that it is not possible to obtain consistent estimates of the break dates for Models A2, B2, and C2, consistency only applies to break fraction estimates. For Models B1 and C1, consistent break dates can be obtained if $T{N}^{-1/2}\to 0$ is assumed. Although there are some contributions in the panel literature that derive consistent break dates estimation procedures—see Bai (2010), Kim (2011, 2014) and Baltagi, Feng and Kao (2016, 2019), among others—these approaches assume that the disturbance term of the model is I(0), whereas here we deal with the novel panel data spurious regression with structural breaks case. Consequently, the framework under which those results are obtained is not comparable to our panel data spurious regression set-up.

It is possible to define an alternative estimation of the break dates that relies on the minimization of the SSR of the model: (15) $$y_{i,t}={\mu }_{i,j,1}^y+{\mu }_{i,j,2}^y\left(t-T_{j-1}\right)+x_{i,t}^{\prime }{\widehat{\beta }}_j+{\upsilon }_{i,t}^{\ast },$$(15) instead of using (12), where ${\widehat{\beta }}_j$ refers to the jth regime parameters in ${\widehat{\beta }}_{\mathit{\text{PCCE}}}(\lambda )$ for a generic $T_B$. The alternative estimation of the break dates is then defined as: (16) $${\tilde{T}}_B=\mathrm{\text{arg}}\underset{\lambda \in \Lambda }{\mathrm{\text{min}}}{\mathit{\text{SSR}}}^{\ast }\left(\lambda \right),$$(16) where ${\tilde{T}}_B={({\tilde{T}}_1,\dots ,{\tilde{T}}_{J^0})}^{\prime },\tilde{\lambda }={\tilde{T}}_B/T,{\mathit{\text{SSR}}}^{\ast }(\lambda )={\sum }_{i=1}^N{\mathit{\text{SSR}}}_i^{\ast }(\lambda ),{\mathit{\text{SSR}}}_i^{\ast }(\lambda )$ being the SSR for the ith panel unit associated with the estimation of (15). It is worth noting that this procedure ignores the common factors in (15), although they are considered in the estimation of ${\widehat{\beta }}_{\mathit{\text{PCCE}}}(\lambda )$. Therefore, this estimator can be seen as an hybrid method in which the common factors are used in the estimation of ${\widehat{\beta }}_{\mathit{\text{PCCE}}}(\lambda )$, but not in the minimization of the SSR of (15). The following theorem summarizes the main features of this alternative estimator.

Theorem 2.

Let $Y_{i,t}$ be a vector of $(1+k)$ stochastic processes with DGP given by (1)–(3) and satisfying Assumptions 1–6. Then, as $(T,N)\to \infty ,(\tilde{\lambda }-{\lambda }^0)=o_p(N^{-1/4}{T}^{-1/2})$ for Models B1 and C1, and $(\tilde{\lambda }-{\lambda }^0)=o_p(N^{-1/3}{T}^{-1})$ for Models A2, B2, and C2.

The proof is given in the companion appendix. As can be seen, the rate of convergence of $\tilde{\lambda }$ is faster than the $\widehat{\lambda }$ one for a given model, which suggests that better performance of the former in finite samples are to be expected. This result is due to the fact that the alternative estimator avoids projecting the observable variables in the model against the common factors, which leads us to obtain a stronger signal around the combination of break dates that minimizes the SSR of (15)—see the appendix. Finally, consistent break dates can be obtained for Models A2, B2, and C2—although at a slow rate of convergence—and for Models B1 and C1 if it is assumed that $T{N}^{-1/2}\to 0$. As mentioned above, the existing panel data structural break estimation techniques that have been recently proposed in the literature assume that the disturbance term of the model is I(0), which leads to consistent estimates of the break dates with higher rates of convergence.

The estimation of the number of structural breaks can be carried out using panel information criteria similar to the ones proposed in Bai and Ng (2002, 2004). To be specific, let us define an information criteria of the form: $$\mathit{\text{IC}}(\widehat{\lambda },J)=\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}{\widehat{\sigma }}^2(\widehat{\lambda },J)+\left(\left(1+J\right)k\right)g\left(N,T\right),$$with ${\widehat{\sigma }}^2(\widehat{\lambda },J)=N^{-1}{T}^{-1}{\sum }_{i=1}^N{\sum }_{t=1}^T\Delta {\widehat{\upsilon }}_{i,t}^2,J=0,1,\dots ,J_{\mathrm{\text{max}}}$, where ${\widehat{\upsilon }}_{i,t}$ denotes the estimated residuals from (12) and $g(N,T)$ is the penalty function—one possibility suggested in Bai and Ng (2002) is the panel BIC that is established with $g(N,T)=\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}(\mathit{\text{NT}}/(N+T))(N+T)/\mathit{\text{NT}}$. Then, $\widehat{J}=\mathrm{\text{arg}}{\mathrm{\text{min}}}_{J=0,\dots ,J_{\mathrm{\text{max}}}}\mathit{\text{IC}}(\widehat{\lambda },J)$. Similarly, we can define $\tilde{J}$ if $\tilde{\lambda }$ is used instead of $\widehat{\lambda }$ in the estimation of (12). The following theorem shows that the suggested information criterion provides consistent estimation of the number of structural breaks.

Theorem 3.

Let $Y_{i,t}$ be a vector of $(1+k)$ stochastic processes with DGP given by (1)–(3) and satisfying Assumptions 1–6. Then, $\underset{N,T\to \infty }{\mathrm{\text{lim}}}\mathrm{\text{Pr}}{(}_J^{°}=J^0)=1$, ${}_J^{°}\in \{\widehat{J},\tilde{J}\}$, if (i) $g(N,T)\to 0$ and (ii) $C_{\mathit{\text{NT}}}^{2}g(N,T)\to \infty$ as $(T,N)\to \infty$, where $C_{\mathit{\text{NT}}}=\mathrm{\text{min}}\{\sqrt{N},\sqrt{T}\}$.

The proof is outlined in the appendix. This result defines a framework general enough to be able to treat all the elements of the model in either an endogenous or exogenous way.

4 PCCE-based Cointegration Test Statistic

Following Pesaran (2007), Holly, Pesaran and Yamagata (2010), and Banerjee and Carrion-i-Silvestre (2017), we design a cointegration statistic that approximates the unobserved common factors using cross-section averages of the observable variables. The simplicity of this approach comes at the cost of restricting the framework defined above since Pesaran (2007) and Pesaran, Smith and Yamagata (2013) assume that both the common factor and the idiosyncratic components are I(1) under the null hypothesis of unit root—that is, $m^0=m_1^0$. In what follows, this section focuses on the most general specification given by Model C2 using $\widehat{\lambda }$, although the procedure is valid for the other model specifications and the $\tilde{\lambda }$ estimator.

Since ${\widehat{\beta }}_{\mathit{\text{PCCE}}}(\widehat{\lambda })$ is a consistent estimate of β—for all models except A1—we define ${\widehat{y}}_{i,t}=y_{i,t}-x_{i,t}^{\prime }{\widehat{\beta }}_{\mathit{\text{PCCE}}}(\widehat{\lambda })$ and base the testing of the null hypothesis of no cointegration for each unit on the cross-section augmented ADF (CADF) type regression equation:³ (17) $$\begin{array}{c}\Delta {\widehat{y}}_{i,t}=\sum _{j=0}^{\widehat{J}}{\theta }_{i,j}D{U}_{j,t}+\sum _{j=1}^{\widehat{J}}{\gamma }_{i,j}D{({\widehat{T}}_j)}_t+\sum _{j=0}^{\widehat{J}}{\vartheta }_{i,j}D{T}_{j,t}+{\alpha }_{i,0}{\widehat{y}}_{i,t-1}\\ +\sum _{l=1}^{p_i}{\alpha }_{i,l}\Delta {\widehat{y}}_{i,t-l}+\sum _{j=0}^{\widehat{J}}{\phi }_{i,j}^{\prime }{\left(D{U}_j\odot \overline{A}\right)}_{t-1}\\ +\sum _{j=0}^{\widehat{J}}\sum _{l=0}^{p_i}{\kappa }_{i,j,l}^{\prime }{\left(D{U}_j\odot \Delta \overline{A}\right)}_{t-l}+{\nu }_{i,t},\end{array}$$(17) with $D{({\widehat{T}}_j)}_t=1$ for $t={\widehat{T}}_j+1$, 0 otherwise, $j=1,\dots ,\widehat{J}$. Note that the order of augmentation p_i in (17) is heterogeneous, which can be selected using the modified information criteria in Ng and Perron (2001)—homogeneous p_i as in Pesaran (2007) can be imposed, if desired. The pseudo t-ratio statistic of ${\widehat{\alpha }}_{i,0}$ in (17), $t_{{\widehat{\alpha }}_{i,0}}(\widehat{\lambda })$, is used to test the null hypothesis of no cointegration.⁴ The number of common factors that is assumed defines ${\overline{A}}_t$, so that ${\overline{A}}_t={\overline{\widehat{y}}}_t$ for $m^0=1$, whereas ${\overline{A}}_t={({\overline{\widehat{y}}}_t,{\overline{x}}_{1,t},\dots ,{\overline{x}}_{k,t})}^{\prime }$ for $m^0=1+k$. For the intermediate cases where $m^0<1+k,{\overline{A}}_t$ will be defined with ${\overline{\widehat{y}}}_t$ and $m^0-1$ elements of ${\overline{x}}_t$. The following theorem provides the limiting distribution of $t_{{\widehat{\alpha }}_{i,0}}(\widehat{\lambda })$ for both the known and unknown breaks cases.

Theorem 4.

Let $Y_{i,t}$ be a vector of $(1+k)$ stochastic processes with DGP given by (1)–(3) and satisfying Assumptions 1–6. The $t_{{\widehat{\alpha }}_{i,0}}{(}_{\lambda }^{°})$ statistic, ${}_{\lambda }^{°}\in \{\widehat{\lambda },\tilde{\lambda }\}$, converges sequentially (if $N\to \infty$ first, then $T\to \infty$, that is, ${\Rightarrow }_{{(N,T)}_{\text{seq}}}$) and jointly (if $(N,T)\stackrel{j}{\to }\infty$ with $\sqrt{T}/N\to 0$, that is, ${\Rightarrow }_{{(N,T)}_j}$) to: $$t_{{\widehat{\alpha }}_{i,0}}{(}_{\lambda }^{°})\Rightarrow \frac{{\int }_0^1{W}_i\left(r\right)d{W}_i\left(r\right)-{\omega }_{\mathit{\text{iF}}}^l{\left({\lambda }^0\right)}^{\prime }{G}_F^l{\left({\lambda }^0\right)}^{-1}{\pi }_{\mathit{\text{iF}}}^l\left({\lambda }^0\right)}{{\left({\int }_0^1{W}_i^2\left(r\right)\mathit{\text{dr}}-{\pi }_{\mathit{\text{iF}}}^l{\left({\lambda }^0\right)}^{\prime }{G}_F^l{\left({\lambda }^0\right)}^{-1}{\pi }_{\mathit{\text{iF}}}^l\left({\lambda }^0\right)\right)}^{1/2}}$$ $$\equiv R\left({\lambda }^0,l\right),$$where $W_i(r)$ denotes a scalar standard Brownian motion, ${\omega }_{\mathit{\text{iF}}}^l({\lambda }^0)$, $G_F^l({\lambda }^0)$ and ${\pi }_{\mathit{\text{iF}}}^l({\lambda }^0)$ are functions of a m⁰-vector of standard Brownian motions associated with the common factors that are defined in the appendix, with $l\in \{A1,A2,B1,B2,C1,C2\}$ for the known breaks case $(\widehat{\lambda }=\tilde{\lambda }={\lambda }^0)$ and, for the unknown breaks case, $l\in \left\{B1,C1\right\}$ when $\widehat{\lambda }$ is used and $l\in \left\{A2,B1,B2,C1,C2\right\}$ when $\tilde{\lambda }$ is applied.

The proof is given in the appendix and some remarks are in order. First, the limiting distribution of $t_{{\widehat{\alpha }}_{i,0}}(\widehat{\lambda })$ and $t_{{\widehat{\alpha }}_{i,0}}(\tilde{\lambda })$ depends on the number and position of structural breaks $({\lambda }^0)$ and on the number of I(1) nonstationary common factors—note that it is assumed that $m^0=m_1^0$. If analysts have knowledge—for example, based on economic theory—about the number of common factors, one can consider such information when ${\overline{A}}_t$ is defined. In this regard, we could follow the strategy in Pesaran et al. (2013) and compute the statistic using all possible combinations of $m^0$ cross-section averages available in the system as a way of obtaining robust conclusions. When m⁰ is unknown, we can follow a conservative strategy and assume that the rank condition is satisfied with equality. The price that we would pay if $m^0<1+k$, but we impose $m^0=1+k$, is to have a test statistic with empirical size smaller than the nominal size accompanied by loss of power. The advantage is to allow us to remain agnostic about the number of integrated stochastic trends driving the data.⁵ Second, the appendix shows that the limiting distribution of $t_{{\widehat{\alpha }}_{i,0}}({\lambda }^0)$ for Models A and B is the same, for each variant of the deterministic component, whereas it is different for Model C.

Third, for the unknown breaks case and Models B1 and C1, $t_{{\widehat{\alpha }}_{i,0}}(\widehat{\lambda })$ converges to the same limiting distribution as for the known structural breaks case. Unfortunately, this does not occur for Models A2, B2, and C2. In these cases, we suggest the implementation of the modified procedure in Kim and Perron (2009) that relies on the use of trimmed data. In brief, the idea is based on establishing a window of observations around a consistent estimate of λ—that is, $(\widehat{\lambda }-{\lambda }^0)=O_p(T^{-a})$ for some $0<a<1$—of length $2\omega (T)$ with $\omega (T)\equiv \kappa T^{\delta }$, $\kappa >0$ and $-1<-a<\delta <0$. The window defines the set of observations compressed between ${\widehat{T}}_l+1\equiv T(\widehat{\lambda }-\omega (T))+1$ and ${\widehat{T}}_h\equiv T(\widehat{\lambda }+\omega (T))$ and has the characteristics of: (i) increasing slowly enough to be asymptotically negligible relative to T and, (ii) increasing fast enough to include the true structural break dates since ${\widehat{T}}_l-{\widehat{T}}_B=T(\widehat{\lambda }-\omega (T))-T{\lambda }^0=(T^{-a-\delta }{T}^a(\widehat{\lambda }-{\lambda }^0)-\kappa )T^{\delta +1}\stackrel{p}{\to }-\infty$ and ${\widehat{T}}_B-{\widehat{T}}_h=T{\lambda }^0-T(\widehat{\lambda }+\omega (T))=(T^{-a-\delta }{T}^a({\lambda }^0-\widehat{\lambda })-\kappa )T^{\delta +1}\stackrel{p}{\to }-\infty$. Once ${\widehat{T}}_l$ and ${\widehat{T}}_h$ have been specified, Kim and Perron (2009) define a new dataset as:⁶ $${\widehat{y}}_{i,t}^n={\widehat{y}}_{i,t+{\widehat{T}}_{j-1,h}-{\widehat{T}}_{j-1,l}}-S({\widehat{y}}_{i,t},j-1)\text{ for }{\widehat{T}}_{j-1}<t\le {\widehat{T}}_j,$$with $S({\widehat{y}}_{i,t},j-1)={\widehat{y}}_{i,{\widehat{T}}_{j-1,h}}-{\widehat{y}}_{i,{\widehat{T}}_{j-1,l}},{\widehat{T}}_{0,h}={\widehat{T}}_{0,l}=0$ and $S({\widehat{y}}_{i,t},0)=0$—the same transformation applies to obtain ${\overline{A}}_t^n$ with $S({\overline{A}}_t,j-1)={\overline{A}}_{{\widehat{T}}_{j-1,h}}-{\overline{A}}_{{\widehat{T}}_{j-1,l}}$. The new dataset removes the data points between ${\widehat{T}}_l+1$ and ${\widehat{T}}_h$ and reconnects the remaining data shifting down the data after the window with the $S(·,j-1)$ function. Then, the pseudo t-ratio statistic of ${\widehat{\alpha }}_{i,0}$ in (17) is computed using the trimmed data ${\widehat{y}}_{i,t}^n$ and ${\overline{A}}_t^n$ with the break dates given by ${\widehat{T}}_l$. The resulting statistic is denoted as $t_{{\widehat{\alpha }}_{i,0}}({\widehat{\lambda }}_{\mathit{\text{tr}}})$, which limiting distribution is given in the following corollary.

Corollary 1.

Let $Y_{i,t}$ be a vector of $(1+k)$ stochastic processes with DGP given by (1)–(3) and satisfying Assumptions 1–6. For the unknown structural breaks case and as either ${(N,T)}_{\text{seq}}$ or ${(N,T)}_j,t_{{\widehat{\alpha }}_{i,0}}({\widehat{\lambda }}_{\mathit{\text{tr}}})\Rightarrow R({\lambda }^0,l),l\in \{A2,B2,C2\}$.

The proof derives from the developments in Theorem 4. As can be seen, the limiting distribution of $t_{{\widehat{\alpha }}_{i,0}}({\widehat{\lambda }}_{\mathit{\text{tr}}})$ equals the limiting distribution that is obtained for the known structural breaks case. Finally, it should be stressed that the limiting distributions in Theorem 4 and Corollary 1 also correspond to the distributions obtained for a panel unit root test with multiple structural breaks that affect the level and/or the slope of the time trend. Therefore, as side contributions, this article (i) extends the panel data unit root tests in Pesaran (2007) and Pesaran, Smith and Yamagata (2013) to the multiple structural breaks case—affecting the deterministic component and/or the factor loadings—and (ii) provides an extension of Perron (1989, 1990) and Kim and Perron (2009) unit root tests to panel data framework.

The limiting distributions in Theorem 4 are approximated by Monte Carlo simulation. A limited set of critical values for $t_{{\widehat{\alpha }}_{i,0}}({\lambda }^0)$ are reported in Tables B.1 to B.6 for Models A1, B1, and C1, and in Tables B.7 to B.12 for Models A2, B2, and C2, with J = 1 and $m^0\le k+1\in \{2,3,4\}$ common factors. Pesaran (2007) and Pesaran, Smith and Yamagata (2013) also propose a truncated version of the $t_{{\widehat{\alpha }}_{i,0}}({\lambda }^0)$ statistic in order to ensure that the statistic has finite moments, although these papers provide evidence that for T > 15 the empirical distributions of the truncated and untruncated statistics are equivalent. Therefore, in this article we do not consider the truncated version of $t_{{\widehat{\alpha }}_{i,0}}({\lambda }^0)$.

The combination of the individual test statistics defines the cross-section augmented ADF panel cointegration statistics $\mathit{\text{CIPS}}{(}_{\lambda }^{°})=N^{-1}{\sum }_{i=1}^N{t}_{{\widehat{\alpha }}_{i,0}}{(}_{\lambda }^{°})$, ${}_{\lambda }^{°}\in \{\widehat{\lambda },\tilde{\lambda },{\widehat{\lambda }}_{\mathit{\text{tr}}}\}$, as proposed in Pesaran (2007). This statistic allows us to test the null hypothesis of no panel data cointegration against the alternative that there is a fraction of panel units for which cointegration holds. The critical values for the $\mathit{\text{CIPS}}({\lambda }^0)$ statistic are presented in Tables B.13 to B.18 for Models A1, B1, and C1, and in Tables B.19–B.24 for Models A2, B2, and C2, with J = 1 and $m^0\le k+1\in \{2,3,4\}$. Interestingly, the critical values seem to be symmetrical around ${\lambda }^0=0.5$ for large T, a characteristic that is more evident for the panel data cointegration statistic critical values. This feature has also been found for the ADF unit root test with one structural break proposed in Perron (1989, 1990).⁷

Westerlund, Hosseinkouchack and Solberger (2016) derive the local-to-unity asymptotic power functions of Pesaran (2007) type-test statistics, and show that these statistics have non-negligible power in a neighborhood of $T^{-1}$. This is a distinctive feature of these statistics when compared to other proposals in the literature, that have been shown to have non-negligible power in the neighborhood given by $N^{-\kappa }{T}^{-1},\kappa >0$—see, for instance, Moon and Perron (2008). This implies that as N increases Pesaran (2007) type-test statistics will tend to be dominated by other panel data tests having non-negligible local power for $\kappa >0$.

5 Finite Sample Performance

Let us consider the DGP defined by: (18) $$\begin{array}{c}{y}_{i,t}={\theta }_{i,1}D{U}_t+{\theta }_{i,2}D{T}_t+x_{i,t}^{\prime }{\varsigma }_0+D{U}_t{x}_{i,t}^{\prime }{\varsigma }_1+F_t^{\prime }{\phi }_{i,0}^y\\ +D{U}_t{F}_t^{\prime }{\phi }_{i,1}^y+u_{i,t}\\ \Delta x_{i,t}=\Delta F_t^{\prime }{\phi }_{i,0}^x+\Delta \left(D{U}_t{F}_t^{\prime }\right){\phi }_{i,1}^x+{\upsilon }_{i,t};\\ F_{j,t}=\rho F_{j,t-1}+w_{j,t};\\ u_{i,t}={\varphi }_i{u}_{i,t-1}+{\epsilon }_{i,t},\end{array}$$(18) where ${\theta }_{i,1}\sim N(10,1),{\theta }_{i,2}\sim N(0.3,1),{\varsigma }_0=1,{\varsigma }_1\in \{1,5\},{\phi }_{i,0}^y\sim U\left[0,1\right],{\phi }_{i,1}^y\sim U\left[0,2\right],{\phi }_{i,0}^x\sim U\left[0,1\right],{\phi }_{i,1}^x\sim U\left[0,3\right],{\upsilon }_{i,t}\sim N(0,1),w_{j,t}\sim N(0,1),j\in \left\{1,2\right\}$, and ${\epsilon }_{i,t}\sim N(0,1)$ are mutually independent groups. Under the null hypothesis of no cointegration ${\varphi }_i=1\forall i$, whereas under the alternative hypothesis of cointegration we set ${\varphi }_i=0.9\forall i$. This section investigates the performance of the proposed statistics considering one unknown structural break—additional simulation results for known structural breaks are available upon request. We distinguish three different cases of interest. Case 1 specifies a DGP in which both the deterministic component and the cointegration vector change, but the loadings remain constant across regimes—that is, ${\varsigma }_1\in \{1,5\}$ and ${\phi }_{i,1}^y={\phi }_{i,1}^x=0$ $\forall i$ in (18). Case 2 uses the Case 1 DGP but the computations are carried out allowing for a nonexistent change in factor loadings. This is done to proxy situations where the investigator may be uncertain about whether there are changes in the factor loadings and wishes to protect herself (in terms of size of the test) against this possibility. Finally, Case 3 deals with the pure structural break case where all parameters (deterministic component, cointegration vector and factor loadings) change across regimes. Note that for each of these cases we consider the two variants of the deterministic component—that is, the constant and the time trend cases. The break date is specified at ${\lambda }^0=0.5$ and we consider up to two common factors, $m^0\in \left\{1,2\right\}$, with $\rho \in \left\{1,0.99,0.95,0.9\right\}$. The time dimension is set at $T\in \left\{50,100,200\right\}$ and the cross-section dimension is $N\in \left\{20,50\right\}$. The nominal size is set at 5% and the critical values tabulated in the previous section are used.

5.1 Constant Case

5.1.1 Performance of the Break Fraction and Pooled Estimators

This section investigates the performance of $\widehat{\lambda }$ and $\tilde{\lambda }$ for Models B and C with the deterministic component given by the constant case—that is, ${\theta }_{i,2}=0\forall i$ in (18). Let us first focus on the $\widehat{\lambda }$ estimator. Figures C.1–C.4 present histograms of $\widehat{\lambda }$ for Model B1 (Case 1) when ${\varsigma }_1=1$ for all possible combinations that consider that the true number of common factors is $m^0\in \left\{1,2\right\}$ and the imposed number of common factors is $m\in \left\{1,2\right\}$. It is worth stressing that consistency of the $\widehat{\lambda }$ and $\tilde{\lambda }$ estimators has been established under the null hypothesis of spurious regression situation, but for completeness we also include the histograms for these estimators under cointegration.

In general and regardless of the order of integration of the common factors, histograms are concentrated around ${\lambda }^0$ with higher probability as either T or N increases. It is interesting to note that allowing for $m>m^0$ does not affect the performance of $\widehat{\lambda }$, whereas the probability mass of $\widehat{\lambda }$ around ${\lambda }^0$ increases with $m^0$. As expected, as the magnitude of the structural change that affect the slope parameters increases, the performance of $\widehat{\lambda }$ shows a substantial improvement, with almost all probability mass located on ${\lambda }^0$, regardless of m⁰, m and the order of integration of the idiosyncratic and the common factor components—see Figures C.5–C.8 with the histograms of $\widehat{\lambda }$ for Model B1 when ${\varsigma }_1=5$. These features are also found when the results for Cases 2 and 3 are analyzed—see Figures C.9–C.24. The use of the $\tilde{\lambda }$ estimator produces better results since now, and even for ${\varsigma }_1=1$, the probability mass of $\tilde{\lambda }$ around ${\lambda }^0$ is higher than the $\widehat{\lambda }$ one—see Figures C.25–C.48 that depict the histograms of $\tilde{\lambda }$ for Cases 1 to 3, $m\in \left\{1,2\right\},m^0\in \left\{1,2\right\}$ and ${\varsigma }_1\in \{1,5\}$. This simulation evidence is demonstration of the statement made in Theorems 1 and 2.

Finally, we have also analyzed the performance of ${\widehat{\beta }}_{\mathit{\text{PCCE}}}(\widehat{\lambda })$ and ${\widehat{\beta }}_{\mathit{\text{PCCE}}}(\tilde{\lambda })$ and found evidence that supports that its application leads to a consistent estimation of β. To be specific, the mean squared error (MSE) decreases as N and/or T increase, and as ρ moves away from one—detailed results are available upon request.

5.1.2 Empirical Size and Power of the Panel Cointegration Statistic

Table C.1 summarizes the performance of the $\mathit{\text{CIPS}}(\widehat{\lambda })$ statistic. When the magnitude of the structural change affecting the slope is ${\varsigma }_1=1$, we can observe size distortion problems for the panel cointegration statistic for Case 1 when the order of integration of both the idiosyncratic and the common factor components is the same, regardless of the (true and assumed) number of common factors. These distortions reduce as the model specification becomes more flexible—that is, Cases 2 and 3—and even disappear in some situations—see the results for $m^0=m=1$. As the magnitude of the structural break increases to ${\varsigma }_1=5$, the empirical size tends to the nominal one in all cases. It is worth highlighting that the misspecification error of either under-specifying ($0<m<m^0$) or over-specifying ($m^0<m$) the number of common factors does not cause high size distortion problems. The panel cointegration $\mathit{\text{CIPS}}(\widehat{\lambda })$ statistic becomes conservative under the null hypothesis of spurious regression if the common factor is I(0). This is a consequence of the violation of the key assumption underlying the design of the statistic—see Section 4. The empirical power of $\mathit{\text{CIPS}}(\widehat{\lambda })$ increases with N and T, and as ρ moves away from one. As expected, the more flexible model specifications—either because the type of model is more general (Case 2) or because the number of common factors is over-specified ($m^0<m$)—show lower empirical power values.

The performance of the statistical inference improves when $\tilde{\lambda }$ is used. Table C.2 evidences a clear enhancement of the empirical size figures when ${\varsigma }_1=1$—the exception is found for Case 3 with $m^0=m=2$ and T = 200, although the size distortion almost disappears for ${\varsigma }_1=5$. However and contrary to what has been found for $\mathit{\text{CIPS}}(\widehat{\lambda })$, the under-specification of the number of common factors causes mild over-size distortions for Model C—see the results for Case 3. This results is somehow to be expected since we are missing the presence of common factors in the model specification. As above, $\mathit{\text{CIPS}}(\tilde{\lambda })$ becomes conservative if the assumption of common order of integration of the idiosyncratic and common components under the null hypothesis of spurious regression is not satisfied. The $\mathit{\text{CIPS}}(\tilde{\lambda })$ statistic encompasses $\mathit{\text{CIPS}}(\widehat{\lambda })$ in term of empirical power, since the former statistic shows similar, if not higher, empirical power values—note that the higher power shown by $\mathit{\text{CIPS}}(\widehat{\lambda })$ when ${\varsigma }_1=1$ are due to the mild size distortions that have been documented.

Table 2 Parameter estimates and CIPS panel cointegration test statistic.

					Panel	CIPS
J	Model	${\widehat{T}}_B$	${\tilde{T}}_B$	$\widehat{\beta }$	BIC	m = 1	m = 2
0				1.23	–6.28	–2.76**	–2.82**
1	B1	1984		${(0.78,0.87)}^{\prime }$	–6.00	–2.87**	–3.39**
1	C1	1991		${(0.62,0.69)}^{\prime }$	–6.14	–1.79	–2.99**
1	B1		1985	${(0.91,0.89)}^{\prime }$	–5.98	–2.50**	–3.09**
1	C1		1986	${(0.27,0.73)}^{\prime }$	–6.29	–2.26	–2.87**
2	B1	${(1984,1993)}^{\prime }$		${(0.77,0.94,0.90)}^{\prime }$	–5.81	–2.58	–2.71
2	C1	${(1988,2004)}^{\prime }$		${(0.43,0.82,0.45)}^{\prime }$	–6.07	–3.15**	–1.99
2	B1		${(1984,2009)}^{\prime }$	${(0.86,0.91,0.77)}^{\prime }$	–5.76	–3.58**	–3.20**
2	C1		${(1986,2003)}^{\prime }$	${(0.27,0.70,0.54)}^{\prime }$	–6.04	–3.12**	–3.02

** and * denote rejection of the null hypothesis at the 5% and 10% significance levels, respectively

To the best of our knowledge, there are no other proposals in the literature that can be used to establish a direct comparison with the test statistics that are proposed in this article. The closest proposal might the PANIC-based panel cointegration test statistic that is developed by Banerjee and Carrion-i-Silvestre (2015). That paper designed a panel data statistic to test the joint null hypothesis of spurious regression without structural breaks against the alternative hypothesis of panel cointegration with structural breaks for large T compared to N. As mentioned in the introduction, the cross-section dependence was introduced in the specification using an approximate common factor model. Tables C.3 and C.4 present the empirical size and power for the panel ADF statistic that is computed for the estimated idiosyncratic component ($Z_c^{\ast }$) and the $M{Q}_c^c$ statistic to estimate the number of I(0) and I(1) common stochastic trends. The simulation results are reported for Cases 1 and 3, considering $m^0\in \left\{1,2\right\}$ and allowing for up to $m_{\mathrm{\text{max}}}=6$ common factors—the number of common factors is estimated using the panel BIC statistic in Bai and Ng (2004). The $Z_c^{\ast }$ statistic shows important size distortions, the higher the value of ${\varsigma }_1$, which in fact can be interpreted in terms of empirical power since this situation is a mere consequence of the violation of the DGP that is assumed under null hypothesis that is assumed in this proposal. Thus, note that the DGP defined in (18) allows for the presence of structural breaks under the null hypothesis, so that the joint hypothesis of spurious regression without structural breaks that is assumed by Banerjee and Carrion-i-Silvestre (2015) is not fulfilled. The behavior of the $M{Q}_c^c$ statistic is also affected by this feature, although it is fair to note that the $M{Q}_c^c$ statistic tries to capture that misspecification error detecting more common factors than really exist—in general, it can be shown that the estimated number of common factors is double the true number.

Finally and following the suggestion in Westerlund (2018), we have also studied the performance of the CIPS panel cointegration test statistic when the structural breaks are ignored. Table C.5 summarizes the empirical size and power of the CIPS proposed in Banerjee and Carrion-i-Silvestre (2017) for Cases 1 and 3, and ${\varsigma }_1\in \{1,5\}$. As can be seen, the CIPS presents size distortions, the higher the value of ${\varsigma }_1$. For ${\varsigma }_1=1$ and $m^0=1$, the size distortions can be reduced if we over-specify the number of common factors, although this desirable feature disappears as ${\varsigma }_1$ increases. The violation of the assumption of common order of integration of the idiosyncratic and common components under the null hypothesis of spurious regression leads to a conservative test statistic, but only when ${\varsigma }_1=1$ and $m>m^0$. For the other cases, the statistic shows over-rejection distortions. These results are in line with the literature that shows that the null hypothesis of panel unit root is rejected when cross-section dependence is ignored—see Banerjee, Marcellino and Osbat (2005). Here the cross-section dependence is considered, but to some extent the omission of common structural breaks might leave one source of cross-section dependence that links the panel units uncaptured. This shows that modelling the presence of structural breaks is required if misleading statistical inference is to be avoided.

5.1.3 Estimation of the Number of Structural Breaks

Table C.10 provides the frequency of the estimated number of structural breaks that derives from the use of the panel BIC statistic, considering up to two structural breaks. The ability of the panel BIC to select the correct number of structural breaks improves as N, T and ${\varsigma }_1$ increase, with the statistic that is based on $\tilde{\lambda }$ the one that produces better results.

Let us first focus on the results for ${\varsigma }_1=1$. Under the spurious regression case with ρ = 1, $J(\widehat{\lambda })$ tends to under-estimate $J$ for small N and T, with correct number of structural breaks detection frequencies that range between 0.26 and 0.48 for Model B1, and between 0.48 and 0.6 for Model C1. As can be seen, $J(\tilde{\lambda })$ outperforms $J(\widehat{\lambda })$ with correct frequencies in the range 0.52-0.86 for Model B1, and 0.57–0.7 for Model C1. Further, note that part of the under-estimation frequency of $J(\widehat{\lambda })$ moves to over-estimation as N increases for a given T—this feature is hardly observed for $J(\tilde{\lambda })$. When $\rho =0.9,J(\widehat{\lambda })$ still shows a tendency to under-estimate J for Model B1, although this behavior is attenuated for Model C1. As for the $J(\tilde{\lambda })$ estimator, the probability of correct detection remains unchanged as N increases, regardless of T, with a tendency of reallocate part of the under-estimation frequency on the over-estimation one as T increases. Under the cointegration scenario, the behavior of both $J(\tilde{\lambda })$ and $J(\widehat{\lambda })$ statistics improve with correct J estimation frequencies that tend to one as T increases, and with $J(\tilde{\lambda })$ outperforming $J(\widehat{\lambda })$ in all cases. Finally, when the magnitude of the structural break increases to ${\varsigma }_1=5$ both estimators provide good results.

Table C.12 summarizes the performance of the panel BIC when there are no structural breaks ($J^0=0$) affecting the parameters of the model. The frequency of correct selection of the number of breaks for both $J(\tilde{\lambda })$ and $J(\widehat{\lambda })$ statistics is similar, and higher than the one obtained for the one structural break case. In general, higher correct selection frequencies are obtained for Model C1. This results indicates that the proposed panel BIC is useful to detect the presence of structural breaks in our model set-up.

5.2 Time Trend Case

5.2.1 Performance of the Break Fraction and Pooled Estimators

This section investigates the performance of $\widehat{\lambda }$ and $\tilde{\lambda }$ for Models B and C with the deterministic component given by the time trend case—that is, ${\theta }_{i,1}\ne {\theta }_{i,2}\ne 0\forall i$ in (18). Figures C.49–C.52 depict the histograms of $\widehat{\lambda }$ for Model B2 (Case 1) when ${\varsigma }_1=1$ for all possible combinations that consider that the true number of common factors is $m^0\in \left\{1,2\right\}$ and the imposed number of common factors is $m\in \left\{1,2\right\}$. In general, we can observe that the probability mass concentrates around ${\lambda }^0$ as T and/or N increase for a given set of ${\varphi }_i$ and ρ parameters. The probability around ${\lambda }^0$ is slightly reduced for Cases 2 and 3, especially for small T, although it tends toward one as T and N get large—see Figures C.53–C.60. Note that when this estimation bias appears, the estimates are located around ${\lambda }^0$. As for the constant case, the use of $\tilde{\lambda }$ produces better results, since the probability mass of $\tilde{\lambda }$ around ${\lambda }^0$ is, in general, similar or larger than the $\widehat{\lambda }$ one—see Figures C.61–C.72.

The simulation experiment has also studied the performance of ${\widehat{\beta }}_{\mathit{\text{PCCE}}}(\widehat{\lambda })$ and has found that the MSE decreases as N and/or $T$ increase, and as ρ moves away from one—detailed results are available upon request. All these results reinforce the theoretical analysis that has been detailed in Theorem 1.

5.2.2 Empirical Size and Power of the Panel Cointegration Statistic

Table C.6 collects the empirical size and power of $\mathit{\text{CIPS}}({\widehat{\lambda }}_{\mathit{\text{tr}}})$. The statistic shows good performance, with an empirical size that is close to 0.05 when ${\varphi }_i=\rho =1$. Under-specification of the number of common factors does not seem to cause major size distortions for Case 1, although this is not the case for Case 3 and, to a lesser extent, for Case 2, for large N and T. As for the constant case, the statistic becomes conservative if the assumption of common order of integration of the idiosyncratic and common factor components is not met. The empirical power increases as T and N increase, and decreases as the model specification gets more complicated—this is something to be expected due to the presence of a higher number of parameters to be estimated. For completeness, we have also included the results that are based on the badly-behaved $\mathit{\text{CIPS}}(\widehat{\lambda })$ statistic, which behaves similar to $\mathit{\text{CIPS}}({\widehat{\lambda }}_{\mathit{\text{tr}}})$. This result might be surprising, but we need to take into account that the term that invalidates the use of $\mathit{\text{CIPS}}(\widehat{\lambda })$ in the limit also depends on the magnitude of the change in the slope of the time trend, so that small values of $\left||{\theta }_{i,2}\right||$ might have little effect on the limiting distribution.

We have also analyzed the performance of the statistics under local-to-zero breaks through the definition of shrinking break magnitudes for all parameters related to the structural instability. To be specific, ${\theta }_{i,1},{\theta }_{i,2},{\varsigma }_1,{\phi }_{i,1}^y$, and ${\phi }_{i,1}^x$ defined above have been rescaled by $T^{-1/2-\delta },\delta >0$—we have used an arbitrarily small positive value of d = 0.01—so that in the limit the effects of the structural break tend to zero.⁸ Table C.8 indicates that the performance of the $\mathit{\text{CIPS}}({\widehat{\lambda }}_{\mathit{\text{tr}}})$ statistic under shrinking breaks is similar to the fixed breaks one—the empirical size is close to the nominal one, whereas in some cases a mild drop in the empirical power is observed, although the empirical power equals one as T gets large.

The finite sample performance of $\mathit{\text{CIPS}}(\tilde{\lambda })$ is investigated in Table C.7, which reveals a test statistic with an empirical size that is close to 0.05 in most cases—with the exception of a mild over-rejection distortion found for Case 3 when $m^0=m=2$ with large T. The values for the empirical power of $\mathit{\text{CIPS}}(\tilde{\lambda })$ are higher than the $\mathit{\text{CIPS}}({\widehat{\lambda }}_{\mathit{\text{tr}}})$ ones—and similar to the $\mathit{\text{CIPS}}(\widehat{\lambda })$ ones—for a given set of ρ, T, and N values. As for the shrinking breaks configuration, we observe that the empirical size of $\mathit{\text{CIPS}}(\tilde{\lambda })$ is close to 0.05, although the empirical power is reduced when compared to the results that are based on fixed break magnitudes. To some extent, this is something to be expected since in this case, although $\tilde{\lambda }$ is consistent, the rate at which $\tilde{\lambda }$ tends to ${\lambda }^0$ is reduced—it has been shown that $(\tilde{\lambda }-{\lambda }^0)=o_p(N^{-1/3}{T}^{-1})$ for the fixed breaks case. From this point of view, it seems that $\mathit{\text{CIPS}}({\widehat{\lambda }}_{\mathit{\text{tr}}})$ outperforms $\mathit{\text{CIPS}}(\tilde{\lambda })$, although the empirical power equals one as T gets large.

Contrary to what has been done for the constant case, here we cannot compare the performance of the statistic with other existing proposal in the literature, since the case of panel cointegration testing with unknown structural changes that affect the slope of the time trend have not been previously addressed in the literature—Banerjee and Carrion-i-Silvestre (2015) only deal with the unknown breaks date situation for the constant case. However, we have investigated the properties of the statistical inference when the structural breaks are ignored with the computation of the panel cointegration statistic in Banerjee and Carrion-i-Silvestre (2017), which results are summarized in Table C.9. Unfortunately, in this case the empirical size is not controlled since, for instance, for Case 1 with $m^0=m=1$, the empirical size tends toward one as T and N increase. The size distortion is less noticeable for Case 3 with $m^0=m=1$, although it is still important in some cases—see results for $m^0=m=2$ and T = 200. In addition, the violation of the assumption of common order of integration of the idiosyncratic and common components under the null hypothesis of spurious regression does not seem to lead to a conservative test statistic. Finally, even in those cases where the size distortions are not very large—that is, Case 1 with $m^0=m=1$, N = 20 and T = 50—the empirical power of the statistic is only slightly above the empirical size—that is, 0.07 versus 0.09. Again, this suggests that empirical analyses should include parameter instabilities in the model specification if there is evidence that structural changes might have affected the model.

5.2.3 Estimation of the Number of Structural Breaks

Table C.11 provides the frequency of the estimated number of structural breaks that derives from the use of the panel BIC statistic, considering up to two structural breaks. The ability of the panel BIC to select the correct number of structural breaks is very good, with a correct detection frequency that range between 0.92 and 1 regardless of N, T, ${\varphi }_i$, and ρ, for both $\widehat{\lambda }$ and $\tilde{\lambda }$ based estimators, and model specifications. Table C.12 presents the frequency of the number of structural breaks estimation for model specifications that do not include structural breaks ($J^0=0$). For Model B2 the frequency of correct classification lies in the ranges [0.73, 0.92] for $J(\widehat{\lambda })$ and [0.69, 0.92] for $J(\tilde{\lambda })$ under the null hypothesis of spurious regression. For Model C2 the classification ranges narrow to [0.93, 0.96] for $J(\widehat{\lambda })$ and [0.90, 0.97] for $J(\tilde{\lambda })$ under the null hypothesis of spurious regression. Under the alternative hypothesis of cointegration, the frequency of correct estimation of the number of structural breaks equals one for both $J(\widehat{\lambda })$ and $J(\tilde{\lambda })$ statistics. This leads us to suggest the use of the panel BIC statistic in empirical applications.

6 Empirical Illustration

A well-functioning housing market has been shown to be very relevant for the proper evolution of credit and financial markets, which in turn has effects on macroeconomic variables such as output, fiscal deficit and unemployment. The empirical evidence that analyses the potential relationship between housing prices and real disposable income per capita is mixed, and depends both on the scope (national or regional) and on the period of analysis. Following Holly, Pesaran and Yamagata (2010), we focus on the US economy considering the 48 contiguous U.S. States and the District of Columbia (N = 49) using annual data between 1975 to 2019 (T = 45)—see Holly, Pesaran and Yamagata (2010) for the sources of the statistical information—with the model: (19) $$h{p}_{i,t}={\alpha }_i+\beta y_{i,t}+u_{i,t},$$(19) where $h{p}_{i,t}$ denotes the logarithm of the real housing prices index and $y_{i,t}$ is the logarithm of the real disposable income per capita. They argue that the use of panel data cointegration analysis can provide better statistical inference given the short time period of the available information. The preliminary analysis conducted in their paper reveals that cross-section dependence is present in the dataset.

As discussed in Holly, Pesaran and Yamagata (2010), the boom house prices started in early 2000 in the United States and accelerated during 2003–2006, something that has been interpreted in the literature as a housing price bubble. Figures C.73 and C.74 depict the variables, which show hump-shaped behavior in some of the house price time series during mid-eighties and in the middle of the first decade of 2000. If this is the case, it is safe to assume that the potential long-run relationship between housing prices and disposable income per capita might have been affected. This in turn implies that parameter instabilities should be accounted for in the model. The empirical specification that is used here generalizes the one in (19) by considering the effects of structural breaks as follows: (20) $$h{p}_{i,t}={\mu }_{i,j,1}+{\beta }_j{y}_{i,t}+{\overline{Y}}_t^{\prime }{\delta }_{i,j}^y+{\upsilon }_{i,t}\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{T}_{j-1}<t\le T_j,$$(20) with ${\overline{Y}}_t={({\overline{\mathit{\text{hp}}}}_t,{\overline{y}}_t)}^{\prime },{\overline{\mathit{\text{hp}}}}_t=N^{-1}{\sum }_{i=1}^{N}h{p}_{i,t},{\overline{y}}_t=N^{-1}{\sum }_{i=1}^N{y}_{i,t}$ and $j=1,\dots ,J$.

Table 2 summarizes the results that are based on both break dates estimation procedures discussed in Section 3, although it has to be born in mind that ${\tilde{T}}_B$ clearly outperforms ${\widehat{T}}_B$, so that more weight should be imposed on the conclusions drawn from the former estimator. The order of augmentation p_i in (17) is selected using the modified Akaike’s information criterion (MAIC) in Ng and Perron (2001) and Perron and Qu (2007) for each panel unit with up to five lags, and the trimming is set at $ϵ=0.2$. Let us first focus on the specification that considers one structural break. As can be seen, there is strong evidence of panel cointegration for Model B1, since the null hypothesis of spurious regression is rejected at the 5% significance level regardless of m. Evidence of cointegration is also found for Model C1, although only when m = 2. This situation illustrates the potential effect on the empirical power of CIPS when $m<m^0$, which can be reduced if a nonstationary common factor is not accounted for. Although the estimated break date—${\widehat{T}}_B=1984$ and ${\tilde{T}}_B=1985$ for Model B1, and ${\widehat{T}}_B=1991$ and ${\tilde{T}}_B=1986$ for Model C1—is not located around the 2003–2006 housing price bubble period mentioned above, it does reflect important changes that experience the U.S. real state market—that is, the deregulations of the mortgage market, the development of a secondary mortgage market and the increasing role of the government-sponsored enterprises that affected the U.S. economy during the first half of the eighties; see Gerardi, Harvey and Willen (2010) and Ahamada and Diaz-Sanchez (2013). It is interesting to note that all estimated elasticities are below one, with the major change observed for Model C1 with ${\tilde{T}}_B$.

The conclusions that are obtained for the two structural breaks depend on the break estimation procedure. Evidence of panel cointegration is quite weak when using ${\widehat{T}}_B$—that is, the null hypothesis of spurious regression is only rejected for Model C1 with m = 1—and is stronger when using ${\tilde{T}}_B$—panel cointegration is found except for Model C1 with m = 2. It is worth noting that for the latter, the first estimated break date detects the policy changes that experienced the U.S. housing market during the eighties, whereas the second break date is close to the global financial crisis (Model B1) or to the housing prices bubble discussed above (Model C1). As for the parameter estimates, income elasticity shows a mild (large) increase from the first to the second regimes Model B1 (Model C1), to experience a decrease from the second to the third regimes for both models.

Finally, we have also reported the value of the panel BIC statistic that has been proposed in the article as a way to estimate the number of structural breaks. The model specification that minimizes the panel BIC is the one given by Model C1 with J = 1 and ${\tilde{T}}_B$. This result reinforces the analysis the overall discussion about the existence of a long-run relationship between the logarithm of the real housing prices index and the logarithm of the real disposable income per capita for the U.S. States, a conclusion that is robust to the accommodation of one unknown structural break.

7 Conclusions

The article has shown that a consistent estimate of the long-run average coefficient can be obtained when cross-section dependence is present among the panel data units. The type of cross-section dependence that is considered in the article is strong, which is accounted for using an approximate common factor model. The model specification is quite flexible and allows for multiple structural breaks that can affect the deterministic component, the cointegrating vector and/or the loadings of the common factors. The estimation procedure that is applied is based on the CCE approach in Pesaran (2006), which approximates the unobserved common factors using cross-section averages of the observable variables of the model. Our result contributes to the literature of nonstationary panel data analysis, where consistent estimation of the parameters of the model is feasible in a spurious regression framework. The article conducts an extensive simulation exercise to study the finite sample performance of the estimator and test statistic that has been proposed in the article.

The application of the procedures that are designed in the article is illustrated with a model that defines a potential relationship between housing prices and real disposable income per capita. The analysis builds upon the use of U.S. regional data that covers a long time period. The main conclusion that is drawn indicates that robust evidence on the existence of a long-run relationship between housing prices and real disposable income per capita can be obtained once the presence of structural breaks that capture relevant events for the U.S. housing market are allowed for in the model specification.

Supplementary Materials

Appendix A. Mathematical appendix: Collects the proofs of the theorems.

Appendix B. Tables of critical values: Presents the tables of critical values for the individual and panel data cointegration test statistics, considering one structural break.

Appendix C. Monte Carlo experiment and empirical illustration: Pooled CCE estimator, break fraction histograms and empirical size and power of the panel cointegration statistics. Plots of the variables used in the empirical application.

Disclosure Statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Carrion-i-Silvestre acknowledges the financial support from the grant PID2020-114646RB-C41 funded by MCIN/AEI/10.13039/501100011033.

Notes

1 A GAUSS computer code is available upon request to estimate the break dates in panel data with up to two structural breaks, although the Bai and Perron (1998) efficient estimation algorithm which allows global minimizers of SSR that is of order $O(T^2)$ for any $J^0\ge 2$ can also be adapted to our panel data framework. In this regard, Ditzen, Karavias and Westerlund (2023) extend the Bai and Perron (1998) methodology to panel data with the development of a toolbox for Stata software using the Bai and Perron (1998) efficient estimation algorithm adapted to panel data.

2 This situation is equivalent to the one found in time series analysis, see Perron (1990).

3 For Models A and B the factor loadings are assume to be constant across regimes so that the sixth and seventh element of the right hand side of (17) are not affected by the structural breaks—that is, ${\sum }_{j=0}^{\widehat{J}}{\phi }_{i,j}^{\prime }{(D{U}_j\odot \overline{A})}_{t-1}$ has to be replaced by ${\phi }_i^{\prime }{\overline{A}}_{t-1}$ and ${\sum }_{j=0}^{\widehat{J}}{\sum }_{l=0}^p{\kappa }_{i,j,l}^{\prime }{(D{U}_j\odot \Delta \overline{A})}_{t-l}$ by ${\sum }_{l=0}^p{\kappa }_{i,l}^{\prime }\Delta {\overline{A}}_{t-l}$.

4 We abuse notation when using the Hadamard product in the expressions $D{U}_j\odot \overline{A}$ and $D{U}_j\odot \Delta \overline{A}$ that appear in (17) since, in general, the involved matrices might have different column dimensions. In this case, it should be understood that for any given $A_{m\times n}$ and $B_{m\times o}$ matrices $A\odot B={[(A_1\odot B_1)\dots (A_1\odot B_o)\dots (A_n\odot B_1)\dots (A_n\odot B_o)]}_{m\times \mathit{\text{no}}}$. This convention applies throughout the article.

5 An alternative strategy in the case of unknown number of factors, as followed by Pesaran et al. (2013), is to undertake the testing for all permissible values of m⁰ (using all combinations of m⁰ cross-section averages for each choice of m⁰). The size properties of such a procedure are not clear nor are the likely conclusions if one accepts the null hypothesis for some values of m⁰ and rejects for others. This is a topic for future research.

6 For the empirical implementation of this procedure we can follow Kim and Perron (2009) and use a trimming window of 6 observations.

7 A GAUSS program, available upon request, allows the computation of critical values for the multiple structural breaks case for both the individual and panel data cointegration statistics.

8 The specification of the shrinking factor $T^{-1/2-\delta },\delta >0$, is motivated by the fact that $(\widehat{\lambda }-{\lambda }^0)=O_p(T^{-1/2})$ for Models A2, B2, and C2.