Common Correlated Effects Estimation for Dynamic Heterogeneous Panels with Non-Stationary Multi-Factor Error Structures

Abstract

In this paper, we consider the estimation of a dynamic panel data model with non-stationary multi-factor error structures. We adopted the common correlated effect (CCE) estimation and established the asymptotic properties of the CCE and common correlated effects mean group (CCEMG) estimators, as N and T tend to infinity. The results show that both the CCE and CCEMG estimators are consistent and the CCEMG estimator is asymptotically normally distributed. The theoretical findings were supported for small samples by an extensive simulation study, showing that the CCE estimators are robust to a wide variety of data generation processes. Empirical findings suggest that the CCE estimation is widely applicable to models with non-stationary factors. The proposed procedure is also illustrated by an empirical application to analyze the U.S. cigar dataset.

1. Introduction

Recently, there has been increased interest in the analysis of panel data models with cross-sectionally dependent errors (also known as unobserved common factors or multi-factor error structures), which are motivated by empirical applications in economics, such as common shocks and the global financial crisis, see Omay and Kan (2010); Bussiere et al. (2013); Eberhardt et al. (2013) and Chudik et al. (2017), etc. The dependencies across the units violate the traditional assumption of independent and identically distributed errors; conventional panel estimation methods (such as fixed effects estimation) could have serious consequences and lead to inconsistent estimations and misleading inferences. Therefore, in econometrics literature, much effort has been devoted to the estimations for panels with cross-sectional dependence, for example, Pesaran (2006); Bai (2009); Zaffaroni (2009); Greenaway-McGrevy et al. (2012); Kao et al. (2012); Chudik and Pesaran (2013); Moon and Weidner (2015, 2017), among others. See also Chudik and Pesaran (2015b) for a survey of recent developments in large panel models with cross-sectional dependence.

Among these studies, a predominant approach of dealing with cross-sectionally dependent errors in panel models is the so-called common correlated effect (CCE) method proposed by Pesaran (2006)1. The basic idea of CCE estimation is to proxy the unobserved common factors using the cross-sectional averages of the observables in the regression. Comparatively, it has several advantages. For instance, it can be computed by least squares to auxiliary regression, and it does not require the knowledge of the number of unobserved factors. The CCE method has been further developed and applies to different types of panel models. To name a few, Chudik and Pesaran (2015a) suggested the CCE approach to analyze dynamic heterogeneous panels with stationary unobserved common factors. Kapetanios et al. (2011) extended the CCE method to static panel data models with non-stationary multi-factor error structures. Westerlund et al. (2019) considered the CCE for short panels, and Zhou and Zhang (2016) extended the CCE for unbalanced panels.

Among the aforementioned works, there is a gap in the CCE estimation for dynamic panels with non-stationary unobservable factors. To fill this gap, in this paper we consider a linear dynamic heterogeneous panel data model with non-stationary unobserved common factors when both the cross-sectional and time dimensions of the dataset grow to infinity. Under these settings, we find that the CCE estimator of the individual coefficient is consistent, and the CCE mean group (CCEMG) estimator is consistent and has a normal limit distribution. The practical implication of this finding is that for inferential purposes of the CCE estimation, one does not necessarily need to test the stationarity of the unobserved common factors in the model. The finite sample properties are examined through Monte Carlo simulations and the simulation results confirm our theoretical findings in the paper. Moreover, the proposed procedure is illustrated by an empirical application, which analyzes the U.S. cigar dataset.

The rest of the paper is organized as follows. Section 2 sets up the basic model and introduces the CCE estimation of the dynamic heterogeneous panel data model with common factors. The asymptotics of the CCE estimation with non-stationary unobserved common factors is provided in Section 3. Monte Carlo simulation results and an empirical application are reported in Section 4 and Section 5, respectively. The concluding remarks are made in Section 6. Proof of the main results is provided in Appendix A.

Notation: The letter K stands for a finite positive constant. All vectors are column vectors represented by bold lower case letters, and matrices are represented by bold capital letters. Let $∥\mathbf{A}∥=\sqrt{tr\left({\mathbf{AA}}^{\prime }\right)}$ denote the Frobenius norm. ${∥\mathbf{A}∥}_1={max}_{1\le j\le n}{\Sigma }_{i=1}^n\left|a_{ij}\right|$, and ${∥\mathbf{A}∥}_{\infty }={max}_{1\le i\le n}{\Sigma }_{j=1}^n\left|a_{ij}\right|$ denote the maximum absolute column and row sum matrix norms, respectively. ${\mathbf{A}}^{+}$ denotes the Moore–Penrose inverse of $\mathbf{A}$, $rank\left(\mathbf{A}\right)$ and $\varrho \left(\mathbf{A}\right)$ denotes the rank and the spectral radius of $\mathbf{A}$, respectively.

2. Dynamic Panel Data Model with Non-Stationary Unobserved Common Factors

2.1. The Model

We assume the scalar dependent variable $y_{it}$ and regressors ${\mathbf{x}}_{it}$ are generated as follows2

$$y_{it}=c_{yi}+{\varphi }_i{y}_{i,t-1}+{\beta }_i^{\prime }{\mathbf{x}}_{it}+{\gamma }_i^{\prime }{\mathbf{f}}_t+{\epsilon }_{it},$$

(1)

and

$${\mathbf{x}}_{it}={\mathbf{c}}_{xi}+{\alpha }_i{y}_{i,t-1}+{\mathbf{\Gamma }}_i^{\prime }{\mathbf{f}}_t+{\mathbf{u}}_{it},$$

(2)

for $i=1,2,\dots ,N$ and $t=1,2,\dots ,T$, where $c_{yi}$ and ${\mathbf{c}}_{xi}$ are individual fixed effects for unit i, ${\mathbf{x}}_{it}$ is a $k\times 1$ vector of the regressors specific to cross-sectional unit i at time t, ${\epsilon }_{it}$ are the individual-specific (idiosyncratic) errors and ${\mathbf{u}}_{it}$ are the individual-specific components of ${\mathbf{x}}_{it}$, ${\gamma }_i$ and ${\mathbf{\Gamma }}_i$ are $m\times 1$ and $m\times k$ factor loading matrices, and the $m\times 1$ vector ${\mathbf{f}}_t$ represents unobserved common factors. In what follows, we maintain the restriction that model (1) is stationary, such that $0<$$\left|{\varphi }_i\right|<1$ for $i=1,2,\dots ,N.$

Models (1)–(2) have been widely studied in the literature; see, for instance, Pesaran (2006), Chudik and Pesaran (2015a), Westerlund et al. (2019), and the references therein. We follow these studies to consider the CCE estimation for ${\varphi }_i$ and ${\beta }_i,$ and reexamine the validity of the CCE estimation when ${\mathbf{f}}_t$ is non-stationary.

2.2. CCE Estimation

Following Chudik and Pesaran (2015a), let ${\mathbf{z}}_{it}={\left(y_{it},{\mathbf{x}}_{it}^{\prime }\right)}^{\prime }$, then (1) and (2) can be compactly written as

$${\mathbf{z}}_{it}={\mathbf{c}}_{zi}+{\mathbf{A}}_i{\mathbf{z}}_{i,t-1}+{\mathbf{A}}_{0i}^{-1}{\mathbf{C}}_i{\mathbf{f}}_t+{\mathbf{e}}_{zit},$$

(3)

where ${\mathbf{c}}_{zi}={\mathbf{A}}_{0i}^{-1}{\mathbf{c}}_i$, ${\mathbf{A}}_i={\mathbf{A}}_{0i}^{-1}{\mathbf{A}}_{1i}$, ${\mathbf{C}}_i={\left({\gamma }_i,{\mathbf{\Gamma }}_i\right)}^{\prime }$, and ${\mathbf{e}}_{zit}={\mathbf{A}}_{0i}^{-1}{\mathbf{e}}_{it}$, with ${\mathbf{c}}_i={(c_{yi},{\mathbf{c}}_{xi}^{\prime })}^{\prime }$, ${\mathbf{e}}_{it}={\left({\epsilon }_{it},{\mathbf{u}}_{it}^{\prime }\right)}^{\prime }$, and

$${\mathbf{A}}_{0i}=\left(\begin{array}{cc}1& -{\beta }_i^{\prime }\\ \mathbf{0}& {\mathbf{I}}_k\end{array}\right),{\mathbf{A}}_{1i}=\left(\begin{array}{cc}{\varphi }_i& 0\\ {\alpha }_i& \mathbf{0}\end{array}\right).$$

(4)

If the support of $\varrho \left({\mathbf{A}}_i\right)$ lies strictly inside the unit circle, then (3) can be rewritten as the following distributed lag form

$${\mathbf{z}}_{it}=\sum _{l=0}^{\infty }{\mathbf{A}}_i^l\left({\mathbf{c}}_{zi}+{\mathbf{A}}_{0i}^{-1}{\mathbf{C}}_i{\mathbf{f}}_{t-l}+{\mathbf{e}}_{zi,t-l}\right),$$

(5)

for $i=1,2,\dots ,N.$

Taking the cross-sectional average of (5) yields

$${\overline{\mathbf{z}}}_t={\overline{\mathbf{c}}}_z+\mathbf{\Lambda }\left(L\right){\mathbf{Cf}}_t+O_p\left(N^{-1/2}\right),$$

where ${\overline{\mathbf{z}}}_t=\frac{1}{N}{\sum }_{i=1}^N{\mathbf{z}}_{it}$ is a $k+1$ dimensional vector of the cross-section average, ${\overline{\mathbf{c}}}_z=\frac{1}{N}{\sum }_{i=1}^N{({\mathbf{I}}_{k+1}-{\mathbf{A}}_i)}^{-1}{\mathbf{c}}_{zi}$, and $\mathbf{C}=E\left({\mathbf{C}}_i\right)={\left(\gamma ,\mathbf{\Gamma }\right)}^{\prime }$, $\mathbf{\Lambda }\left(L\right)={\sum }_{l=0}^{\infty }{\mathbf{\Lambda }}_l{L}^l$ with ${\mathbf{\Lambda }}_l=E\left({\mathbf{A}}_i^l{\mathbf{A}}_{0i}^{-1}\right)$, and L being the lag operator. Furthermore, if $\mathbf{\Lambda }\left(L\right)$ is invertible (see Assumption 4 below), then we have

$${\mathbf{Cf}}_t={\mathbf{\Lambda }}^{-1}\left(L\right)({\overline{\mathbf{z}}}_t-{\overline{\mathbf{c}}}_z)+O_p\left(N^{-1/2}\right).$$

When the $(k+1)\times m$ matrix $\mathbf{C}$ has the full column rank, i.e., the rank condition

$$rank\left(\mathbf{C}\right)=m\le k+1,$$

(6)

holds, we have

$${\mathbf{f}}_t=\mathbf{B}\left(L\right)({\overline{\mathbf{z}}}_t-{\overline{\mathbf{c}}}_z)+O_p\left(N^{-1/2}\right),$$

(7)

where $\mathbf{B}\left(L\right)={\left({\mathbf{C}}^{\prime }\mathbf{C}\right)}^{-1}{\mathbf{C}}^{\prime }{\mathbf{\Lambda }}^{-1}\left(L\right)$. This suggests that the contemporary and lagged value of ${\overline{\mathbf{z}}}_t={\left({\overline{y}}_t,{\overline{\mathbf{x}}}_t^{\prime }\right)}^{\prime }$ can be used as observable proxies for the unobserved common factors ${\mathbf{f}}_t$.

Substituting the observed proxies of the unobserved common factors (7) into (1) yields the following augmented regression

$$y_{it}=c_{yi}^{*}+{\varphi }_i{y}_{i,t-1}+{\beta }_i^{\prime }{\mathbf{x}}_{it}+{\delta }_i^{\prime }\left(L\right){\overline{\mathbf{z}}}_t+{\epsilon }_{it}+O_p\left(N^{-1/2}\right),$$

or

$$y_{it}=c_{yi}^{*}+{\varphi }_i{y}_{i,t-1}+{\beta }_i^{\prime }{\mathbf{x}}_{it}+\sum _{l=0}^{p_T}{\delta }_{il}^{\prime }{\overline{\mathbf{z}}}_{t-l}+w_{it},$$

(8)

for $t=p_T+1,p_T+2,\dots ,T$, where $c_{yi}^{*}=c_{yi}-{\delta }_i^{\prime }\left(1\right){\overline{\mathbf{c}}}_z$, ${\delta }_i\left(L\right)={\mathbf{B}}^{\prime }\left(L\right){\gamma }_i={\sum }_{l=0}^{\infty }{\delta }_{il}{L}^l$, $p_T$ is the number of lags used to truncate the infinite polynomial distributed lag function ${\delta }_i\left(L\right)$,3 and the composite error $w_{it}$ has the form of

$$w_{it}={\epsilon }_{it}+\sum _{l=p_T+1}^{\infty }{\delta }_{il}^{\prime }{\overline{\mathbf{z}}}_{t-l}+O_p\left(N^{-1/2}\right).$$

For notational simplicity, let ${\mathbf{y}}_i={\left(y_{i,p_T+1},y_{i,p_T+2},\dots ,y_{iT}\right)}^{\prime }$, ${\mathbf{\Xi }}_i=({\mathbf{y}}_{i,-1},{\mathbf{X}}_i)$, with ${\mathbf{y}}_{i,-1}={\left(y_{i,p_T},y_{i,p_T+1},\dots ,y_{i,T-1}\right)}^{\prime }$ and ${\mathbf{X}}_i={({\mathbf{x}}_{i,p_T+1},{\mathbf{x}}_{i,p_T+2},\dots ,{\mathbf{x}}_{iT})}^{\prime }$, and ${\mathbf{w}}_i={\left(w_{i,p_T+1},w_{i,p_T+2},\dots ,w_{iT}\right)}^{\prime }$, then the augmented regression (8) can be expressed in vector form as

$${\mathbf{y}}_i={\mathbf{\Xi }}_i{\pi }_i+\overline{\mathbf{Q}}{\mathbf{d}}_i+{\mathbf{w}}_i,$$

(9)

where ${\pi }_i={\left({\varphi }_i,{\beta }_i^{\prime }\right)}^{\prime }$ are the parameters of interest, ${\mathbf{d}}_i={(c_{yi}^{*},{\delta }_{i0}^{\prime },{\delta }_{i1}^{\prime },\dots ,{\delta }_{ip}^{\prime })}^{\prime }$ are nuisance parameters, and4

$$\overline{\mathbf{Q}}=\left(\begin{array}{ccccc}1& {\overline{\mathbf{z}}}_{p_T+1}^{\prime }& {\overline{\mathbf{z}}}_{p_T}^{\prime }& \cdots & {\overline{\mathbf{z}}}_1^{\prime }\\ 1& {\overline{\mathbf{z}}}_{p_T+2}^{\prime }& {\overline{\mathbf{z}}}_{p_T+1}^{\prime }& \cdots & {\overline{\mathbf{z}}}_2^{\prime }\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\overline{\mathbf{z}}}_T^{\prime }& {\overline{\mathbf{z}}}_{T-1}^{\prime }& \cdots & {\overline{\mathbf{z}}}_{T-p_T}^{\prime }\end{array}\right).$$

(10)

Based on the cross-sectionally augmented regression model (9) and by the formula for partitioned regression, the CCE estimator of the individual coefficients ${\pi }_i$ is given by

$${\widehat{\pi }}_i={\left({\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\mathbf{\Xi }}_i\right)}^{-1}{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\mathbf{y}}_i,$$

(11)

which is an ordinary least squares estimate, where ${\mathbf{M}}_q={\mathbf{I}}_{T-p_T}-\overline{\mathbf{Q}}{\left({\overline{\mathbf{Q}}}^{\prime }\overline{\mathbf{Q}}\right)}^{+}{\overline{\mathbf{Q}}}^{\prime }$ is an orthogonal projection matrix, with ${\mathbf{I}}_{T-p_T}$ a $(T-p_T)$-dimensional identity matrix. In panel models with N large, the primary parameters of interest are the means of the individual–specific coefficients, $E\left({\pi }_i\right)=\pi$, which can be estimated by the common correlated effects mean group (CCEMG) estimator

$${\widehat{\pi }}_{MG}=\frac{1}{N}\sum _{i=1}^N{\widehat{\pi }}_i.$$

(12)

3. Asymptotics of CCE Estimators with Non-Stationary Factors

3.1. Assumptions

When the unobserved common factors ${\mathbf{f}}_t$ are stationary processes, Chudik and Pesaran (2015a) showed that the CCE estimator (11) of the individual coefficient is consistent, and the CCEMG estimator (12) is consistent and asymptotically normal. However, in practice, the common factors ${\mathbf{f}}_t$ may follow a non-stationary process (see Bai and Ng 2004, 2010; Pesaran 2007; Pesaran et al. 2013, among others). In this scenario, the validity of CCE estimators and their asymptotic properties need to be re-examined.

Following Kapetanios et al. (2011), we assume the unobserved common factors follow the multivariate unit root process

$${\mathbf{f}}_t={\mathbf{f}}_{t-1}+{\varsigma }_t.$$

(13)

To derive the asymptotic properties of the CCE type estimators (11) and (12) when ${\mathbf{f}}_t$ follows (13), we make the following assumptions.

Assumption 1.

(Individual-specific errors). (i) The individual–specific errors${\epsilon }_{it}$follow a linear stationary process with uniformly-bounded positive variance,${sup}_i{\sigma }_i^2<K$, for some constant K, and uniformly-bounded fourth-order cumulants.${\mathbf{u}}_{it}$follows a linear stationary process with absolute summable auto-covariances (uniformly in i), with covariance matrices,${\Sigma }_{{\mathbf{u}}_i}$, which are non-singular and satisfy${sup}_i∥{\Sigma }_{{\mathbf{u}}_i}∥<K$, and have uniformly-bounded fourth-order cumulants. (ii)${\epsilon }_{it}$are independently distributed of${\mathbf{u}}_{j{t}^{\prime }}$for all$i,j,t$, and$t^{\prime }$. For each i,$e_{it}={({\epsilon }_{it},{\mathbf{u}}_{it}^{\prime })}^{\prime }$is an$(k+1)\times 1$vector of$L_{2+\delta },\delta >0$, stationary near epoch dependent processes of size$2\delta /(2\delta +4)$on the$\alpha$-mixing process of size$-(2+\delta )/\delta$, and for$i=1,2,\dots ,N$,$Var\left(e_{it}\right)={\Sigma }_{{\mathbf{e}}_i}$, which is a non-singular matrix and satisfies${sup}_i∥{\Sigma }_{{\mathbf{e}}_i}∥<K$.

Assumption 2.

(Factor loadings). The factor loadings ${\gamma }_i$ and ${\mathbf{\Gamma }}_i$ are independently and identically distributed ($IID$) across i, and of the common factors ${\mathbf{f}}_t$, for all i and t, with means $\gamma$ and $\mathbf{\Gamma }$, respectively, and the bounded second moments. In particular,

$${\gamma }_i=\gamma +{\eta }_{\gamma i},\mathrm{with}{\eta }_{\gamma i}\sim IID(\mathbf{0},{\Sigma }_{\gamma }),\mathrm{for}i=1,2,\dots ,N,$$

and

$$vec\left({\mathbf{\Gamma }}_i\right)=vec\left(\mathbf{\Gamma }\right)+{\eta }_{\mathbf{\Gamma }i},\mathrm{with}{\eta }_{\mathbf{\Gamma }i}\sim IID(\mathbf{0},{\Sigma }_{\mathbf{\Gamma }}),\mathrm{for}i=1,2,\dots ,N,$$

where ${\Sigma }_{\gamma }$ and ${\Sigma }_{\mathbf{\Gamma }}$ are $m\times m$ and $mk\times mk$ symmetric nonnegative definite matrices, $∥\gamma ∥<K,∥{\Sigma }_{\gamma }∥<K,∥\mathbf{\Gamma }∥<K$ and $∥{\Sigma }_{\mathbf{\Gamma }}∥<K$, for some constant $K.$

Assumption 3.

(Heterogeneous coefficients). The slope coefficients ${\pi }_i={\left({\varphi }_i,{\beta }_i^{\prime }\right)}^{\prime }$ follow the random coefficient model

$${\pi }_i=\pi +{\upsilon }_{\pi i},{\upsilon }_{\pi i}\sim IID(\mathbf{0},{\Sigma }_{\pi }),\mathrm{for}i=1,2,\dots ,N,$$

where $\pi =E\left({\pi }_i\right)={\left(\varphi ,{\beta }^{\prime }\right)}^{\prime }$, $∥\pi ∥<K,∥{\Sigma }_{\pi }∥<K$, ${\Sigma }_{\pi }$ is the $(k+1)\times (k+1)$ symmetric nonnegative definite matrix and the random deviations ${\upsilon }_{\pi i}$ are distributed independently of ${\gamma }_j,{\mathbf{\Gamma }}_j,{\epsilon }_{jt},{\mathbf{u}}_{jt}$ and ${\varsigma }_t$ for $i,j$, and t. Furthermore, the support of ${\varphi }_i$ lies strictly inside the unit circle, and $E∥{\mathbf{c}}_i∥<K$, $E∥{\alpha }_i∥<K$ for all i, where ${\mathbf{c}}_i={(c_{yi},{\mathbf{c}}_{xi}^{\prime })}^{\prime }$.

Assumption 4.

(Exogenous regressors). Regressors ${\mathbf{x}}_{it}$ are either strictly exogenous and generated according to the canonical factor model (2) with ${\alpha }_i=\mathbf{0}$, or weakly exogenous and generated according to (2) with ${\alpha }_i$, for $i=1,2,\dots ,N$, $IID$ across i, and independently distributed of ${\upsilon }_{\pi j}$, ${\gamma }_j,{\mathbf{\Gamma }}_j,{\epsilon }_{jt},{\mathbf{u}}_{jt}$ and ${\mathbf{f}}_t$ for all $i,j$, and t. In the case where the regressors are weakly exogenous, we also assume:

(i) The support of $\varrho \left({\mathbf{A}}_i\right)$ lies strictly inside the unit circle, for $i=1,2,\dots ,N$, where ${\mathbf{A}}_i={\mathbf{A}}_{0i}^{-1}{\mathbf{A}}_{1i}$ with ${\mathbf{A}}_{0i}$ and ${\mathbf{A}}_{1i}$ are defined in (4).

(ii) The inverse of polynomial $\mathbf{\Lambda }\left(L\right)={\sum }_{l=0}^{\infty }{\mathbf{\Lambda }}_l{L}_l$ exists and has exponentially decaying coefficients, where ${\mathbf{\Lambda }}_l=E\left({\mathbf{A}}_i^l{\mathbf{A}}_{0i}^{-1}\right)$.

Assumption 5.

(Rank condition). The $(k+1)\times m$ matrix $\mathbf{C}$ has a full column rank, such that

$$rank\left(\mathbf{C}\right)=m\le k+1,$$

where $\mathbf{C}=E\left({\mathbf{C}}_i\right)=E{\left({\gamma }_i,{\mathbf{\Gamma }}_i\right)}^{\prime }.$

Assumption 6.

(i) As $\left(N,T\right)\to \infty ,$ the $(k+1)\times (k+1)$ matrices ${\mathbf{\Psi }}_{iT}^{-1}={({\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\mathbf{\Xi }}_i/T)}^{-1},$ ${\mathbf{\Psi }}_{ih}^{-1}={({\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_h{\mathbf{\Xi }}_i/T)}^{-1}$ and ${\mathbf{\Psi }}_{ig}^{-1}={({\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_g{\mathbf{\Xi }}_i/T)}^{-1}$ exist for all i, and ${\mathbf{\Psi }}_{iT}^{-1},$ ${\mathbf{\Psi }}_{ih}^{-1}$ and ${\mathbf{\Psi }}_{ig}^{-1}$ have finite second-order moments for all $i,$ where ${\mathbf{M}}_h={\mathbf{I}}_{T-p_T}-\mathbf{H}{\left({\mathbf{H}}^{\prime }\mathbf{H}\right)}^{+}{\mathbf{H}}^{\prime }$ and ${\mathbf{M}}_g={\mathbf{I}}_{T-p_T}-\mathbf{G}{\left({\mathbf{G}}^{\prime }\mathbf{G}\right)}^{+}{\mathbf{G}}^{\prime }$ are projection matrices, where ${\mathbf{A}}^{+}$ denotes the Moore–Penrose generalized inverse of $\mathbf{A},$ $\mathbf{H}$, and $\mathbf{G}$, defined as

$$\mathbf{H}=\left(\begin{array}{cccc}1& {\mathbf{h}}_{p_T+1}^{\prime }& \cdots & {\mathbf{h}}_1^{\prime }\\ 1& {\mathbf{h}}_{p_T+2}^{\prime }& \cdots & {\mathbf{h}}_2^{\prime }\\ ⋮& ⋮& \ddots & ⋮\\ 1& {\mathbf{h}}_T^{\prime }& \cdots & {\mathbf{h}}_{T-p_T}^{\prime }\end{array}\right),\mathbf{G}=(\tau ,\tilde{\mathbf{F}})=\left(\begin{array}{cccc}1& {\mathbf{f}}_{p_{{}_T}+1}^{\prime }& \cdots & {\mathbf{f}}_1^{\prime }\\ 1& {\mathbf{f}}_{p_{{}_T}+2}^{\prime }& \cdots & {\mathbf{f}}_2^{\prime }\\ ⋮& ⋮& \ddots & ⋮\\ 1& {\mathbf{f}}_T^{\prime }& \cdots & {\mathbf{f}}_{T-P_T}^{\prime }\end{array}\right),$$

where ${\mathbf{h}}_t=\mathbf{\Psi }\left(L\right){\mathbf{f}}_t+{\overline{\mathbf{c}}}_z$ with $\mathbf{\Psi }\left(L\right)=\frac{1}{N}{\sum }_{i=1}^N{({\mathbf{I}}_{k+1}-{\mathbf{A}}_{i}L)}^{-1}{\mathbf{A}}_{0i}^{-1}{\mathbf{C}}_i$ and ${\overline{\mathbf{c}}}_z=\frac{1}{N}{\sum }_{i=1}^N{({\mathbf{I}}_{k+1}-{\mathbf{A}}_i)}^{-1}{\mathbf{c}}_{zi}$, $\tau ={(1,1,\dots 1)}^{\prime }$ is a $(T-p_T)\times 1$ vector of ones.

(ii) The matrix ${\mathbf{\Psi }}^{*}={lim}_{N\to \infty }\frac{1}{N}{\sum }_{i=1}^N{\Sigma }_{{\Omega }_i}$ is non-singular, where ${\Sigma }_{\Omega i}={\mathbf{\Psi }}_{\xi i}\left(\begin{array}{cc}{\sigma }_i^2& \mathbf{0}\\ \mathbf{0}& {\Sigma }_{{\mathbf{u}}_i}\end{array}\right){\mathbf{\Psi }}_{\xi i}^{\prime }$ with ${\mathbf{\Psi }}_{\xi i}=\left(\begin{array}{c}{\mathbf{S}}_y^{\prime }{({\mathbf{I}}_{k+1}-{\mathbf{A}}_{i}L)}^{-1}L\\ {\mathbf{S}}_x^{\prime }{({\mathbf{I}}_{k+1}-{\mathbf{A}}_{i}L)}^{-1}\end{array}\right){\mathbf{A}}_{0i}^{-1}$ a $(k+1)\times (k+1)$ matrix, ${\mathbf{S}}_y^{\prime }={\left(\begin{array}{cc}1& \mathbf{0}\end{array}\right)}_{1\times \left(k+1\right)}$, and ${\mathbf{S}}_x^{\prime }={\left(\begin{array}{cc}\mathbf{0}& {\mathbf{I}}_k\end{array}\right)}_{k\times \left(k+1\right)}$.

Assumption 7.

${\varsigma }_t$ in (13) is an $m\times 1$ vector of $L_{2+\delta },\delta >0,$ stationary near epoch-dependent processes of size $1/2$, on an $\alpha$-mixing process of size $-(2+\delta )/\delta$, and is distributed independently of the idiosyncratic errors ${\epsilon }_{it}$ and ${\mathbf{u}}_{it}$ for all i and t.

Several remarks can be made for these assumptions. Assumptions 1–3 are quite standard in the literature for (dynamic) panel models with cross-sectional dependence, for example, see Pesaran (2006) and Kapetanios et al. (2011) and the references therein. Assumption 4 is also made on Chudik and Pesaran (2015a) for exogenous regressors and stationarity conditions for dynamic panels. Assumption 5 is a common condition for the implementation of the CCE estimation (e.g., Pesaran (2006) and Chudik and Pesaran (2015a), etc.), which implies that there are more included regressors than the unobserved factors in the model. See Juodis et al. (2021) for a detailed discussion of the validity of the rank condition and the resulting asymptotics for the CCE estimation. Assumption 6 is a common assumption for the CCE estimation and it is imposed for the partition regression in augmented regression for the dynamic panels (e.g., Chudik and Pesaran 2015a). Assumption 7 requires that the error structures in the unit root process ${\mathbf{f}}_t$ are stationary.

3.2. Asymptotics

Under these assumptions, we can establish the asymptotic properties of CCE estimators (11) and (12) when ${\mathbf{f}}_t$ is non-stationary. To begin with, we note that for the original model (1), it can be rewritten as in the vector form

$${\mathbf{y}}_i=c_{yi}\tau +{\varphi }_i{\mathbf{y}}_{i,-1}+{\mathbf{X}}_i{\beta }_i+\mathbf{F}{\gamma }_i+{\epsilon }_i,$$

(14)

or more compactly as

$${\mathbf{y}}_i=c_{yi}\tau +{\mathbf{\Xi }}_i{\pi }_i+\mathbf{F}{\gamma }_i+{\epsilon }_i,$$

(15)

for $i=1,2,\dots ,N$, where ${\mathbf{y}}_i={\left(y_{i,p_T+1},y_{i,p_T+2},\dots ,y_{iT}\right)}^{\prime }$, ${\mathbf{X}}_i={({\mathbf{x}}_{i,p_T+1},\dots ,{\mathbf{x}}_{iT})}^{\prime }$, ${\mathbf{y}}_{i,-1}={\left(y_{i,p_T},y_{i,p_T+1},\dots ,y_{i,T-1}\right)}^{\prime }$, ${\mathbf{\Xi }}_i=({\mathbf{y}}_{i,-1},{\mathbf{X}}_i)$, $\mathbf{F}={({\mathbf{f}}_{p_T+1},\dots ,{\mathbf{f}}_T)}^{\prime }$, and ${\epsilon }_i={({\epsilon }_{i,p_T+1},\dots ,{\epsilon }_{iT})}^{\prime }$.

Using the CCE estimator (11) into (15), we have

$${\widehat{\pi }}_i-{\pi }_i={\left(\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\mathbf{\Xi }}_i}{T}\right)}^{-1}\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q\mathbf{F}{\gamma }_i}{T}+{\left(\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\mathbf{\Xi }}_i}{T}\right)}^{-1}\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\epsilon }_i}{T},$$

(16)

which shows that the asymptotics of ${\widehat{\pi }}_i$ depends on the unobserved factors through $T^{-1}{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q\mathbf{F}$.

Using the results in Lemma A2, A5, and A6 in Appendix A, we obtain

$$\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q\mathbf{F}}{T}=O_p\left(\frac{1}{N}\right)+O_p\left(\frac{1}{\sqrt{NT}}\right),\mathrm{uniformly}\mathrm{over}i,$$

(17)

and, thus,

$${\widehat{\pi }}_i-{\pi }_i={\left(\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\mathbf{\Xi }}_i}{T}\right)}^{-1}\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\epsilon }_i}{T}+O_p\left(\frac{1}{N}\right)+O_p\left(\frac{1}{\sqrt{NT}}\right),$$

(18)

when the rank condition (6) is satisfied. The above results are summarized in the following theorem, establishing the consistency of the CCE estimator of individual coefficients of interest.

Theorem 1.

Consider the panel models (1) and (2), suppose Assumption 1–7 hold, then, as $(N,T,p_T)\stackrel{j}{\to }\infty$, such that $p_T^3/T\to \lambda ,$$0<\lambda <\infty$, we have

$${\widehat{\pi }}_i-{\pi }_i\stackrel{p}{\to }0.$$

See the Appendix A for the proof.

Remark 1.

The above theorem suggests that the CCE estimator of the individual slope coefficient is consistent even if the common factors are non-stationary. When the rank condition (6) is not satisfied, the CCE estimator of the individual slope coefficients would be inconsistent due to the correlation of ${\mathbf{x}}_{it}$ and ${\mathbf{f}}_t$. See Juodis et al. (2021) for more discussions on the validity of the CCE estimator when the rank condition does not hold.

Next, we establish the asymptotic properties of the CCEMG estimator of the mean group coefficients, $\pi =E\left({\pi }_i\right)$. We have

$$\begin{array}{ccc}\hfill \sqrt{N}({\widehat{\pi }}_{MG}-\pi )& =& \frac{1}{\sqrt{N}}\sum _{i=1}^N{\upsilon }_{\pi i}+\frac{1}{\sqrt{N}}\sum _{i=1}^N{\left(\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\mathbf{\Xi }}_i}{T}\right)}^{-1}\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q\mathbf{F}}{T}{\gamma }_i\hfill \\ & & +\frac{1}{\sqrt{N}}\sum _{i=1}^N{\left(\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\mathbf{\Xi }}_i}{T}\right)}^{-1}\frac{{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q{\epsilon }_i}{T}.\hfill \end{array}$$

When the rank condition is satisfied, by (17), we have

$$\frac{\sqrt{N}{\mathbf{\Xi }}_i^{\prime }{\mathbf{M}}_q\mathbf{F}}{T}=O_p\left(\frac{1}{\sqrt{N}}\right)+O_p\left(\frac{1}{\sqrt{T}}\right),$$

hence, we can obtain

$$\sqrt{N}({\widehat{\pi }}_{MG}-\pi )=\frac{1}{\sqrt{N}}\sum _{i=1}^N{\upsilon }_{\pi i}+O_p\left(\frac{1}{\sqrt{N}}\right)+O_p\left(\frac{1}{\sqrt{T}}\right)\stackrel{d}{\sim }\frac{1}{\sqrt{N}}\sum _{i=1}^N{\upsilon }_{\pi i},$$

and, thus,

Theorem 2.

Consider the panel models (1) and (2), suppose Assumptions 1–7 hold, as $(N,T,p_T)\stackrel{j}{\to }\infty$, such that $p_T^3/T\to \lambda ,$$0<\lambda <\infty$, then we have

$${\widehat{\pi }}_{MG}-\pi \stackrel{p}{\to }0.$$

If it is further assumed that $N/T\to \phi ,0<\phi <\infty$, then

$$\sqrt{N}({\widehat{\pi }}_{MG}-\pi )\stackrel{d}{\to }N(0,{\Sigma }_{MG}).$$

The asymptotic variance of ${\widehat{\pi }}_{MG}$ can be consistently estimated nonparametrically by

$${\widehat{\Sigma }}_{MG}=\frac{1}{N-1}\sum _{i=1}^N({\widehat{\pi }}_i-{\widehat{\pi }}_{MG}){({\widehat{\pi }}_i-{\widehat{\pi }}_{MG})}^{\prime }.$$

For the results in both Theorems 1 and 2, we find that, for models with non-stationary common factors, although the intermediate results needed for deriving the asymptotic properties of the common correlated effects estimators significantly differ from the stationary case, as in Chudik and Pesaran (2015a), the final results are surprisingly similar. This is in direct contrast to the usual phenomenon where distributional results of $I\left(1\right)$ processes are radically different from those of $I\left(0\right)$ processes.

Remark 2.

For the consistency of ${\widehat{\pi }}_i$ and ${\widehat{\pi }}_{MG}$, no restrictions on the relative expansion rates of N and T to infinity are required. However, they require $N/T\to \phi ,0<\phi <\infty$ for the derivation of the asymptotic distribution of ${\widehat{\pi }}_{MG}$ due to the time series bias, which arises from the presence of lagged values of the dependent variable; therefore, it is unsuitable for panels with T being small relative to N.

Including a lagged dependent variable as the regressor in the model could induce the estimators with time series bias of order $O\left(T^{-1}\right)$. When T is not large, the bias is non-negligible; hence, a certain bias correction approach should be considered. In the simulations below, we consider the Jackknife bias-corrected method for bias reduction (e.g., see (21) below), which is used extensively in the relevant literature (e.g., Hahn and Newey 2004).

4. Monte Carlo Simulation

In this section, we investigate the finite sample properties of the CCEMG estimation for dynamic heterogeneous panels with non-stationary common factors. We consider the following data-generating processes5

$$y_{it}=c_{yi}+{\varphi }_i{y}_{i,t-1}+{\beta }_{0i}{x}_{it}+{\beta }_{1i}{x}_{i,t-1}+{\gamma }_i^{\prime }{\mathbf{f}}_t+{\epsilon }_{it},$$

(19)

and

$$x_{it}=c_{xi}+{\alpha }_{xi}{y}_{i,t-1}+{\gamma }_{xi}^{\prime }{\mathbf{f}}_t+u_{it}.$$

(20)

for $i=1,2,\dots ,N,$ and $t=-99,\dots ,0,\dots ,T.$ Let ${\varphi }_i\sim IIDU\left(0,0.8\right)$, ${\beta }_{0i}\sim IIDU\left(0.5,1\right)$, ${\beta }_{1i}=-0.5$ and ${\alpha }_{xi}\sim IIDU\left(0,0.35\right)$, $c_{yi}\sim IIDN\left(1,1\right)$, $c_{xi}=c_{yi}+{ϵ}_{c_{xi}}$ with ${ϵ}_{c_{xi}}\sim IIDN\left(0,1\right)$. The main purpose of this paper is to illustrate the validity of the CCEMG estimator in the case of non-stationary unobserved common factors; hence, for the unobserved common factors ${\mathbf{f}}_t$, we consider the following three different non-stationary DGPs:

DGP 1. Two non-stationary unobserved common factors $\left(m=2\right)$, $f_{lt}=f_{l,t-1}+{\varsigma }_{flt},$${\varsigma }_{flt}\sim IIDN\left(0,{\sigma }_{fl}^2\right),$ where ${\sigma }_{fl}=0.2$, for $l=1,2$, and $t=-99,\dots ,0,\dots ,T$.

DGP 2. One non-stationary unobserved common factor and a stationary common factor $\left(m=2\right)$, $f_{1t}=f_{1,t-1}+{\varsigma }_{flt},$$f_{2t}=0.6f_{2,t-1}+{\varsigma }_{flt},$${\varsigma }_{flt}\sim IIDN\left(0,{\sigma }_{fl}^2\right),$ where ${\sigma }_{fl}=0.5$, for $l=1,2$, and $t=-99,\dots ,0,\dots ,T$.

DGP 3. Cointegrated unobserved common factors $\left(m=2\right)$, $f_{1t}=f_{1,t-1}+{\varsigma }_{f1t},$$f_{2t}=0.5f_{1t}+{\varsigma }_{f2t},$${\varsigma }_{flt}\sim IIDN\left(0,{\sigma }_{fl}^2\right),$ where ${\sigma }_{fl}=1$, for $l=1,2$, and $t=-99,\dots ,0,\dots ,T$.

For the above DGPs, the starting values are $f_{l,-100}=0$, for $l=1,2$; the first 100 observations are discarded.

Correspondingly, the factor loadings are generated independently across replications as

$${\gamma }_{il}={\gamma }_l+{\eta }_{i,\gamma l},{\eta }_{i,\gamma l}\sim IIDN\left(0,{\sigma }_{\gamma l}^2\right),$$

and

$${\gamma }_{xil}={\gamma }_{xl}+{\eta }_{i,\gamma xl},{\eta }_{i,\gamma xl}\sim IIDN\left(0,{\sigma }_{\gamma xl}^2\right),$$

for $i=1,2,\dots ,N$ and $l=1,2$, where ${\sigma }_{\gamma l}^2={0.2}^2$, ${\sigma }_{\gamma xl}^2={0.3}^2$, and ${\gamma }_l=\sqrt{b_{\gamma l}}$, ${\gamma }_{xl}=\sqrt{l{b}_{xl}}$ for $l=1,2$, where $b_{\gamma l}=1/2-{\sigma }_{\gamma l}^2$ and $b_{xl}=1/2-{\sigma }_{\gamma xl}^2$ for $l=1,2$.

For the idiosyncratic errors, ${\epsilon }_{it}\sim IIDN\left(0,1\right)$ for all i and t, and the unit-specific components $u_{it}$ are generated as independent stationary AR(1) processes:

$$u_{it}={\rho }_{xi}{u}_{i,t-1}+{ϵ}_{xit},{\rho }_{xi}\sim IIDU\left(0,0.95\right),{ϵ}_{xit}\sim IIDN\left(0,1\right),$$

for $i=1,2,\dots ,N$ and $t=-99,\dots ,0,\dots ,T$ with the starting values $u_{i,-100}=0$. The first 100 observations are discarded.

We consider the combination of $N=50,100,200$, and $T=50,100,150,200$. The number of replications is set at 2000 times. In what follows, we focus on the lagged coefficient $\varphi$ (the cross-section mean of ${\varphi }_i$), as well as ${\beta }_0$ (the cross-section mean of ${\beta }_{0i}$). To save space, we only report the results of ${\beta }_0$ since the results for ${\beta }_1$ are very similar to that of ${\beta }_0$ and they are available upon request.

Two estimators are considered in the simulation. The first is the main result of the CCEMG estimator ${\widehat{\pi }}_{MG}$ given in (12), in which, the lag order $p_T$ is selected to satisfy $p_T^3/T\to \lambda ,$ as $T\to \infty$, for some $0<\lambda <\infty$; that is, $p_T=\left[T^{1/3}\right]$, which works well in our Monte Carlo design6. The second is the Jackknife bias-corrected CCEMG estimator, which is constructed as

$${\widehat{\pi }}_{Jack-MG}=2{\widehat{\pi }}_{MG}-\frac{1}{2}\left({\widehat{\pi }}_{MG}^a+{\widehat{\pi }}_{MG}^b\right),$$

(21)

where ${\widehat{\pi }}_{MG}^a$ is the CCEMG estimator calculated using the first two-thirds of the available time period, namely over the period $t=1,2,\dots ,[2T/3]$, and ${\widehat{\pi }}_{MG}^b$ denotes the CCEMG estimator computed using the observations over the period $t=[T/3],[T/3]+1,\dots ,T$, where $[T/3]$ denotes the integer part of $T/3$. Note that a new strategy is applied to improve the performance of the Jackknife estimator, i.e., the whole time period is divided into three parts, the first two-thirds of the available period is applied to calculate the first estimator and another one is computed from the last two-thirds of the period. We find that, in our settings, this division strategy performs better than the half-panel Jackknife method discussed in Chudik and Pesaran (2015a).

We used the statistical software MATLAB to conduct the Monte Carlo experiments; the simulation results are summarized in

Table 1. Estimation results for DGP 1.

	Bias					RMSE
Parameter	$(\mathit{N},\mathit{T})$	50	100	150	200	50	100	150	200
$\varphi$	CCEMG estimation
	50	−0.1065	−0.0393	−0.0163	−0.0004	0.1131	0.0530	0.0392	0.0371
	100	−0.1085	−0.0392	−0.0156	−0.0004	0.1120	0.0476	0.0311	0.0284
	200	−0.1105	−0.0402	−0.0163	−0.0012	0.1116	0.0450	0.0265	0.0220
	Jackknife bias-corrected CCEMG estimation
	50	−0.0508	−0.0126	−0.0071	−0.0043	0.0747	0.0417	0.0415	0.0440
	100	−0.0411	−0.0124	−0.0036	0.0034	0.0689	0.0324	0.0309	0.0365
	200	−0.0417	−0.0127	−0.0042	−0.0039	0.0664	0.0273	0.0253	0.0309
${\beta }_0$	CCEMG estimation
	50	0.0136	0.0071	0.0032	0.0012	0.0461	0.0332	0.0282	0.0275
	100	0.0129	0.0058	0.0029	0.0008	0.0341	0.0232	0.0200	0.0192
	200	0.0119	0.0049	0.0024	0.0003	0.0252	0.0169	0.0150	0.0139
	Jackknife bias-corrected CCEMG estimation
	50	0.0112	0.0043	0.0011	−0.0007	0.0550	0.0361	0.0307	0.0289
	100	0.0098	0.0030	0.0003	−0.0015	0.0397	0.0251	0.0215	0.0206
	200	0.0091	0.0020	0.0000	0.0017	0.0281	0.0180	0.0160	0.0150

,

Table 2. Estimation results for DGP 2.

	Bias					RMSE
Parameter	$(\mathit{N},\mathit{T})$	50	100	150	200	50	100	150	200
$\varphi$	CCEMG estimation
	50	−0.0983	−0.0384	−0.0188	−0.0093	0.1053	0.0513	0.0386	0.0348
	100	−0.1004	−0.0389	−0.0200	−0.0104	0.1040	0.0461	0.0312	0.0259
	200	−0.1015	−0.0395	−0.0193	−0.0102	0.1036	0.0434	0.0271	0.0203
	Jackknife bias-corrected CCEMG estimation
	50	−0.0506	−0.0179	−0.0069	−0.0014	0.0840	0.0410	0.0360	0.0346
	100	−0.0491	−0.0172	−0.0072	0.0015	0.0768	0.0316	0.0262	0.0245
	200	−0.0457	−0.0167	−0.0071	−0.0013	0.0704	0.0257	0.0194	0.0177
${\beta }_0$	CCEMG estimation
	50	0.0124	0.0077	0.0048	0.0036	0.0451	0.0334	0.0285	0.0275
	100	0.0122	0.0063	0.0042	0.0033	0.0335	0.0235	0.0202	0.0193
	200	0.0112	0.0056	0.0039	0.0028	0.0248	0.0171	0.0151	0.0138
	Jackknife bias-corrected CCEMG estimation
	50	0.0109	0.0061	0.0037	0.0027	0.0535	0.0362	0.0306	0.0284
	100	0.0104	0.0045	0.0026	0.0020	0.0396	0.0252	0.0212	0.0200
	200	0.0090	0.0036	0.0024	0.0017	0.0279	0.0179	0.0156	0.0145

and

Table 3. Estimation results for DGP 3.

	Bias					RMSE
Parameter	$(\mathit{N},\mathit{T})$	50	100	150	200	50	100	150	200
$\varphi$	CCEMG estimation
	50	−0.0649	−0.0330	−0.0154	−0.0076	0.0733	0.0475	0.0363	0.0342
	100	−0.0760	−0.0370	−0.0159	−0.0119	0.0801	0.0440	0.0313	0.0259
	200	−0.0789	−0.0378	−0.0179	−0.0127	0.0814	0.0434	0.0284	0.0222
	Jackknife bias−corrected CCEMG estimation
	50	−0.0195	−0.0073	0.0041	0.0010	0.0425	0.0368	0.0340	0.0335
	100	−0.0246	−0.0098	−0.0030	0.0001	0.0410	0.0272	0.0252	0.0235
	200	−0.0309	−0.0091	−0.0059	−0.0026	0.0395	0.0224	0.0182	0.0171
${\beta }_0$	CCEMG estimation
	50	0.0094	0.0063	0.0043	0.0045	0.0412	0.0329	0.0299	0.0276
	100	0.0092	0.0060	0.0045	0.0039	0.0312	0.0237	0.0205	0.0198
	200	0.0092	0.0062	0.0043	0.0038	0.0224	0.0174	0.0148	0.0141
	Jackknife bias-corrected CCEMG estimation
	50	0.0069	0.0045	0.0032	0.0039	0.0441	0.0341	0.0307	0.0284
	100	0.0068	0.0041	0.0035	0.0027	0.0330	0.0243	0.0212	0.0202
	200	0.0058	0.0040	0.0029	0.0021	0.0230	0.0175	0.0148	0.0143

for DGPs 1–3, respectively.

From Table 1, we note that for the estimation of $\varphi$, the CCEMG performs well in terms of bias and RMSE, with the bias diminishing as T is increased, and the associated RMSEs fall steadily when T increases, which implies that the CCEMG estimator is consistent. However, it still suffers from the time series bias when T is small. While the Jackknife bias-corrected CCEMG estimator is quite effective at reducing the time series bias of the CCEMG estimator, the bias has been significantly reduced compared with the original CCEMG estimator when T was not large, and the RMSE also decreased with the increase of either N or T. Similar findings can be observed for ${\beta }_0.$

In order to evaluate the robustness of various estimators, we considered additional results in Table 2 and Table 3 for DGPs with both stationary and non-stationary factors or cointegrated factors. Similar to the case with non-stationary factors in DGP1, we find that the CCEMG estimator still performs well regardless of the number of common factors and the non-stationary type, and it can be improved by the Jackknife bias-corrected for the estimation of the autoregressive coefficient $\varphi$, the CCEMG estimator of the slope coefficient ${\beta }_0$ performs very well in almost all cases.

Overall, the findings of our Monte Carlo simulations show that, if the parameter of interest is the mean coefficient of the regressors, ${\beta }_0$, the CCEMG estimator performs well even if N and T are not large. For the mean coefficient of the lagged dependent, $\varphi$, the CCEMG estimator is still consistent, but it suffers from the time series bias unless T is sufficiently large and, thus, the Jackknife bias-corrected CCEMG estimator is proposed, it helps to mitigate the time series bias.

5. Empirical Study

In this section, we illustrate our method by considering the U.S. Cigar dataset, which is frequently used in the literature on panel models (e.g., Baltagi and Li 2004; Bada and Liebl 2014). The panel contains the per capita cigarette consumption of $N=46$ American states from 1963 to 1992 ($T=30$) as well as data on the income per capita and cigarette prices; the dataset can be obtained from the R package phtt.

To test the cross-sectional dependence in the panel data, following Pesaran (2015) and Bailey et al. (2016), we compute the $CD$ statistic and the $\alpha$ statistic for the variables of interest in

Table 4. Exponent of the cross-sectional dependence of variables.

Variable	$\mathit{CD}$	p-Value	$\mathit{\alpha }$	${\mathit{\alpha }}_{0.025}$	${\mathit{\alpha }}_{0.975}$
consumption	101.519	0.000	0.975	0.887	1.064
income	166.270	0.000	1.004	−0.635	2.644
price	154.142	0.000	1.004	0.620	1.389

. As can be seen from the table, the $CD$ statistics turn out to be $101.519$, $166.27$, and $154.142$ for consumption, income, and price, respectively; these are highly significant and reject the null hypothesis of weak cross-sectional dependence for all three variables. Additionally, the estimates of $\alpha$ together with their $95\%$ confidence bands further confirm the above results. As a result, we can conclude that there is an obvious cross-sectional dependence for these three variables.

To investigate the relationship between the per capita cigarette consumption and the income per capita as well as cigarette prices, following Baltagi and Li (2004), we consider the panel model

$$y_{it}=c_i+{\varphi }_i{y}_{i,t-1}+{\beta }_{1i}{x}_{1it}+{\beta }_{2i}{x}_{2it}+e_{it},$$

(22)

where $y_{it},$$x_{1it}$, and $x_{2it}$ denote the per capita cigarette consumption, the income per capita, and cigarette price for the ith state at time t, respectively, and the idiosyncratic error has the multi-factor structure

$$e_{it}={\gamma }_i^{\prime }{\mathbf{f}}_t+{\epsilon }_{it}.$$

(23)

The proposed dynamic CCE approach is applied to estimate the coefficients in model (22), and the augmented equation to be estimated can be written as

$$y_{it}=c_i+{\varphi }_i{y}_{i,t-1}+{\beta }_{1i}{x}_{1it}+{\beta }_{2i}{x}_{2it}+\sum _{l=0}^{p_T}{\delta }_{il}^{\prime }{\overline{\mathbf{z}}}_{t-l}+w_{it},$$

(24)

where the number of lags $p_T=\left[\sqrt[3]{T}\right]=3$, and ${\overline{\mathbf{z}}}_t={\left({\overline{y}}_{t-1},{\overline{x}}_{1t},{\overline{x}}_{2t}\right)}^{\prime }$. We focus on the CCEMG estimators and the results are presented in

Table 5. Estimation results (Jackknife bias-corrected CCEMG).

Variable	coef.	Std.Err.	p-Value	${\mathit{CI}}_{0.025}$	${\mathit{CI}}_{0.975}$
$L.$ consumption	0.368	0.091	0.000	0.190	0.545
income	0.936	0.387	0.016	0.177	1.695
price	−0.629	0.115	0.000	−0.854	−0.404

Note: L. consumption denotes the first lag of cigarette consumption.

.

The following conclusions can be drawn from Table 5. On the one hand, the income per capita has a positive effect on the per capita cigarette consumption, while the increase in cigarette price will restrain cigarette consumption to a certain extent, and both are significant. These results are consistent with the conclusions of Bada and Liebl (2014). On the other hand, the lagged explained variable is highly significant, indicating that it is appropriate to use dynamic models for the per capita cigarette consumption.

To illustrate the heterogeneous slopes across states, we display both the CCE and the CCEMG estimators in Figure 1, which clearly show that the estimates of coefficients vary from state to state, reflecting the heterogeneity among states. Moreover, to illustrate the potential non-stationarity of unobservable common factors in (23), we consider the method proposed by Bada and Kneip (2014) to select the number of unobservable common factors and estimate the selected common factors. The results are given in Figure 2, where the top panel shows the estimated common factors and the bottom panel shows the estimated time-varying individual effects of $N=46$ states. As can be seen from the figure, five common factors have been selected, among which the first and second common factors have obvious tendencies and violate the stationarity condition.

6. Conclusions

In this paper, we re-examined the CCE type estimator for dynamic heterogeneous panel regression models with non-stationary common factors. Asymptotic properties of CCE estimators are established when both N and T are large. It is shown that, under certain conditions, the main results of Pesaran (2006) and Chudik and Pesaran (2015a) hold for a dynamic panel with non-stationary factors. Monte Carlo simulations were conducted to investigate the finite sample properties of the CCE estimation for the panel with non-stationary factors. An empirical application to the U.S. cigarette consumption dataset shows that the real data may have cross-sectional dependence as well as dynamic and non-stationary common factors (at the same time). Based on the findings of this paper, together with the results by Pesaran (2006); Kapetanios et al. (2011), and Chudik and Pesaran (2015a), we can conclude that the CCE method can be widely used to deal with panel models with error cross-sectional dependence, regardless of whether the model is static or dynamic, and whether the unobservable common factors are stationary.