A Robust Approach to Heteroscedasticity, Error Serial Correlation and Slope Heterogeneity in Linear Models with Interactive Effects for Large Panel Data

Abstract

In this article, we propose a robust approach against heteroscedasticity, error serial correlation and slope heterogeneity in linear models with interactive effects for large panel data. First, consistency and asymptotic normality of the pooled iterated principal component (IPC) estimator for random coefficient and homogeneous slope models are established. Then, we prove the asymptotic validity of the associated Wald test for slope parameter restrictions based on the panel heteroscedasticity and autocorrelation consistent (PHAC) variance matrix estimator for both random coefficient and homogeneous slope models, which does not require the Newey-West type time-series parameter truncation. These results asymptotically justify the use of the same pooled IPC estimator and the PHAC standard error for both homogeneous-slope and heterogeneous-slope models. This robust approach can significantly reduce the model selection uncertainty for applied researchers. In addition, we propose a Lagrange Multiplier (LM) test for correlated random coefficients with covariates. This test has nontrivial power against correlated random coefficients, but not for random coefficients and homogeneous slopes. The LM test is important because the IPC estimator becomes inconsistent with correlated random coefficients. The finite sample evidence and an empirical application support the reliability and the usefulness of our robust approach.

Keywords:

• Correlated random coefficients
• Interactive effects
• Panel data
• Slope heterogeneity

1 Introduction

In recent years, the increasing availability of panel data with both a large cross-sectional dimension N and a large time series dimension T has necessitated the development of statistical methods to exploit the rich information they contain, but it has also created technical challenges. In particular, the control of cross-sectional dependence, parameter heterogeneity and serial dependence has been the main focus of the literature. Among other things, the modeling of the cross-sectional dependence by interactive effects has been widely adopted in recent years.

In particular, the past decade has seen significant development in estimation and inferential methods for linear models with interactive effects for large panel data. Two estimation approaches are particularly popular. The first involves eliminating the interactive effects from both the error term and the regressors. Representative methods for static models include the Common Correlated Effects (CCE) approach proposed by Pesaran (2006), the principal component (PC) approach investigated by Westerlund and Urbain (2015), the maximum likelihood (ML) approach of Bai and Li (2014), and the two-step instrumental variable (2SIV) approach proposed by Cui et al. (2021), among others. The second approach asymptotically eliminates the interactive effects from the error term only. The representative method is the iterative principal components (IPC) estimator of Bai (2009), further developed by Moon and Weidner (2015), and Bai and Liao (2017), among many others.

Concerning the $\sqrt{NT}$-consistent estimation and the associated inference against error serial correlation and heteroscedasticity, the CCE, PC, 2SIV, and IPC approaches typically permit idiosyncratic errors to be serially correlated and heteroscedastic (see, Remark 1 in Westerlund and Urbain 2015, Assumptions 2.1 and 2.2 in Cui et al. 2021, Remark 8 in Bai 2009), whereas the ML approach generally does not allow such a degree of robustness (see Assumption B in Bai and Li 2014). The CCE, PC, and IPC estimators, however, will have asymptotic biases, essentially due to the endogeneity arising from estimating the interactive effects. Note that the properties of the asymptotic biases in the CCE and PC estimators and those of the IPC estimator are different, because of the difference in their approaches for eliminating the interactive effects. These asymptotic biases should be controlled, particularly for inference, either analytically or numerically by, for example, jackknife subsampling. For the latter, see Fernández and Weidner (2016) and Westerlund (2018).

For testing linear restrictions on the slope coefficients, most existing work routinely employs the Newey and West (1987) type heteroscedasticity and autocorrelation consistent (HAC) variance-covariance estimator, which requires time-series parameter truncation. This includes the CCE and the PC estimators (see Westerlund 2018), as well as the IPC estimator (see Remark 8 in Bai 2009). The 2SIV estimator of Cui et al. (2021) and other related papers, including Greenaway-McGrevy, Han, and Sul (2012), Bonhomme and Manresa (2015), and Fernández and Weidner (2016), employ the Hansen (2007) type panel HAC (PHAC) variance-covariance estimator without the parameter truncation. In this article we employ this PHAC estimator.¹

One of the attractive features of the CCE, PC, and 2SIV approaches is that they permit estimation and inference for the models with (cross-sectionally) heterogeneous slopes. For such models, the parameter of interest is often the population average of slope coefficients, which is typically estimated as an average of estimators of cross-section specific slopes. This is called a mean-group (MG) estimator. For the MG estimators based on the CCE, PC, and 2SIV approaches, see Pesaran (2006) and Cui et al. (2021). Interestingly, as shown by Pesaran (2006) and Reese and Westerlund (2018), if the cross-sectional variation of slopes is independent of the regressors (i.e., if slopes are random coefficients), the aforementioned pooled CCE and PC estimators are $\sqrt{N}$-consistent to the population average of slope coefficients. Pesaran (2006) proposes a nonparametric estimator of the variance-covariance matrix of the $\sqrt{N}$-consistent pooled CCE estimator. This estimator uses the sample variance-covariance matrix of cross-section specific slope estimates. The evidence therein has shown that this estimator behaves very well in finite samples.²

On the other hand, the IPC approaches typically do not permit models with heterogeneous slopes.³ In this article we investigate asymptotic properties of the pooled IPC estimator for the models with heterogeneous slopes assuming the regressors are subject to a factor structure. To the best of our knowledge, this is new to the literature. In this article, we show that the pooled IPC estimator is $\sqrt{N}$-consistent for the models with heterogeneous slopes (and confirm it is $\sqrt{NT}$-consistent for the models with homogeneous slopes), when the regressors have the factor structure. We also prove that the use of the PHAC variance-covariance matrix estimator for the pooled IPC estimator is asymptotically justified for the models with heterogeneous slopes, as well as for the models with homogeneous slopes. Our findings essentially imply that the inference based on the pooled IPC estimator and the associated PHAC variance estimator is asymptotically valid for both homogeneous slope and heterogeneous slope models.

Building upon these novel asymptotic results, we propose a robust approach against heteroscedasticity, error serial correlation and slope heterogeneity in linear models with interactive effects for large panel data. The basic idea behind this approach is classical. For the models with homogeneous slopes and homoscedastic and uncorrelated idiosyncratic errors, the pooled estimator is expected to be most efficient. Otherwise, an alternative estimator, such as a suitable generalized least squares (GLS) estimator, can be more efficient; see Bai and Liao (2017) for the extension of the IPC estimator in this direction. This potential efficiency gain hinges on the correct specification of, in our context, error serial correlation, heteroscedasticity and slope heterogeneity, and achieving such specifications simultaneously does not seem easy in practice. In this kind of situation, the use of the pooled estimator, together with a robust variance-covariance matrix estimator for inference, is very likely to have practical appeal.⁴

Another main contribution of this article is a novel test for correlation and dependency between the random coefficients and covariates (hereafter called correlated random coefficients). For the proposed robust approach to work, the slope heterogeneity, if any, has to be essentially independent of the covariates (i.e., random coefficients). The proposed test is designed to detect departures from this independence and slope homogeneity. Note that the existing slope homogeneity tests for large panel data, such as Pesaran and Yamagata (2008) and Su and Chen (2013), among others, have power against random coefficients as well as correlated random coefficients, and these are not suitable for our purpose.

Here, we briefly discuss a related study, Galvao and Kato (2014), which considers estimation and inference of pseudo-true parameters under some model misspecifications, such as inconsistent estimators in the fixed effect models with correlated random coefficients (see sec. 5.1.2 in Galvao and Kato 2014). Their theoretical results are closely related to ours in that the convergence rate of the pooled estimator changes from $\sqrt{NT}$ to $\sqrt{N}$ under correlated random coefficients. It is important to emphasize that this article deals with consistent estimation of true parameters of interest and the associated inference in the models with interactive effects, which is not explicitly considered in Galvao and Kato (2014). In addition, in contrast to Galvao and Kato (2014), our approach advocates a method that enables empirical researchers to statistically identify misspecified models (i.e., correlated random coefficient models) that can result in inconsistent estimation.

We have examined the finite sample performance of the estimators, the tests of parameter restrictions using the PHAC variance estimator, and the LM test for correlated random coefficients. The results show that the size of the proposed robust Wald test with the bias-corrected IPC estimator is sufficiently close to the nominal level in both slope homogeneity and slope heterogeneity, and that the LM test for correlated random coefficients has correct size under both slope homogeneity and random coefficients while exhibiting high power for correlated random coefficients. The practical usefulness of our robust approach is illustrated by applying the proposed methods to analyze the Feldstein–Horioka puzzle, presented in Feldstein and Horioka (1980).

The article is organized as follows. The asymptotic properties of the IPC estimator and the associated robust Wald test statistic based on the PHAC variance-covariance estimator are investigated in Section 2. In Section 3, we introduce the LM test statistic for correlated random coefficients and study its limiting distribution. The finite sample performance of the IPC estimators, the robust Wald test and the LM test for correlated random coefficients is investigated using the Monte Carlo method in Section 4. An empirical application is provided in Section 5. Section 6 contains some concluding remarks. Proofs of the main results, additional discussions and experimental results are relegated to the supplementary materials.

Notations: we denote the largest eigenvalues of the N × N matrix $\mathbf{A}=(a_{ij})$ by ${\mu }_{\text{max}}(\mathbf{A})$, its trace by $\text{tr}(\mathbf{A})={\sum }_{i=1}^N{a}_{ii}$, its Frobenius norm by $\left|\right|\mathbf{A}\left|\right|=\sqrt{\text{tr}({\mathbf{A}}^{\prime }\mathbf{A})}$. The projection matrix on A is ${\mathbf{P}}_{\mathbf{A}}=\mathbf{A}{({\mathbf{A}}^{\prime }\mathbf{A})}^{-1}{\mathbf{A}}^{\prime }$ and ${\mathbf{M}}_{\mathbf{A}}=\mathbf{I}-{\mathbf{P}}_{\mathbf{A}}$. Δ is a generic positive constant large enough, ${\delta }_{N\mathrm{}T}^2=\min \{N,T\}$. We use $N,T\to \infty$ to denote that N and T pass to infinity jointly. We use $\lesssim$ ($\gtrsim$) to represent $\le$ ($\ge$) up to a positive constant factor. For any positive sequences a_n and b_n that converge to some points or diverge as $n\to \infty$, we write $a_n\asymp b_n$ if $a_n\lesssim b_n$ and $a_n\gtrsim b_n$.

2 Panel Data Models with Interactive Effects

We consider the following panel data model with possibly heterogeneous coefficients,(1) $$\begin{array}{c}{y}_{it}={\mathbf{x}}_{it}^{\prime }{\mathit{\beta }}_i+{\mathit{\lambda }}_i^{0\prime }{\mathbf{f}}_t^0+{\epsilon }_{it},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}(i=1,2,\dots ,N ;t=1,2,\dots ,T ),\end{array}$$(1) (2) $${\mathbf{x}}_{it}={\mathbf{\Gamma }}_i^{0\prime }{\mathbf{g}}_t^0+{\mathbf{v}}_{it},$$(2) where ${\mathbf{x}}_{it}$ is a $k\times 1$ vector of regressors, ${\mathbf{f}}_t^0$ and ${\mathbf{g}}_t^0$ denote $r_1\times 1$ and $r_2\times 1$ vectors of latent factors, respectively. Correspondingly, their factor loadings are ${\mathit{\lambda }}_i^0$ and ${\mathbf{\Gamma }}_i^0=({\mathit{\gamma }}_{1i}^0,\dots ,{\mathit{\gamma }}_{ki}^0)$. ${\epsilon }_{it}$ and ${\mathbf{v}}_{it}$ are the idiosyncratic disturbance terms.

If we stack the Equations (1)(1) $$\begin{array}{c}{y}_{it}={\mathbf{x}}_{it}^{\prime }{\mathit{\beta }}_i+{\mathit{\lambda }}_i^{0\prime }{\mathbf{f}}_t^0+{\epsilon }_{it},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}(i=1,2,\dots ,N ;t=1,2,\dots ,T ),\end{array}$$(1) and (2)(2) $${\mathbf{x}}_{it}={\mathbf{\Gamma }}_i^{0\prime }{\mathbf{g}}_t^0+{\mathbf{v}}_{it},$$(2) over t, we have(3) $${\mathbf{y}}_i={\mathbf{X}}_i{\mathit{\beta }}_i+{\mathbf{F}}^0{\mathit{\lambda }}_i^0+{\mathit{\epsilon }}_i ,$$(3) (4) $${\mathbf{X}}_i={\mathbf{G}}^0{\mathbf{\Gamma }}_i^0+{\mathbf{V}}_i,$$(4) where ${\mathbf{y}}_i={(y_{i1},\dots ,y_{iT})}^{\prime },{\mathbf{X}}_i={({\mathbf{x}}_{i1},\dots ,{\mathbf{x}}_{iT})}^{\prime },{\mathbf{F}}^0={({\mathbf{f}}_1^0,\dots ,{\mathbf{f}}_T^0)}^{\prime },{\mathbf{G}}^0={({\mathbf{g}}_1^0,\dots ,{\mathbf{g}}_T^0)}^{\prime }$, ${\mathit{\epsilon }}_i={({\epsilon }_{i1},\dots ,{\epsilon }_{iT})}^{\prime }$ and ${\mathbf{V}}_i={({\mathbf{v}}_{i1},\dots ,{\mathbf{v}}_{iT})}^{\prime }$. Supposing random coefficients(5) $${\mathit{\beta }}_i={\mathit{\beta }}^0+{\mathit{\eta }}_i,$$(5) we can rewrite the model (3) using the process (4) as(6) $${\mathbf{y}}_i={\mathbf{X}}_i{\mathit{\beta }}^0+{\mathbf{W}}^0{\mathit{\vartheta }}_i^0+{\mathbf{V}}_i{\mathit{\eta }}_i+{\mathit{\epsilon }}_i,$$(6) where ${\mathbf{W}}^0=({\mathbf{G}}^0,{\mathbf{F}}^0)$ and ${\mathit{\vartheta }}_i^0={({\mathit{\eta }}_i^{\prime }{\mathbf{\Gamma }}_i^{0\prime },{\mathit{\lambda }}_i^{0\prime })}^{\prime }$. The results in this article will hold when ${\mathbf{F}}^0$ and ${\mathbf{G}}^0$ share common factors (say, ${\mathbf{W}}^0={\mathbf{G}}^0\cup {\mathbf{F}}^0$ and ${\mathbf{F}}^0\cap {\mathbf{G}}^0\ne \mathbf{0}$). However, for simplicity, we suppose ${\mathbf{F}}^0\cap {\mathbf{G}}^0=\mathbf{0}$ and $\text{rank}({\mathbf{W}}^0)=r_1+r_2$.

Noting that for the IPC estimator the interactive effects ${\mathbf{F}}^0{\mathit{\lambda }}_i^0$ for the models with homogeneous slopes and ${\mathbf{W}}^0{\vartheta }_i^0$ for the models with heterogeneous slopes are extracted from the error vector ${\mathbf{y}}_i-{\mathbf{X}}_i{\mathit{\beta }}_i^0$, we employ the following common expression to describe both homogeneous-slope and heterogeneous-slope models:(7) $${\mathbf{y}}_i={\mathbf{X}}_i{\mathit{\beta }}^0+{\mathbf{u}}_i,\mathrm{ }\text{with}\mathrm{ }{\mathbf{u}}_i={\mathbf{H}}^0{\mathit{\varphi }}_i^0+{\mathbf{e}}_i,$$(7) where $\{{\mathbf{H}}^0,{\mathit{\varphi }}_i^0,{\mathbf{e}}_i\}=\{{\mathbf{F}}^0,{\mathit{\lambda }}_i^0,{\mathit{\epsilon }}_i\}$ and $r=r_1$ for the models with homogeneous slopes and $\{{\mathbf{H}}^0,{\mathit{\varphi }}_i^0,{\mathbf{e}}_i\}=\{{\mathbf{W}}^0,{\mathit{\vartheta }}_i^0,{\mathbf{V}}_i{\mathit{\eta }}_i+{\mathit{\epsilon }}_i\}$ and $r=r_1+r_2$ for the models with heterogeneous slopes. Hereafter, we use Equation (7)(7) $${\mathbf{y}}_i={\mathbf{X}}_i{\mathit{\beta }}^0+{\mathbf{u}}_i,\mathrm{ }\text{with}\mathrm{ }{\mathbf{u}}_i={\mathbf{H}}^0{\mathit{\varphi }}_i^0+{\mathbf{e}}_i,$$(7) to describe both homogeneous-slope and homogeneous-slope models.

Controlling the interactive effects, ${\mathbf{H}}^0{\mathit{\varphi }}_i^0$, typically results in asymptotic biases of the $\sqrt{NT}$-consistent CCE, PC, and IPC estimators, as shown by Pesaran (2006), Bai (2009), and Westerlund and Urbain (2015), among others.

As discussed in the introduction section, in this article we focus on the iterative principal component (IPC) estimator proposed by Bai (2009), which estimates the interactive effects in ${\mathbf{u}}_i$ by the PC method and the slope coefficient ${\mathit{\beta }}^0$ by the least-squares method, iteratively.

In particular, define the least squares objective function$$\text{SSR}(\mathit{\beta },\mathbf{H})=\frac{1}{NT}\sum _{i=1}^N{({\mathbf{y}}_i-{\mathbf{X}}_i\mathit{\beta }-\mathbf{H}{\mathit{\varphi }}_i)}^{\prime }({\mathbf{y}}_i-{\mathbf{X}}_i\mathit{\beta }-\mathbf{H}{\mathit{\varphi }}_i)$$subject to the constraints ${\mathbf{H}}^{\prime }\mathbf{H}/T={\mathbf{I}}_r$ and ${\sum }_{i=1}^N{\mathit{\varphi }}_i{\mathit{\varphi }}_i^{\prime }$ being diagonal.

The least squares estimator for $(\mathit{\beta },\mathbf{H})$, denoted by $(\widehat{\mathit{\beta }},\widehat{\mathbf{H}})$, is the IPC estimator, which is the solution to the following nonlinear equations:(8) $$\begin{array}{c}\widehat{\mathit{\beta }}={\left(\sum _{i=1}^N{\mathbf{X}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\mathbf{X}}_i\right)}^{-1}\left(\sum _{i=1}^N{\mathbf{X}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\mathbf{y}}_i\right)\\ \left[\frac{1}{NT}\sum _{i=1}^N\left({\mathbf{y}}_i-{\mathbf{X}}_i\widehat{\mathit{\beta }}\right){\left({\mathbf{y}}_i-{\mathbf{X}}_i\widehat{\mathit{\beta }}\right)}^{\prime }\right]\widehat{\mathbf{H}}=\widehat{\mathbf{H}}{\mathcal{V}}_{NT}\end{array}$$(8) where ${\mathcal{V}}_{NT}$ is a diagonal matrix that contains the r largest eigenvalues of the matrix in the square brackets in decreasing order. Given $(\widehat{\mathit{\beta }},\widehat{\mathbf{H}}),$ we can estimate ${\mathit{\varphi }}_i$ by$${\widehat{\mathit{\varphi }}}_i=\frac{1}{T}{\widehat{\mathbf{H}}}^{\prime }({\mathbf{y}}_i-{\mathbf{X}}_i\widehat{\mathit{\beta }}).$$

The number of factors, r, and an initial estimator of ${\mathit{\beta }}^0$, should be chosen to obtain the IPC estimator. In this section and the next, we suppose that r and the initial estimator of ${\mathit{\beta }}^0$ are given, and the procedure to choose them is discussed in Section 4.1.

Remark 1.

Permitting different sets of interactive effects in ${\mathbf{X}}_i$ and ${\mathbf{u}}_i$ is important, because the IPC estimator extracts the interactive effects in ${\mathbf{u}}_i$ only, which does not necessarily project out the factors in ${\mathbf{X}}_i$. Nevertheless, our results in this article hold when ${\mathbf{X}}_i$ and ${\mathbf{u}}_i$ share the same set of interactive effects. See Norkutė et al. (2021) and Cui et al. (2021) for further discussions.

Remark 2.

The factor model for the regressors specified in (2) has been widely employed in the literature: see Pesaran (2006), Bai and Li (2014), Westerlund and Urbain (2015), Li, Cui, and Lu (2020), Norkutė et al. (2021), among many others. For the IPC estimators, typically no factor structure in the regressors is imposed. However, only the model with homogeneous slopes is considered; see Bai (2009) and Moon and Weidner (2015), among others. Importantly, our new results shown below will reveal that, the pooled IPC estimator remains consistent for the random coefficient models with interactive effects if the regressors have a factor structure as in (2).

For the asymptotic analysis of the estimator, we impose the following assumptions.

Assumption 1

(idiosyncratic error in y). (i) ${\epsilon }_{it}$ is distributed independently across i; (ii) $\mathbb{E}({\epsilon }_{it})=0$ and $\mathbb{E}|{\epsilon }_{it}{|}^8\le \Delta$; (iii) $T^{-1}{\sum }_{s=1}^T{\sum }_{t=1}^T\mathbb{E}\left|{\epsilon }_{is}{\epsilon }_{it}\right|\le \Delta$; (iv) $\mathbb{E}|N^{-1/2}{\sum }_{i=1}^N$ $\left[{\epsilon }_{is}{\epsilon }_{it}-\mathbb{E}({\epsilon }_{is}{\epsilon }_{it})\right]{|}^4\le \Delta$ for each (s, t); (v) $N^{-1}{T}^{-2}{\sum }_{i=1}^N$ ${\sum }_{t=1}^T{\sum }_{s=1}^T{\sum }_{r=1}^T{\sum }_{w=1}^T\left|\text{cov}({\epsilon }_{is}{\epsilon }_{it},{\epsilon }_{ir}{\epsilon }_{iw})\right|\le \Delta$; (vi) ${\mathbf{\Omega }}_{\epsilon ,i}=\mathbb{E}({\mathit{\epsilon }}_i{\mathit{\epsilon }}_i^{\prime })$ is positive definite and its largest eigenvalue is bounded, uniformly every i and T.

Assumption 2

(idiosyncratic error in x). Let $v_{\mathcal{l}it}$ be the $\mathcal{l}$th element of ${\mathbf{v}}_{it}$ and ${\mathbf{v}}_{\mathcal{l},i}={(v_{\mathcal{l}i1},\dots ,v_{\mathcal{l}iT})}^{\prime }$. Then we assume that (i) ${\mathbf{v}}_{it}$ is independently distributed across i and independent of $\left\{{\epsilon }_{js}\right\}$ for $1\le j\le N$ and $1\le s\le T$; (ii) $\mathbb{E}({\mathbf{v}}_{it})=\mathbf{0}$ and $\mathbb{E}\left|\right|{\mathbf{v}}_{it}|{|}^8\le \Delta$; (iii) $T^{-1}{\sum }_{s=1}^T{\sum }_{t=1}^T\mathbb{E}\left|v_{\mathcal{l}is}{v}_{\mathcal{l}it}\right|\le \Delta$; (iv) $\mathbb{E}|N^{-1/2}{\sum }_{i=1}^N\left[v_{\mathcal{l}is}{v}_{\mathcal{l}it}-\mathbb{E}(v_{\mathcal{l}is}{v}_{\mathcal{l}it})\right]{|}^4\le \Delta$ for every $\mathcal{l}$, s, and t; (v) $N^{-1}{T}^{-2}{\sum }_{i=1}^N{\sum }_{t=1}^T{\sum }_{s=1}^T{\sum }_{r=1}^T{\sum }_{w=1}^T$ $|\text{cov}(v_{\mathcal{l}is}{v}_{\mathcal{l}it},v_{\mathcal{l}ir}{v}_{\mathcal{l}iw})|\le \Delta$; (vi) the largest eigenvalue of $\mathbb{E}({\mathbf{v}}_{\mathcal{l},i}{\mathbf{v}}_{\mathcal{l},i}^{\prime })$ is bounded uniformly for every $\mathcal{l}$, i, and T.

Assumption 3

(random coefficients). $$ is independent across i, and is independent of ${\epsilon }_{it}$, ${\mathbf{v}}_{it},{\mathit{\lambda }}_i^0$ and ${\mathbf{\Gamma }}_i^0$; $\mathbb{E}({\mathit{\eta }}_i)=\mathbf{0}$ and $\mathbb{E}\left|\right|{\mathit{\eta }}_i|{|}^4\le \Delta$; $\mathbb{E}({\mathit{\eta }}_i{\mathit{\eta }}_i^{\prime })={\mathbf{\Omega }}_{\eta ,i}$ is a fixed positive definite matrix uniformly over i.

Assumption 4

(factor components). (i) Let ${\mathbf{w}}_t^0={({\mathbf{f}}_t^{0\prime },{\mathbf{g}}_t^{0\prime })}^{\prime }$. $\mathbb{E}\left|\right|{\mathbf{w}}_t^0|{|}^4\le \Delta$ and $T^{-1}{\sum }_{t=1}^T{\mathbf{w}}_t^0{\mathbf{w}}_t^{0\prime }\stackrel{p}{\to }{\mathbf{\Sigma }}_w=\left[\begin{array}{cc}{\mathbf{\Sigma }}_f& {\mathbf{\Sigma }}_{fg}\\ {\mathbf{\Sigma }}_{gf}& {\mathbf{\Sigma }}_g\end{array}\right]$, which is positive definite and ${\mathbf{w}}_t^0$ is independent of ${\mathbf{v}}_{it}$ and ${\epsilon }_{it}$; (ii) $\mathbb{E}\left|\right|{\mathbf{\Gamma }}_i^0|{|}^4\le \Delta$, $\mathbb{E}\left|\right|{\mathit{\lambda }}_i^0|{|}^4\le \Delta$, $N^{-1}{\sum }_{i=1}^N{\mathbf{\Gamma }}_i^0{\mathbf{\Gamma }}_i^{0\prime }\stackrel{p}{\to }{\Sigma }_{\Gamma },N^{-1}{\sum }_{i=1}^N{\mathbf{\Gamma }}_i^0{\mathit{\eta }}_i{\mathit{\eta }}_i^{\prime }{\mathbf{\Gamma }}_i^{0\prime }\stackrel{p}{\to }{\mathbf{\Sigma }}_{\gamma \eta }$ and $N^{-1}{\sum }_{i=1}^N{\mathit{\lambda }}_i^0{\mathit{\lambda }}_i^{0\prime }\stackrel{p}{\to }{\mathbf{\Sigma }}_{\lambda }$, where ${\mathbf{\Sigma }}_{\Gamma },{\mathbf{\Sigma }}_{\gamma \eta }$, and ${\mathbf{\Sigma }}_{\lambda }$ are positive definite matrices. ${\mathbf{\Gamma }}_i^0$ and ${\mathit{\lambda }}_i^0$ are independent of ${\mathbf{v}}_{it}$ and ${\epsilon }_{it}$.

Assumption 5

(identification and variance matrices). We assume that the four matrices$$\begin{array}{c}{\mathbf{A}}_0=\underset{N,T\to \infty }{\lim }\frac{1}{NT}\sum _{i=1}^N\mathbb{E}({\mathbf{V}}_i^{\prime }{\mathbf{V}}_i),\\ {\mathbf{A}}_1=\underset{N,T\to \infty }{\text{plim}}\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i,\\ {\mathbf{C}}_0=\underset{N,T\to \infty }{\lim }\frac{1}{N{T}^2}\sum _{i=1}^N\mathbb{E}({\mathbf{V}}_i^{\prime }{\mathbf{V}}_i{\mathit{\eta }}_i{\mathit{\eta }}_i^{\prime }{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i),\\ {\mathbf{B}}_1=\underset{N,T\to \infty }{\text{plim}}\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{\Omega }}_{\epsilon ,i}{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\end{array}$$are fixed and positive definite, where ${\mathbf{Z}}_i={\mathbf{X}}_i-N^{-1}{\sum }_{j=1}^N{\mathit{\varphi }}_i^{0\prime }$ ${\mathbf{\Upsilon }}_{{\varphi }^0}^{-1}{\mathit{\varphi }}_j^0{\mathbf{X}}_j$, ${\mathbf{\Upsilon }}_{{\varphi }^0}=N^{-1}{\sum }_{i=1}^N{\mathit{\varphi }}_i^0{\mathit{\varphi }}_i^{0\prime }$.

Assumption 6

(Central Limit Theorem). When $N,T\to \infty$such that $T/N\to \rho \in (0,\Delta ]$, $\frac{1}{\sqrt{NT}}{\sum }_{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathit{\epsilon }}_i\stackrel{d}{\to }N(\mathbf{0},{\mathbf{B}}_1)$.

Assumptions 1 and 2 are very similar to Assumption in Bai (2009), though, ${\epsilon }_{it}$ and ${\mathbf{v}}_{it}$ are independent of each other and independent over i. Assumption 3 is a typical assumption for random coefficient models with interactive effects (see, e.g., Assumption HET in Reese and Westerlund 2018). Assumption 4 implies that the factor loadings ${\mathbf{\Gamma }}_i^0$ and ${\mathit{\lambda }}_i^0$ are allowed to be correlated with each other. Assumption 5 is for identification of the parameters and for estimation of the variance-covariance matrices. Assumption 6 is imposed following Bai (2009) for the asymptotic normality of the pooled IPC estimator for the models with homogeneous slopes.

In what follows, we investigate asymptotic properties of the IPC estimator for the models with homogeneous and heterogeneous slopes, assuming that the regressors have a factor structure. As discussed in the introduction section, the results for the models with heterogeneous slopes are new additions to the literature. We first derive asymptotic expressions of the pooled IPC estimator, $\widehat{\mathit{\beta }}$, for the models with homogeneous and heterogeneous slopes in Proposition 1. These confirm the presence of asymptotic bias in $\widehat{\mathit{\beta }}$ for the models with homogeneous slopes, as shown by Bai (2009) and Cui et al. (2021), among others. After introducing an asymptotic expression for the bias-corrected IPC estimator, which is common to the models with homogeneous and heterogeneous slopes in Proposition 2, we present asymptotic normality of the bias-corrected IPC estimator for both models in Theorem 1. Finally, in Theorem 2, we show the asymptotic validity of the robust Wald test for linear constraints on ${\mathit{\beta }}^0$, which is based on the same statistic for both homogeneous-slope and heterogeneous-slope models.

A similar line of discussion in Bai (2009) and Cui et al. (2021) for the model with homogeneous slopes, and in Pesaran (2006) for the case of heterogeneous slopes provides the following asymptotic expressions:

Proposition 1.

Suppose that Assumptions 1–5 hold.

1. For the model with homogeneous slopes, we have(9) $$\begin{array}{c}\widehat{\mathit{\beta }}-{\mathit{\beta }}^0={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathit{\epsilon }}_i\\ + \frac{1}{N}{\mathit{\xi }}_{NT}+\frac{1}{T}{\mathit{\zeta }}_{NT}+O_p({\delta }_{N\mathrm{}T}^{-3}),\end{array}$$(9)

   where the bias terms are given by$$\begin{array}{c}{\mathit{\xi }}_{NT}=-{\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\\ \frac{1}{NT}\sum _{i=1}^N\sum _{t=1}^T{\mathbf{Z}}_i^{\prime }{\mathbf{H}}^0{({\mathbf{H}}^{0\prime }{\mathbf{H}}^0)}^{-1}{\mathbf{\Upsilon }}_{{\varphi }^0}^{-1}{\mathit{\varphi }}_i^0\mathbb{E}(e_{it}^2),\\ {\mathit{\zeta }}_{NT}=-{\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\\ \frac{1}{N^2}\sum _{i=1}^N\sum _{j=1}^N{\mathbf{X}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}\mathbb{E}({\mathbf{e}}_j{\mathbf{e}}_j^{\prime }){\mathbf{H}}^0{({\mathbf{H}}^{0\prime }{\mathbf{H}}^0)}^{-1}{\mathbf{\Upsilon }}_{{\varphi }^0}^{-1}{\mathit{\varphi }}_i^0.\end{array}$$

2. For the model with heterogeneous slopes, we have(10) $$\widehat{\mathit{\beta }}-{\mathit{\beta }}^0={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i\right)}^{-1}\frac{1}{NT}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i{\mathit{\eta }}_i+O_p({\delta }_{NT}^{-2}).$$(10)

A few remarks regarding Proposition 1 are warranted. Proposition 1(a) is essentially the same result for the models with the homogeneous slopes shown in Bai (2009) and Cui et al. (2021). Because of the estimation of the interactive effects, ${\mathbf{X}}_i$ has to be replaced by ${\mathbf{Z}}_i$, and the bias terms ${\mathit{\xi }}_{NT}$ and ${\mathit{\zeta }}_{NT}$ are present. Under our assumptions ${\mathit{\xi }}_{NT}$ and ${\mathit{\zeta }}_{NT}$ are $O_p(1)$, because ${\mathbf{Z}}_i$ can be correlated with ${\mathit{\varphi }}_i^0$ and $\mathbb{E}({\mathbf{e}}_j{\mathbf{e}}_j^{\prime })$ is not necessarily proportional to ${\mathbf{I}}_T$. Proposition 1(b) reveals that the pooled IPC estimator is $\sqrt{N}$-consistent and asymptotically free from the effects of interactive effects. See a similar result for pooled CCE and PC estimators in Pesaran (2006) and Reese and Westerlund (2018).

To derive an asymptotically unbiased estimator of ${\mathit{\beta }}^0$ in the case of homogeneous slope, we should estimate the bias terms. Following Bai (2009) and Moon and Weidner (2015), we propose an analytically bias corrected estimator, which is defined by(11) $$\tilde{\mathit{\beta }}=\widehat{\mathit{\beta }}-\frac{1}{N}{\widehat{\mathit{\xi }}}_{NT}-\frac{1}{T}{\widehat{\mathit{\zeta }}}_{NT},$$(11) where$$\begin{array}{c}{\widehat{\mathit{\xi }}}_{NT}=-{\left(\frac{1}{NT}\sum _{i=1}^N{\widehat{\mathbf{Z}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\widehat{\mathbf{Z}}}_i\right)}^{-1}\frac{1}{N{T}^2}\sum _{i=1}^N\sum _{t=1}^T{\widehat{\mathbf{Z}}}_i^{\prime }\widehat{\mathbf{H}}{\widehat{\mathbf{\Upsilon }}}_{\varphi }^{-1}{\widehat{\mathit{\varphi }}}_i{\widehat{e}}_{it}^2\\ {\widehat{\mathit{\zeta }}}_{NT}=-\left(\frac{1}{NT}\sum _{i=1}^N{\widehat{\mathbf{Z}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\widehat{\mathbf{Z}}}_i\right)\frac{1}{N}\sum _{i=1}^N\frac{\widehat{{\mathbf{X}}_i^{\prime }{\mathbf{M}}_{\mathbf{H}}\mathbf{\Omega }\mathbf{H}}}{T}{\widehat{\mathbf{\Upsilon }}}_{\varphi }^{-1}{\widehat{\mathit{\varphi }}}_i\end{array}$$ with ${\widehat{\mathbf{Z}}}_i={\mathbf{X}}_i-N^{-1}{\sum }_{j=1}^N{\widehat{\mathit{\varphi }}}_i^{\prime }{\widehat{\mathbf{\Upsilon }}}_{\varphi }^{-1}{\widehat{\mathit{\varphi }}}_j{\mathbf{X}}_j,{\widehat{\mathbf{\Upsilon }}}_{\varphi }=N^{-1}{\sum }_{i=1}^N{\widehat{\mathit{\varphi }}}_i{\widehat{\mathit{\varphi }}}_i^{\prime }$,(12) $$\begin{array}{c}\frac{\widehat{{\mathbf{X}}_i^{\prime }{\mathbf{M}}_{\mathbf{H}}\mathbf{\Omega }\mathbf{H}}}{T}=\frac{1}{TN}\sum _{j=1}^N[\sum _{t=1}^T{\widehat{e}}_{jt}^2{\widehat{\mathbf{x}}}_{it}{\widehat{\mathbf{h}}}_t^{\prime }+\sum _{s=1}^S\sum _{t=s+1}^T\left(1-\frac{s}{S+1}\right)\\ {\widehat{e}}_{jt}{\widehat{e}}_{j,t-s}\left({\widehat{\mathbf{x}}}_{it}{\widehat{\mathbf{h}}}_{t-s}^{\prime }+{\widehat{\mathbf{x}}}_{i,t-s}{\widehat{\mathbf{h}}}_t^{\prime }\right)]\end{array}$$(12) and ${\widehat{e}}_{it}$ being the tth element of ${\widehat{\mathbf{e}}}_i={\mathbf{y}}_i-{\mathbf{X}}_i\widehat{\mathit{\beta }}-\widehat{\mathbf{H}}{\widehat{\mathit{\varphi }}}_i$.

Following Bai’s (2009) suggestion we employ the Newey and West (1987) type estimator of ${\mathit{\zeta }}_{NT}$ as shown in Equation (12)(12) $$\begin{array}{c}\frac{\widehat{{\mathbf{X}}_i^{\prime }{\mathbf{M}}_{\mathbf{H}}\mathbf{\Omega }\mathbf{H}}}{T}=\frac{1}{TN}\sum _{j=1}^N[\sum _{t=1}^T{\widehat{e}}_{jt}^2{\widehat{\mathbf{x}}}_{it}{\widehat{\mathbf{h}}}_t^{\prime }+\sum _{s=1}^S\sum _{t=s+1}^T\left(1-\frac{s}{S+1}\right)\\ {\widehat{e}}_{jt}{\widehat{e}}_{j,t-s}\left({\widehat{\mathbf{x}}}_{it}{\widehat{\mathbf{h}}}_{t-s}^{\prime }+{\widehat{\mathbf{x}}}_{i,t-s}{\widehat{\mathbf{h}}}_t^{\prime }\right)]\end{array}$$(12) . In the Monte Carlo experiment in Section 4, we set $S=\lfloor T^{1/4}\rfloor$.

Remark 3.

Recently, jackknife bias-correction, a numerical way of asymptotically eliminating the asymptotic biases, was proposed and applied in Fernández and Weidner (2016) and Westerlund (2018), among others. The bias terms in Equation (11)(11) $$\tilde{\mathit{\beta }}=\widehat{\mathit{\beta }}-\frac{1}{N}{\widehat{\mathit{\xi }}}_{NT}-\frac{1}{T}{\widehat{\mathit{\zeta }}}_{NT},$$(11) are $O_p(N^{-1})$ and $O_p(T^{-1})$, which indicates that the jackknife bias-correction of Dhaene and Jochmans (2015) by splitting the time-series and cross-section dimensions would be applicable for $\widehat{\mathit{\beta }}$. We provide the procedure for this jackknife bias-correction and examine its finite sample performance in Section 4.

Now consider the robust approach against slope heterogeneity. Such robustness may be achieved when the leading terms in Equations (9)(9) $$\begin{array}{c}\widehat{\mathit{\beta }}-{\mathit{\beta }}^0={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathit{\epsilon }}_i\\ + \frac{1}{N}{\mathit{\xi }}_{NT}+\frac{1}{T}{\mathit{\zeta }}_{NT}+O_p({\delta }_{N\mathrm{}T}^{-3}),\end{array}$$(9) and (10)(10) $$\widehat{\mathit{\beta }}-{\mathit{\beta }}^0={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i\right)}^{-1}\frac{1}{NT}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i{\mathit{\eta }}_i+O_p({\delta }_{NT}^{-2}).$$(10) share a common expression, and the estimated asymptotic bias terms are bounded for the models with heterogeneous slopes. Indeed, we prove that for the models with heterogeneous slopes $\frac{1}{NT}{\sum }_{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i=\frac{1}{NT}{\sum }_{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i+o_p(1)$, $\frac{1}{\sqrt{N}T}{\sum }_{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i=\frac{1}{\sqrt{N}T}{\sum }_{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i{\mathit{\eta }}_i+o_p(1)$, ${\widehat{\mathit{\zeta }}}_{NT}=O_p(1)$ and ${\widehat{\mathit{\zeta }}}_{NT}=O_p(1)$, and for the models with homogeneous slopes $\frac{1}{\sqrt{NT}}{\sum }_{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i=\frac{1}{\sqrt{NT}}{\sum }_{i=1}^N$ ${\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathit{\epsilon }}_i$ under Assumptions 1–5; see the supplementary materials.

The following proposition shows a common stochastic representation of the bias-corrected IPC estimator for the models with homogeneous and heterogeneous slopes.

Proposition 2.

Suppose that Assumptions 1–5 hold. Then, we have(13) $$\begin{array}{c}\tilde{\mathit{\beta }}-{\mathit{\beta }}^0={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\\ \times \frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i+O_p({\delta }_{NT}^{-ϵ})\end{array}$$(13) where ϵ = 3 for the models with homogeneous slopes and ϵ = 2 for the models with heterogeneous slopes, as $(N,T)\to \infty$ and $T/N\to \rho \in (0,\Delta ]$.

Equipped with Assumption 6, we are ready to present the asymptotic normality of the bias-corrected estimator $\tilde{\mathit{\beta }}$.

Theorem 1.

Suppose that Assumptions 1–6 hold and $T/N\to \rho \in (0,\Delta ]$.

1. For the model with homogeneous slopes, we have$$\sqrt{NT}(\tilde{\mathit{\beta }}-{\mathit{\beta }}^0)\stackrel{d}{\to }N(\mathbf{0},{\mathbf{A}}_1^{-1}{\mathbf{B}}_1{\mathbf{A}}_1^{-1}).$$

2. For the model with heterogeneous slopes, we have$$\sqrt{N}(\tilde{\mathit{\beta }}-{\mathit{\beta }}^0)\stackrel{d}{\to }N(\mathbf{0},{\mathbf{A}}_0^{-1}{\mathbf{C}}_0{\mathbf{A}}_0^{-1}),$$

   where ${\mathbf{A}}_0,{\mathbf{A}}_1,{\mathbf{B}}_1$, and ${\mathbf{C}}_0$ are defined in Assumption 5.

Let us now turn our attention to the estimation of the asymptotic variance-covariance matrix of $\tilde{\mathit{\beta }}$. Define the panel heteroscedasticity and autocorrelation consistent (PHAC) variance-covariance matrix estimator by(14) $$\begin{array}{c}{\widehat{\mathbf{\Sigma }}}_{\tilde{\beta }}={\left(\sum _{i=1}^N{\widehat{\mathbf{Z}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\widehat{\mathbf{Z}}}_i\right)}^{-1}\left(\sum _{i=1}^N{\widehat{\mathbf{Z}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\widehat{\mathbf{u}}}_i{\widehat{\mathbf{u}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\widehat{\mathbf{Z}}}_i\right)\\ \times {\left(\sum _{i=1}^N{\widehat{\mathbf{Z}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\widehat{\mathbf{Z}}}_i\right)}^{-1},\end{array}$$(14) where ${\widehat{\mathbf{u}}}_i={\mathbf{y}}_i-{\mathbf{X}}_i\widehat{\mathit{\beta }}$. Note that this PHAC estimator does not require Newey and West (1987) type time-series parameter truncation.

One of the highlights of this article is the revelation of the asymptotic validity of the PHAC variance estimator for both homogeneous-slope and heterogeneous-slope models with interactive effects.⁵ This may be intuitively understood by the following discussion, with the help of Propositions 1 and 2 and Theorem 1. For the models with homogeneous slopes, Equations (9)(9) $$\begin{array}{c}\widehat{\mathit{\beta }}-{\mathit{\beta }}^0={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathit{\epsilon }}_i\\ + \frac{1}{N}{\mathit{\xi }}_{NT}+\frac{1}{T}{\mathit{\zeta }}_{NT}+O_p({\delta }_{N\mathrm{}T}^{-3}),\end{array}$$(9) and (13)(13) $$\begin{array}{c}\tilde{\mathit{\beta }}-{\mathit{\beta }}^0={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\\ \times \frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i+O_p({\delta }_{NT}^{-ϵ})\end{array}$$(13) yield$$\begin{array}{c}\sqrt{NT}(\tilde{\mathit{\beta }}-{\mathit{\beta }}^0)={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\\ \times \frac{1}{\sqrt{NT}}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i+o_p(1)\\ ={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\\ \times \frac{1}{\sqrt{NT}}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathit{\epsilon }}_i+o_p(1).\end{array}$$

From this, we can easily infer that $NT{\widehat{\mathbf{\Sigma }}}_{\tilde{\beta }}\stackrel{p}{\to }{\mathbf{A}}_1^{-1}{\mathbf{B}}_1{\mathbf{A}}_1^{-1}$. For the models with heterogeneous slopes, Equations (10)(10) $$\widehat{\mathit{\beta }}-{\mathit{\beta }}^0={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i\right)}^{-1}\frac{1}{NT}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i{\mathit{\eta }}_i+O_p({\delta }_{NT}^{-2}).$$(10) and (13)(13) $$\begin{array}{c}\tilde{\mathit{\beta }}-{\mathit{\beta }}^0={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\\ \times \frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i+O_p({\delta }_{NT}^{-ϵ})\end{array}$$(13) give(15) $$\begin{array}{c}\sqrt{N}(\tilde{\mathit{\beta }}-{\mathit{\beta }}^0)={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\\ \times \frac{1}{\sqrt{N}T}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i+o_p(1)\\ ={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i\right)}^{-1}\\ \times \frac{1}{\sqrt{N}T}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i{\mathit{\eta }}_i+o_p(1),\end{array}$$(15) from which, we can conjecture that $N{\widehat{\mathbf{\Sigma }}}_{\tilde{\beta }}\stackrel{p}{\to }{\mathbf{A}}_0^{-1}{\mathbf{C}}_0{\mathbf{A}}_0^{-1}$.

In the following theorem, we formally establish the limit distribution of the Wald test statistic for linear restrictions on ${\mathit{\beta }}^0$, based on the bias-corrected IPC estimator $\tilde{\mathit{\beta }}$ and the PHAC variance estimator ${\widehat{\mathbf{\Sigma }}}_{\tilde{\beta }}$.

Theorem 2.

Consider testing m linearly independent restrictions on ${\mathit{\beta }}^{\mathit{0}}$, $H_0:\mathbf{ℛ}{\mathit{\beta }}^{\mathit{0}}=\mathit{r}$ against $H_1:\mathbf{ℛ}{\mathit{\beta }}^{\mathit{0}}\ne \mathit{r}$, where $\mathbf{ℛ}$ is an m × k fixed matrix of full row rank. Consider the model (7) and the Wald statistic(16) $${\tilde{W}}_{NT}={(\mathbf{ℛ}\tilde{\mathit{\beta }}-\mathit{r})}^{\prime }{\left(\mathbf{ℛ}{\widehat{\mathbf{\Sigma }}}_{\tilde{\beta }}{\mathbf{ℛ}}^{\prime }\right)}^{-1}(\mathbf{ℛ}\tilde{\mathit{\beta }}-\mathit{r})$$(16) where $\tilde{\mathit{\beta }}$ and ${\widehat{\mathbf{\Sigma }}}_{\overline{\beta }}$ are defined by (11) and (14), respectively. Suppose that Assumptions 1–6 hold. Then, under the null hypothesis, for the models with homogeneous slopes (${\mathit{\eta }}_i=\mathbf{0}$ for all i) and for the models with heterogeneous slopes, ${\tilde{W}}_{NT}\stackrel{d}{\to }{\chi }_m^2$, as $(N,T)\to \infty$ and $T/N\to \rho \in (0,\Delta ]$.

To our knowledge, this is the first article which shows the consistency of the PHAC estimator in the random coefficients models with interactive effects for large panel data. As discussed in the introduction section, the use of the PHAC variance estimator for the models with homogeneous slopes with large N and T is proposed by Hansen (2007) for the fixed effects models, and considered by Greenaway-McGrevy, Han, and Sul (2012), Vogelsang (2012), Fernández and Weidner (2016), and Cui et al. (2021) for the models with homogeneous slopes, among others.

Remark 4.

Our approach is also robust against mixtures of homogeneous and heterogeneous slopes.⁶ To illustrate this, consider the model without common components and the case in which the k slopes are partitioned in such a way that $k=k_1+k_2$, where ${\mathit{\beta }}_i={({\mathit{\beta }}_{1i}^{\prime },{\mathit{\beta }}_2^{0\prime })}^{\prime },{\mathit{\beta }}_{1i}={\mathit{\beta }}_1^0+{\mathit{\eta }}_{1i},\mathbb{E}({\mathit{\eta }}_{1i})=\mathbf{0}$ and $\text{var}({\mathit{\eta }}_{1i})={\mathbf{\Omega }}_{1i}$, with ${\mathit{\beta }}^0\mathit{=}{({\mathit{\beta }}_1^{0\prime },{\mathit{\beta }}_2^{0\prime })}^{\prime }$. The expansion of the pooled ordinary least squares (OLS) estimator $\widehat{\mathit{\beta }}=({\widehat{\mathit{\beta }}}_1^{\prime },{\widehat{\mathit{\beta }}}_2^{\prime })\prime$ gives$$\begin{array}{c}\sqrt{N}\left(\widehat{\mathit{\beta }}-\mathit{\beta }\right)={\left(\frac{\sum _{i=1}^N{\mathbf{X}}_i^{\prime }{\mathbf{X}}_i}{NT}\right)}^{-1}\frac{1}{\sqrt{N}T}\sum _{i=1}^N{\mathbf{X}}_i^{\prime }{\mathbf{u}}_i\\ ={\left(\frac{\sum _{i=1}^N{\mathbf{X}}_i^{\prime }{\mathbf{X}}_i}{NT}\right)}^{-1}\frac{1}{\sqrt{N}}\sum _{i=1}^N\left(\frac{{\mathbf{X}}_i^{\prime }{\mathbf{X}}_{1i}}{T}\right){\mathit{\eta }}_{1i}\\ + O_p\left(1/\sqrt{T}\right)\\ \stackrel{d}{\to }N\left(\mathbf{0},{\mathbf{A}}^{-1}\mathbf{C}{\mathbf{A}}^{-1}\right)\end{array}$$where $\mathbf{A}=\underset{N,T\to \infty }{\text{plim}}{(NT)}^{-1}{\sum }_{i=1}^N{\mathbf{X}}_i^{\prime }{\mathbf{X}}_i$, which is accordingly partitioned as$$\begin{array}{c}\mathbf{A}=\left[\begin{array}{cc}{\mathbf{A}}_{11}& {\mathbf{A}}_{12}\\ {\mathbf{A}}_{21}& {\mathbf{A}}_{22}\end{array}\right],\hspace{1em}\text{and}\\ \mathbf{C}=\underset{N,T\to \infty }{\text{plim}}\frac{1}{N}\sum _{i=1}^N\left(\frac{{\mathbf{X}}_i^{\prime }{\mathbf{X}}_{1i}}{T}\right){\mathit{\eta }}_{1i}{\mathit{\eta }}_{1i}^{\prime }\left(\frac{{\mathbf{X}}_{1i}^{\prime }{\mathbf{X}}_i}{T}\right),\end{array}$$ which are assumed to be fixed and positive definite. Also note that the convergence rate of (${\widehat{\mathit{\beta }}}_1$ and) ${\widehat{\mathit{\beta }}}_2$ is $\sqrt{N}$, as the variation of $\widehat{\mathit{\beta }}$ is dominated by that of ${\mathit{\eta }}_{1i}$.

Now consider a special case in which ${\text{plim}}_{T\to \infty }{\mathbf{X}}_{1i}^{\prime }{\mathbf{X}}_{2i}=\mathbf{0}$. Then, a similar discussion gives $\sqrt{N}({\widehat{\mathit{\beta }}}_1-{\mathit{\beta }}_1^0)\mathrm{}\stackrel{d}{\to }N$ $(\mathbf{0},{\mathbf{A}}_{11}^{-1}{\mathbf{C}}_{11}{\mathbf{A}}_{11}^{-1})$ and $\sqrt{NT}(\mathrm{}{\widehat{\mathit{\beta }}}_2-{\mathit{\beta }}_2^0)\stackrel{d}{\to }N(\mathbf{0},{\mathbf{A}}_{22}^{-1}{\mathbf{B}}_{22}{\mathbf{A}}_{22}^{-1})$, where ${\mathbf{C}}_{11}=\text{plim}N,T\to \infty {(N{T}^2)}^{-1}{\sum }_{i=1}^N{\mathbf{X}}_{1i}^{\prime }{\mathbf{X}}_{1i}{\mathit{\eta }}_{1i}{\mathit{\eta }}_{1i}^{\prime }$ ${\mathbf{X}}_{1i}^{\prime }{\mathbf{X}}_{1i}$ and ${\mathbf{B}}_{22}={\text{plim}}_{N,T\to \infty }{(NT)}^{-1}{\sum }_{i=1}^N{\mathbf{X}}_{2i}^{\prime }{\mathit{\epsilon }}_i{\mathit{\epsilon }}_i^{\prime }{\mathbf{X}}_{2i}$. Therefore, the asymptotic normality of $\widehat{\mathit{\beta }}$, the consistency of the PHAC estimator and the asymptotic validity of Wald test hold with mixtures of homogeneous and heterogeneous slopes.

3 LM Test for Correlated Random Coefficients with Covariates

As discussed earlier, the proposed robust approach is asymptotically justified for the models with homogeneous slopes and for random coefficient models, in which the slope heterogeneity is generally independent of the variations of the regressors (see Assumption 3). In particular, when the heterogeneous coefficients are correlated with the defactored regressor matrix, ${\mathbf{M}}_{{\mathbf{H}}^{\mathbf{0}}}{\mathbf{X}}_i$, the approach may not work. To see this, consider the bias-corrected IPC estimator $\tilde{\mathit{\beta }}$ defined by Equation (11)(11) $$\tilde{\mathit{\beta }}=\widehat{\mathit{\beta }}-\frac{1}{N}{\widehat{\mathit{\xi }}}_{NT}-\frac{1}{T}{\widehat{\mathit{\zeta }}}_{NT},$$(11) for the models with heterogeneous slopes. Equation (15)(15) $$\begin{array}{c}\sqrt{N}(\tilde{\mathit{\beta }}-{\mathit{\beta }}^0)={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\\ \times \frac{1}{\sqrt{N}T}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i+o_p(1)\\ ={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i\right)}^{-1}\\ \times \frac{1}{\sqrt{N}T}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i{\mathit{\eta }}_i+o_p(1),\end{array}$$(15) gives $\tilde{\mathit{\beta }}-{\mathit{\beta }}^0\approx {({\sum }_{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i)}^{-1}{\sum }_{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i{\mathit{\eta }}_i$. If $\mathbb{E}({\mathit{\eta }}_i|{\mathbf{V}}_i)\ne \mathbf{0}$, which violates the random coefficient assumption (in Assumption 3), $\mathbb{E}\left[{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i\mathbb{E}({\mathit{\eta }}_i|{\mathbf{V}}_i)\right]$ is not necessarily a zero vector, and this may render $\tilde{\mathit{\beta }}$ inconsistent, in general. Note that even if we had nonzero correlations between ${\mathit{\eta }}_i$ and the factor loadings, the slope estimator would remain consistent so long as the regressor matrix is defactored, because the leading term in Equation (15)(15) $$\begin{array}{c}\sqrt{N}(\tilde{\mathit{\beta }}-{\mathit{\beta }}^0)={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i\right)}^{-1}\\ \times \frac{1}{\sqrt{N}T}\sum _{i=1}^N{\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i+o_p(1)\\ ={\left(\frac{1}{NT}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i\right)}^{-1}\\ \times \frac{1}{\sqrt{N}T}\sum _{i=1}^N{\mathbf{V}}_i^{\prime }{\mathbf{V}}_i{\mathit{\eta }}_i+o_p(1),\end{array}$$(15) does not contain the factor loadings. Therefore, we focus on the correlated random coefficients caused by $\mathbb{E}({\mathit{\eta }}_i|{\mathbf{V}}_i)\ne \mathbf{0}$.

There are studies which propose tests to detect slope heterogeneity for large panel data. However, most of them have focused only on testing the slope heterogeneity, not the dependence of the slope heterogeneity on the regressors. Such examples are found in Pesaran and Yamagata (2008) and Su and Chen (2013), among others. These existing tests are not suitable for our robust approach, because they have power against the random coefficient models as well as the correlated random coefficient models.

In view of this, we propose a novel Lagrange Multiplier (LM) or Score test of dependence of random coefficients on covariates. The proposed test is designed to detect departures from the assumption of independence between slope heterogeneity and the regressors, or departures from the assumption of slope homogeneity toward slope heterogeneity dependent on the regressors. Hereafter, for simplicity, we call it the LM test for correlated random coefficients (CRC). Specifically, we extend the test of correlated random effects proposed in (Wooldridge 2010, chap. 11.7.4) for short panels in the fixed effects models. Therein, using a similar approach taken by Mundlak (1978) and Chamberlain (1980), the correlated heterogeneity is modeled as a linear function of within average of the regressors, ${\mathit{\beta }}_i=\mathit{\mu }+\theta {\overline{\mathbf{x}}}_i+{\mathit{\omega }}_i$, where ${\overline{\mathbf{x}}}_i=T^{-1}{\sum }_{t=1}^T{\mathbf{x}}_{it}$ and ${\mathit{\omega }}_i$ is a mean-zero random vector which is uncorrelated with ${\overline{\mathbf{x}}}_i$. From equations (4)(4) $${\mathbf{X}}_i={\mathbf{G}}^0{\mathbf{\Gamma }}_i^0+{\mathbf{V}}_i,$$(4) and (5)(5) $${\mathit{\beta }}_i={\mathit{\beta }}^0+{\mathit{\eta }}_i,$$(5) , it is easily seen that the error term ${\mathbf{u}}_i$ for the model of omitted slope heterogeneity will contain the term ${\mathbf{X}}_i{\mathit{\beta }}_i$ and this motivates his test of correlated random coefficients, which examines the significance of ${\mathbf{X}}_i{\overline{\mathbf{x}}}_i$. As pointed out by (Wooldridge 2010, p. 386), this test does not have power against uncorrelated slope heterogeneity, which is the desirable property for our approach. We extend the test by considering the models with interactive effects for large N and T and by permitting the alternatives, $\mathbb{E}({\mathit{\eta }}_i|{\mathbf{V}}_i)$, to be a nonlinear function of ${\mathbf{V}}_i$.

To motivate the LM test for correlated random coefficients, suppose that the random part of the coefficients, ${\mathit{\eta }}_i$ (see Equation (5)(5) $${\mathit{\beta }}_i={\mathit{\beta }}^0+{\mathit{\eta }}_i,$$(5) ), is decomposed into correlated and uncorrelated parts:(17) $${\mathit{\eta }}_i=\mathit{\psi }\left({\mathbf{V}}_i\right)+{\mathit{\omega }}_i,$$(17) where $\mathbb{E}[\mathit{\psi }({\mathbf{V}}_i)|{\mathbf{V}}_i]\ne \mathbf{0}$ but $\mathbb{E}[{\mathit{\omega }}_i|{\mathbf{V}}_i]=\mathbf{0}$. Various functional forms of $\mathit{\psi }(.)$ can be entertained. As an illustration, suppose the $\mathcal{l}$th element of $\mathit{\psi }({\mathbf{V}}_i)$ to be ${\psi }_{\mathcal{l}i}=\frac{1}{\sqrt{T}}{\sum }_{t=1}^T\{f(v_{\mathcal{l}it})-\mathbb{E}[f(v_{\mathcal{l}it})]\}$, where the function $f(.)$ satisfies certain regularity conditions. Then, the qth-order Maclaurin series expansion of $f(v_{\mathcal{l}it})-\mathbb{E}[f(v_{\mathcal{l}it})]$ gives $\mathrm{(}\text{Tex}\mathrm{ }\text{translation}\mathrm{ }\text{failed}\mathrm{)}$, where ${\theta }_p=\frac{1}{p!}{f}^{(p)}(0),f^{(p)}(.)$ is the pth-order derivative of $f(.)$, and R_q is the reminder term. For testing purposes, examining whether the first few terms of the expansion, say $g(\le q)$ terms, are zero or not may suffice. This motivates our test of correlated random coefficients (defined below) for the null hypothesis,(18) $$H_0:\mathit{\theta }=\mathbf{0}\mathrm{ }\text{against}\mathrm{ }{H}_1:\mathit{\theta }\ne \mathbf{0},$$(18) where $\mathit{\theta }=({\theta }_1,\dots ,{\theta }_g)\prime$.⁷

Having discussed a process of the correlated random coefficients, let us turn our attention to the test statistic. Following on from the above discussion, we define the infeasible score of the LM test for the models with homogeneous (respectively, heterogeneous) slopes as $\frac{1}{\sqrt{NT}}{\sum }_{i=1}^N{\mathbf{L}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i$ (resp. $\frac{1}{\sqrt{N}T}{\sum }_{i=1}^N{\mathbf{L}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{u}}_i$), where ${\mathbf{L}}_i={\mathbf{X}}_i\sqrt{T}({\mathbf{\Xi }}_i-\mathbf{\Xi })$,(19) $${\mathbf{\Xi }}_i=\left({\overline{\mathbf{v}}}_i^{\left(1\right)},{\overline{\mathbf{v}}}_i^{\left(2\right)},\dots ,{\overline{\mathbf{v}}}_i^{\left(g\right)}\right),$$(19)

$\mathbf{\Xi }=\underset{N\to \infty }{\lim }{N}^{-1}{\sum }_{i=1}^N\mathbb{E}[{\mathbf{\Xi }}_i]$, with ${\overline{\mathbf{v}}}_i^{(p)}=({\overline{v}}_{1i}^{(p)},{\overline{v}}_{2i}^{(p)},\dots ,$ ${{\overline{v}}_{ki}^{(p)})}^{\prime },{\overline{v}}_{\mathcal{l}i}^{(p)}=T^{-1}{\sum }_{t=1}^T{v}_{\mathcal{l}it}^p,v_{\mathcal{l}it}^p$ is the pth power of the $(t,\mathcal{l})$ element of ${\mathbf{V}}_i,p=1,\dots ,g$. Observe that the direct alternative model of $\mathit{\psi }({\mathbf{V}}_i)$ for this score is $\sqrt{T}({\mathbf{\Xi }}_i-\mathbf{\Xi })\mathit{\theta }$, where $\mathit{\theta }$ is the $g\times 1$ constant vector.

To construct the associated feasible score, ${\mathbf{L}}_i$, ${\mathbf{H}}^0$, and ${\mathbf{u}}_i$ should be replaced by their sample counterparts, which are ${\widehat{\mathbf{L}}}_i={\mathbf{X}}_i\sqrt{T}({\mathbf{\Xi }}_i-\overline{\mathbf{\Xi }})$, $\widehat{\mathbf{H}}$, and $\tilde{\mathbf{u}}={\mathbf{y}}_i-{\mathbf{X}}_i\tilde{\mathit{\beta }}$, respectively, where ${\widehat{\mathbf{\Xi }}}_i=({\widehat{\overline{\mathbf{v}}}}_i^{(1)},{\widehat{\overline{\mathbf{v}}}}_i^{(2)},\dots ,{\widehat{\overline{\mathbf{v}}}}_i^{(g)}),{\widehat{\overline{\mathbf{v}}}}_i^{(p)}={({\widehat{\overline{v}}}_{1i}^{(p)},{\widehat{\overline{v}}}_{2i}^{(p)},\dots ,{\widehat{\overline{v}}}_{ki}^{(p)})}^{\prime },{\widehat{\overline{v}}}_{\mathcal{l}i}^{(p)}=T^{-1}{\sum }_{t=1}^T{\widehat{v}}_{\mathcal{l}it}^p$ with ${\widehat{\mathbf{V}}}_i={\mathbf{M}}_{\widehat{\mathbf{H}}}{\mathbf{X}}_i=({\widehat{v}}_{\mathcal{l}it}),p=1,\dots ,g$. The resulting feasible score for the models with homogeneous (resp. heterogeneous) slopes is $\frac{1}{\sqrt{NT}}{\sum }_{i=1}^N{\widehat{\mathbf{L}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\tilde{\mathbf{u}}}_i$ (respectively, $\frac{1}{\sqrt{N}T}{\sum }_{i=1}^N{\widehat{\mathbf{L}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\tilde{\mathbf{u}}}_i$). Let us sketch the derivation of the LM test statistic next. For the models with homogeneous slopes, a familiar expansion gives ${\tilde{\mathbf{u}}}_i={\mathbf{u}}_i-{\mathbf{X}}_i(\tilde{\mathit{\beta }}-{\mathit{\beta }}^0)$, which yields$$\begin{array}{c}\frac{1}{\sqrt{NT}}\sum _{i=1}^N{\widehat{\mathbf{L}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\tilde{\mathbf{u}}}_i=\frac{1}{\sqrt{NT}}\sum _{i=1}^N{\widehat{\mathbf{L}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\mathbf{u}}_i\\ +\frac{1}{NT}\sum _{i=1}^N{\widehat{\mathbf{L}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\mathbf{X}}_i\sqrt{NT}\left(\tilde{\mathit{\beta }}-{\mathit{\beta }}^0\right).\end{array}$$

For a reason similar to why $\widehat{\mathit{\beta }}$ has asymptotic bias terms, the first term of the right hand side of the above equation will have asymptotic bias terms. For the second term, demeaning ${\mathbf{L}}_i$ and ${\mathbf{X}}_i$ is required for a reason similar to why the asymptotic variance of $\widehat{\mathit{\beta }}$ involves ${\mathbf{Z}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{Z}}_i$ rather than ${\mathbf{X}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{X}}_i$. Define such a demeaned ${\mathbf{L}}_i$, ${\mathfrak{L}}_i={\mathbf{L}}_i-N^{-1}{\sum }_{j=1}^N{\mathit{\varphi }}_i^{0\prime }{\mathbf{\Upsilon }}_{{\varphi }^0}^{-1}{\mathit{\varphi }}_j^0{\mathbf{L}}_j$. For the LM test, we impose the following assumption, which is similar to Assumption 6:

Assumption 7.

When $N,T\to \infty$such that $T/N\to \rho \in (0,\Delta ]$, $\frac{1}{\sqrt{NT}}{\sum }_{i=1}^N{\mathfrak{L}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathit{\epsilon }}_i\stackrel{d}{\to }N(\mathbf{0},{\mathbf{B}}_2)$, where ${\mathbf{B}}_2={\text{plim}}_{N,T\to \infty }\frac{1}{NT}{\sum }_{i=1}^N{\mathfrak{L}}_i^{\prime }{\mathbf{M}}_{{\mathbf{H}}^0}{\mathbf{\Omega }}_{\epsilon ,i}{\mathbf{M}}_{{\mathbf{H}}^0}{\mathfrak{L}}_i$, which is positive definite.

For the slope heterogeneous case, the use of the same score scaled by $\frac{1}{\sqrt{T}}$ will be justified. For more details, see proof of Theorem 3 in the supplementary materials.

Consequently, the LM test statistic for correlated random coefficients, which is common to the model with homogeneous slopes and the random coefficient model, is defined by(20) $${\text{LM}}_{\text{CRC}}^{\left(g\right)}={\widehat{\mathbf{s}}}^{†\prime }{\widehat{\mathbf{\Sigma }}}_{ss}^{-1}{\widehat{\mathbf{s}}}^{†}$$(20) where ${\widehat{\mathbf{s}}}^{†}=\widehat{\mathbf{s}}-N^{-1}{\widehat{\mathit{\xi }}}_{NT}^{†}-T^{-1}{\widehat{\mathit{\zeta }}}_{NT}^{†}$ with $\widehat{\mathbf{s}}={\sum }_{i=1}^N{\widehat{\mathbf{L}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\tilde{\mathbf{u}}}_i$,$$\begin{array}{c}{\mathit{\xi }}_{NT}^{†}=-\frac{1}{N{T}^2}\sum _{i=1}^N\sum _{t=1}^T{\widehat{\mathfrak{L}}}_i^{\prime }\widehat{\mathbf{H}}{\widehat{\mathbf{\Upsilon }}}_{\varphi }^{-1}{\widehat{\mathit{\varphi }}}_i{\widehat{e}}_{it}^2,\mathrm{ }{\mathit{\zeta }}_{NT}^{†}\\ =-\frac{1}{N}\sum _{i=1}^N\frac{\widehat{{\mathfrak{L}}_i^{\prime }{\mathbf{M}}_{\mathbf{H}}\mathbf{\Omega }\mathbf{H}}}{T}{\widehat{\mathbf{\Upsilon }}}_{\varphi }^{-1}{\widehat{\mathit{\varphi }}}_i,\end{array}$$

${\widehat{\mathfrak{L}}}_i={\widehat{\mathbf{L}}}_i-N^{-1}{\sum }_{j=1}^N{\widehat{\mathit{\varphi }}}_i^{\prime }{\widehat{\mathbf{\Upsilon }}}_{\varphi }^{-1}{\widehat{\mathit{\varphi }}}_j{\widehat{\mathbf{L}}}_j$, $\frac{\widehat{{\mathfrak{L}}_i^{\prime }{\mathbf{M}}_{\mathbf{H}}\mathbf{\Omega }\mathbf{H}}}{T}$ is computed analogously to the computation of $\frac{\widehat{{\mathbf{X}}_i^{\prime }{\mathbf{M}}_{\mathbf{H}}\mathbf{\Omega }\mathbf{H}}}{T}$ (see Equation (12)(12) $$\begin{array}{c}\frac{\widehat{{\mathbf{X}}_i^{\prime }{\mathbf{M}}_{\mathbf{H}}\mathbf{\Omega }\mathbf{H}}}{T}=\frac{1}{TN}\sum _{j=1}^N[\sum _{t=1}^T{\widehat{e}}_{jt}^2{\widehat{\mathbf{x}}}_{it}{\widehat{\mathbf{h}}}_t^{\prime }+\sum _{s=1}^S\sum _{t=s+1}^T\left(1-\frac{s}{S+1}\right)\\ {\widehat{e}}_{jt}{\widehat{e}}_{j,t-s}\left({\widehat{\mathbf{x}}}_{it}{\widehat{\mathbf{h}}}_{t-s}^{\prime }+{\widehat{\mathbf{x}}}_{i,t-s}{\widehat{\mathbf{h}}}_t^{\prime }\right)]\end{array}$$(12) ), and ${\widehat{\mathbf{\Sigma }}}_{ss}={\sum }_{i=1}^N{\widehat{\mathbf{K}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\tilde{\mathbf{u}}}_i{\tilde{\mathbf{u}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\widehat{\mathbf{K}}}_i$ with$${\widehat{\mathbf{K}}}_i={\widehat{\mathfrak{L}}}_i-{\widehat{\mathbf{Z}}}_i{\left(\sum _{i=1}^N{\widehat{\mathbf{Z}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\widehat{\mathbf{Z}}}_i\right)}^{-1}\sum _{i=1}^N{\widehat{\mathbf{Z}}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{H}}}{\widehat{\mathfrak{L}}}_i.$$

Now we are ready to present the results of the limiting distribution of the LM test statistic.

Theorem 3.

Consider the model (7). Suppose that Assumptions 1–7 hold. Then, under the null hypothesis specified in (18), ${\text{LM}}_{\text{CRC}}^{(g)}\stackrel{d}{\to }{\chi }_g^2$, as $N,T\to \infty$ and $T/N\to \rho \in (0,\Delta ]$, for both homogeneous-slope and heterogeneous-slope models.

In practice, the degree of polynomial, g, for the statistic ${\text{LM}}_{\text{CRC}}^{(g)}$, should be chosen. In order to approximate the function $\mathit{\psi }({\mathbf{V}}_i)$ to the extent that the LM test has sufficient power, a small value of g, such as two, may be sufficient. The experimental results reported below support the choice of g = 2. In addition, the asymptotic local power analysis, which is provided in the supplementary materials, suggests that the LM test has nontrivial power when $\left|\right|\mathit{\theta }\left|\right|\asymp \frac{1}{\sqrt{NT}}$ and that the LM test associated with the pth column of the test matrix ${\mathbf{L}}_i$ has nontrivial power only against the pth term in the expansion of $\mathit{\psi }({\mathbf{V}}_i)$, that is, $H_1:{\theta }_p\ne 0$. The experimental results shown below support these analytical results.

When the null hypothesis of random coefficient models is rejected in favor of the alternatives, it is preferable to employ estimators which are consistent to ${\mathit{\beta }}^0$. For the estimation of the models with interactive effects for large panel data, the CCE and 2SIV mean group estimators and the ML estimators proposed by Pesaran (2006), Cui et al. (2021), and Li, Cui, and Lu (2020), among others, would be possible choices.

4 Monte Carlo Experiments

In this section we investigate the finite sample performance of our robust approach against slope heterogeneity, error serial correlation and heteroscedasticity. We consider the performance of the (analytically) bias-corrected IPC estimator defined by (11), which is denoted as BC-IPC. In addition, as discussed in Remark 3, we introduce a jackknife bias-corrected IPC estimator, which is denoted as JK-IPC. To evaluate accuracy of the bias-correction, the IPC estimator without bias-corrections is considered as well. In particular, we examine biases, standard deviations and root mean square errors (RMSEs) of these estimators, and empirical size of the (Wald) test for linear restrictions of slope coefficients. Furthermore, performance of the LM test for correlated random coefficients is examined.

Before introducing our experimental design, we discuss how to choose the initial estimator and the number of factors r for the IPC estimation in Section 4.1, then the procedure for computing the jackknife bias-corrected IPC estimator in Section 4.2.

4.1 Initial Value for the IPC Estimator and the Number of Factors

In practice, to compute the IPC estimator, an initial (consistent) estimator of ${\mathit{\beta }}^0$ and the number of factors in ${\mathbf{u}}_i$ must be given. Exploiting the factor structure in ${\mathbf{X}}_i$ and ${\mathbf{u}}_i$, we propose to employ the consistent PC estimator of Westerlund and Urbain (2015) and Reese and Westerlund (2018) as the initial estimator of ${\mathit{\beta }}^0$, which is defined as ${\widehat{\mathit{\beta }}}_{\text{PC}}={({\sum }_{i=1}^N{\mathbf{X}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{W}}}{\mathbf{X}}_i)}^{-1}{\sum }_{i=1}^N{\mathbf{X}}_i^{\prime }{\mathbf{M}}_{\widehat{\mathbf{W}}}{\mathbf{y}}_i$, where $\widehat{\mathbf{W}}$ is a $T\times \widehat{r_1+r_2}$ PC estimator of a span of ${\mathbf{W}}^0=({\mathbf{F}}^0,{\mathbf{G}}^0)$. $\widehat{r_1+r_2}$ and $\widehat{\mathbf{W}}$ are estimated using ${\{{\mathbf{X}}_i,{\mathbf{y}}_i\}}_{i=1}^N$, applying the Eigenvalue Ratio (ER) or the Growth Ratio (GR) estimators proposed by Ahn and Horenstein (2013), and the PC method, respectively. The estimate of the number of factors in ${\mathbf{u}}_i,\widehat{r}$, is obtained using the residual, ${\widehat{\mathbf{u}}}_{\text{PC}i}={\mathbf{y}}_i-{\mathbf{X}}_i{\widehat{\mathit{\beta }}}_{\text{PC}}$, applying the ER or GR methods.⁸ In the experiment below, we set the maximum number of factors to six to estimate $r_1+r_2$ and r. Alternatively, the CCE or 2SIV estimator of Pesaran (2006) and Cui et al. (2021) can be used as an initial estimator of ${\mathit{\beta }}^0$, which is not pursued in this article.

4.2 Jackknife Bias Correction

In Sections 2 and 3, the analytical bias correction of the IPC estimator $\widehat{\mathit{\beta }}$ is considered. Instead, as discussed in Remark 3, subsampling methods can be employed for bias corrections. Specifically, we consider the split panel jackknife in both the cross-section and time dimensions, which is proposed in Fernández and Weidner (2016).

To define the jackknife bias-corrected IPC estimator, assuming that N and T are even numbers for simplicity, consider the index set, $\mathcal{N}=\left\{1,\dots ,N\right\}$, which is divided into two subgroups, ${\mathcal{N}}_1=\left\{1,\dots ,N/2\right\}$ and ${\mathcal{N}}_2=\left\{(N/2)+1,\dots ,N\right\}$. Similarly define $\mathcal{T}=\left\{1,\dots ,T\right\},{\mathcal{T}}_1=\left\{1,\dots ,T/2\right\}$ and ${\mathcal{T}}_2=\left\{(T/2)+1,\dots ,T\right\}$. Let us denote the IPC estimator as ${\widehat{\mathit{\beta }}}_{\mathcal{N},\mathcal{T}}(=\widehat{\mathit{\beta }})$, in order to explicitly express the dependence on the (sub)sample, $(i,t)\in (\mathcal{N},\mathcal{T})$. Then, we define the jackknife bias-corrected IPC estimator as ${\widehat{\mathit{\beta }}}_{JK}=({\widehat{\mathit{\beta }}}_{TS}+{\widehat{\mathit{\beta }}}_{CS})-{\widehat{\mathit{\beta }}}_{\mathcal{N},\mathcal{T}}$, where ${\widehat{\mathit{\beta }}}_{TS}=2{\widehat{\mathit{\beta }}}_{\mathcal{N},\mathcal{T}}-\frac{{\widehat{\mathit{\beta }}}_{\mathcal{N},{\mathcal{T}}_1}+{\widehat{\mathit{\beta }}}_{\mathcal{N},{\mathcal{T}}_2}}{2}$ and ${\widehat{\mathit{\beta }}}_{CS}=2{\widehat{\mathit{\beta }}}_{\mathcal{N},\mathcal{T}}-\frac{{\widehat{\mathit{\beta }}}_{{\mathcal{N}}_1,\mathcal{T}}+{\widehat{\mathit{\beta }}}_{{\mathcal{N}}_2,\mathcal{T}}}{2}$. For more details, see the supplementary materials and Section 3.2 in Fernández and Weidner (2016).

4.3 Experimental Design

Consider the following data generating process:(21) $$\begin{array}{c}{y}_{it}=\sum _{\mathcal{l}=1}^k{x}_{\mathcal{l}it}{\beta }_{\mathcal{l}i}^0+f_{1t}^0{\lambda }_{1i}^0+f_{2t}^0{\lambda }_{2i}^0+{\sigma }_{\epsilon ,it}{\epsilon }_{it},\\ i=1,2,\dots ,N;t=1,2,\dots ,T\end{array}$$(21) where ${\lambda }_{si}^0\sim \text{iid}N(0,1),f_{st}^0={\rho }_f{f}_{s,t-1}^0+\sqrt{1-{\rho }_f^2}{ϵ}_{f,st},{ϵ}_{f,st}\sim \text{iid}N(0,1)$ with $f_{s0}^0\sim \text{iid}N(0,1)$ for $s=1,2,3,{\epsilon }_{it}={\rho }_{\epsilon }{\epsilon }_{it-1}+\sqrt{1-{\rho }_{\epsilon }^2}{ϵ}_{\epsilon ,it},{ϵ}_{\epsilon ,it}\sim \text{iid}N(0,1)$ with ${\epsilon }_{i0}\sim \text{iid}N(0,1)$, and(22) $${\sigma }_{\epsilon ,it}={\left({\kappa }_{\epsilon ,i}{\kappa }_{\epsilon ,t}\right)}^{1/2},\mathrm{ }{\kappa }_{\epsilon ,i}\sim \text{iid}U\left(0.5,1.5\right)\mathrm{ }\text{and}\mathrm{ }{\kappa }_{\epsilon ,t}=0.5+t/T.$$(22)

The regressors $x_{\mathcal{l}it},\mathcal{l}=1,2,\dots ,k$, are generated as(23) $$x_{\mathcal{l}it}=f_{1t}^0{\gamma }_{\mathcal{l}1i}^0+f_{3t}^0{\gamma }_{\mathcal{l}3i}^0+0.3{\sigma }_{v,it}{v}_{\mathcal{l}it},$$(23) where $v_{\mathcal{l}it}={\rho }_v{v}_{\mathcal{l}it-1}+\sqrt{1-{\rho }_v^2}{\varpi }_{\mathcal{l}it}$. Note that $f_{1t}^0$ is common to both y_it and $x_{\mathcal{l}it},$ whereas $f_{2t}^0$ and $f_{3t}^0$ are unique to y_it and $x_{\mathcal{l}it},$ respectively. ${\varpi }_{\mathcal{l}it}$ is generated as ${\varpi }_{\mathcal{l}it}=({\varpi }_{\mathcal{l}it}^{\ast }-c)/\sqrt{2c},{\varpi }_{\mathcal{l}it}^{\ast }\sim \text{iid}{\chi }_c^2$ and $v_{\mathcal{l}i0}=(v_{\mathcal{l}i0}^{\ast }-c)/\sqrt{2c},v_{\mathcal{l}i0}^{\ast }\sim \text{iid}{\chi }_c^2$ with $c=6.$⁹ The factor loadings in $x_{\mathcal{l}it}$ are generated as(24) $${\gamma }_{\mathcal{l}si}^0=0.7{\lambda }_{si}^0+{\left(1-{0.7}^2\right)}^{1/2}{\phi }_{\mathcal{l}si},$$(24)

${\phi }_{\mathcal{l}si}\sim \text{iid}N(0,1)$ for $\mathcal{l}=1,\dots ,k$ and $s=1,3,$ so that they are correlated with factor loadings in y_it .(25) $${\sigma }_{v,it}={\left({\kappa }_{v,i}{\kappa }_{v,t}\right)}^{1/2},{\kappa }_{v,i}\sim \text{iid}U\left(0.5,1.5\right)\mathrm{ }\text{and}\mathrm{ }{\kappa }_{v,t}=4.5+t/T.$$(25)

Finally we have(26) $${\beta }_{\mathcal{l}i}={\beta }_{\mathcal{l}}^0+{\sigma }_{\eta }{\eta }_{\mathcal{l}i},$$(26) (27) $${\eta }_{\mathcal{l}i}={\rho }_{x\eta }{\psi }_{\mathcal{l}i}+\sqrt{1-{\rho }_{x\eta }^2}{\omega }_{\mathcal{l}i}$$(27)

${\omega }_{\mathcal{l}i}\sim \text{iid}N(0,1)$ for $\mathcal{l}=1,\dots ,k$ and ${\psi }_{\mathcal{l}i}^{(\mathcal{P})}={\sum }_{p=1}^q{\theta }_p{\tilde{z}}_{\mathcal{l}i,p}/$ $\sqrt{{\sum }_{p=1}^q{\theta }_p^2}$, where $\mathcal{P}=\{p:{\theta }_p\ne 0,p\in I\}$ with $I=\{1,2,\dots ,q\}$ and ${\theta }_p\in \{0,1\}$,$${\tilde{z}}_{\mathcal{l}i,p}=\frac{z_{\mathcal{l}i,p}-{\overline{z}}_{\mathcal{l},p}}{s_{z\mathcal{l},p}}$$with ${\overline{z}}_{\mathcal{l},p}=N^{-1}{\sum }_{i=1}^N{z}_{\mathcal{l}i,p},s_{z\mathcal{l},p}^2={(N-1)}^{-1}{\sum }_{i=1}^N$ ${(z_{\mathcal{l}i,p}-{\overline{z}}_{\mathcal{l},p})}^2$ and $z_{\mathcal{l}i,p}=T^{-1}{\sum }_{t=1}^T{v}_{\mathcal{l}it}^p$. Observe that the variance of ${\beta }_{\mathcal{l}i}$ is controlled to be ${\sigma }_{\eta }^2$. θ_p selects the term ${\tilde{z}}_{\mathcal{l}i,p}$. Specifically, in the experiment, we consider ${\psi }_{\mathcal{l}i}^{(1,2,3,4)}={\sum }_{p=1}^4{\tilde{z}}_{\mathcal{l}i,p}/\sqrt{4}$ and ${\psi }_{\mathcal{l}i}^{(p)}={\tilde{z}}_{\mathcal{l}i,p},p=1,\dots ,4$. It is interesting to examine the impact of such linear and highly nonlinear relationships on the estimators and the power of the LM tests.

We set k = 2 (two regressors) for all experiments. Before the estimation, the data is within-transformed and cross-sectionally demeaned, to make the results invariant to the inclusion of (additive) individual effects and time effects. For parameter values, we set $({\beta }_1^0,{\beta }_2^0)=(1,1)$.

We consider seven experimental designs as outlined in Table 1. Designs 1–7 are considered to investigate the finite sample behavior of the slope estimator, the associated robust Wald test and the LM test for correlated random coefficients.

Table 1 Experimental designs.

Design	Model	r	${\sigma }_{\eta }$	${\rho }_{x\eta }$	${\psi }_{\mathcal{l}i}^{(\mathcal{P})}$
1	Homogeneous slopes	2	0	0	–
2	Heterogeneous slopes uncorrelated with $x_{\mathcal{l}it}$	3	0.50	0	–
3	Heterogeneous slopes correlated with $x_{\mathcal{l}it}$	3	0.50	0.50	${\psi }_{\mathcal{l}i}^{(1,2,3,4)}$
4	Heterogeneous slopes correlated with $x_{\mathcal{l}it}$	3	0.50	0.50	${\psi }_{\mathcal{l}i}^{(1)}$
5	Heterogeneous slopes correlated with $x_{\mathcal{l}it}$	3	0.50	0.50	${\psi }_{\mathcal{l}i}^{(2)}$
6	Heterogeneous slopes correlated with $x_{\ell it}$	3	0.50	0.50	${\psi }_{\mathcal{l}i}^{(3)}$
7	Heterogeneous slopes correlated with $x_{\mathcal{l}it}$	3	0.50	0.50	${\psi }_{\mathcal{l}i}^{(4)}$

NOTE: r is the number of factors in the error term. ${\sigma }_{\eta }$ is standard deviation of ${\beta }_{\mathcal{l}i}-{\beta }_{\mathcal{l}}^0(={\sigma }_{\eta }{\eta }_{\mathcal{l}i})$, and ${\rho }_{x\eta }$ is the correlation between ${\eta }_{\mathcal{l}i}$ and a function of the defactored regressor ${\psi }_{\mathcal{l}i}^{(\mathcal{P})}$, where ${\psi }_{\mathcal{l}i}^{(1,2,3,4)}={\sum }_{p=1}^4{\tilde{z}}_{\mathcal{l}i,p}/\sqrt{4},{\psi }_{\mathcal{l}i}^{(p)}={\tilde{z}}_{\mathcal{l}i,p}$ for $p=1,\dots ,4$, and ${\tilde{z}}_{\mathcal{l}i,p}$ is the cross-sectionally standardized within average of the defactored regressor, $T^{-1}{\sum }_{t=1}^T{v}_{\mathcal{l}it}^p$, for $\mathcal{l}=1,\dots ,k$. In the article we report the results of designs 1–3, with the rest in the supplementary materials.

Design 1 corresponds to the model with homogeneous slopes (${\sigma }_{\eta }=0$) and Design 2 to the model with random coefficients (${\sigma }_{\eta }=0.5$ and ${\rho }_{x\eta }=0$). Designs 3-7 correspond to models with heterogeneous slopes that depend on the regressors (${\sigma }_{\eta }=0.5$ and ${\rho }_{x\eta }=0.5$) in different ways. These designs are considered in order to investigate the impact of such dependence on the estimators and the power properties of the LM test. In Design 3, the slope heterogeneity contains a fourth-degree polynomial of the defactored regressors (${\psi }_{\mathcal{l}i}^{(1,2,3,4)}$), which is highly nonlinear. In Designs 4–7, the slope heterogeneity contains ${\psi }_{\mathcal{l}i}^{(1)},{\psi }_{\mathcal{l}i}^{(2)},{\psi }_{\mathcal{l}i}^{(3)}$, and ${\psi }_{\mathcal{l}i}^{(4)}$, respectively. In the article we present the results for Designs 1 to 3, with the rest in the supplementary materials.

We consider all combinations of N = 50, 100, 200 and $T=25,50,100,200$. Throughout the experiment, we set ${\rho }_f=0.5,{\rho }_{\epsilon }=0.5$ and ${\rho }_v=0.5$. All tests are conducted at 5% significance level. All experimental results are based on 2000 replications.

4.4 Results

In order to compute the IPC estimator, it is necessary to estimate the number of factors in the error term, ${\mathbf{u}}_i$. The frequencies of the number of factors estimated by the ER and GR methods proposed by Ahn and Horenstein (2013) are reported in Table 2; see Section 4.1 for the estimation procedure. Note that, as summarized in Table 1, r is two in the models with homogeneous slopes (Design 1) and three in the models with random coefficients (Design 2). The results in Table 2 suggest that ER and GR estimate r precisely enough for all designs, though the ER estimator is slightly more accurate than the GR, especially when T is not large. Based on this result, the experimental results reported hereafter are based on the r estimated by the ER method.

Table 2 Frequency of the estimated number of factors in the error term by the ER and GR methods.

Design 1: homogeneous slopes (r = 2)							Design 2: random coefficients (r = 3)
Freq	ER			GR			ER			GR
T / N	50	100	200	50	100	200	50	100	200	50	100	200
$\widehat{r}=1$
25	1.0	0.3	0.2	3.6	1.5	1.0	3.2	1.4	0.7	13.5	7.9	5.0
50	0.0	0.0	0.0	0.2	0.0	0.0	0.4	0.0	0.0	3.1	0.6	0.1
100	0.0	0.0	0.0	0.0	0.0	0.0	0.1	0.0	0.0	0.6	0.1	0.0
200	0.0	0.0	0.0	0.0	0.0	0.0	0.1	0.0	0.0	0.3	0.0	0.0
$\widehat{r}=2$
25	98.9	99.7	99.8	96.5	98.5	99.1	16.6	10.3	6.7	23.4	17.9	14.4
50	100.0	100.0	100.0	99.9	100.0	100.0	6.5	1.4	0.5	11.7	4.3	1.7
100	100.0	100.0	100.0	100.0	100.0	100.0	2.3	0.1	0.0	6.4	0.6	0.0
200	100.0	100.0	100.0	100.0	100.0	100.0	1.7	0.0	0.0	4.2	0.1	0.0
$\widehat{r}=3$
25	0.1	0.0	0.0	0.0	0.0	0.0	80.1	88.3	92.7	63.1	74.3	80.7
50	0.0	0.0	0.0	0.0	0.0	0.0	93.2	98.7	99.6	85.2	95.2	98.3
100	0.0	0.0	0.0	0.0	0.0	0.0	97.7	100.0	100.0	93.1	99.4	100.0
200	0.0	0.0	0.0	0.0	0.0	0.0	98.3	100.0	100.0	95.5	99.9	100.0
$\widehat{r}=4$
25	0.0	0.0	0.0	0.0	0.0	0.0	0.2	0.0	0.0	0.0	0.0	0.0
50	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
100	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
200	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

NOTE: We set the maximum number of factors to six to estimate r (and $r_1+r_2$). The frequencies of $\widehat{r}=5$ and 6 were zero. ER and GR stands for the Eigenvalue Ratio and Growth Ratio methods proposed by Ahn and Horenstein (2013). See Table 1 for a summary of experimental designs.

Given the estimated number of factors, we can proceed to the estimation of ${\mathit{\beta }}^0$ and associated inferences based on the robust PHAC standard errors. In this experiment, we consider three different estimators of ${\mathit{\beta }}^0$: the IPC estimator without bias correction (IPC), the analytically bias-corrected IPC estimator (BC-IPC), and the jackknife bias-corrected IPC estimator (JK-IPC). Table 3 summarizes the performance of these IPC estimators for Designs 1–3. The “ER” signifies the dependence of the estimators on the ER estimation of r. In the table, the bias (Bias), standard deviation (SD), root mean square error (RMSE) of the estimates for ${\beta }_1^0$, and the size (in percent) of the robust Wald test for $H_0:{\beta }_1^0=1$ (Size), are reported. Bias, SD, and RMSE are all multiplied by 100. The results for ${\beta }_2^0$ are not reported because they are qualitatively very similar.

Table 3 Finite sample properties of the IPC estimators of ${\beta }_1^0$ and the associated robust Wald test with interactive effects, heteroscedastic and serially correlated idiosyncratic errors.

Design 1: homogeneous slopes $({\sigma }_{\eta =0},{\rho }_{x\eta }=0)$
	Bias($\times 100$)			SD($\times 100$)			RMSE($\times 100$)			Size(%)
T / N	50	100	200	50	100	200	50	100	200	50	100	200
IPC(ER)
25	–3.599	–1.674	–0.739	6.556	4.418	3.074	7.478	4.725	3.161	49.4	44.8	39.0
50	–4.069	–1.978	–0.992	4.485	2.878	1.974	6.056	3.492	2.209	57.5	46.6	39.6
100	–4.315	–2.085	–1.034	3.122	1.887	1.250	5.326	2.812	1.622	71.6	54.4	41.2
200	–4.721	–2.274	–1.113	2.323	1.424	0.996	5.262	2.683	1.494	87.9	68.5	51.7
BC-IPC(ER)
25	0.228	0.160	0.099	3.042	2.083	1.441	3.050	2.089	1.445	8.2	7.6	6.0
50	0.047	0.033	0.016	1.958	1.367	0.939	1.959	1.368	0.939	7.4	6.9	5.6
100	–0.003	–0.006	0.008	1.497	1.010	0.710	1.497	1.010	0.710	7.7	6.7	5.7
200	0.006	0.000	0.001	1.218	0.825	0.573	1.218	0.825	0.573	7.0	5.9	5.4
JK-IPC(ER)
25	0.096	0.093	0.038	3.506	2.274	1.580	3.507	2.276	1.580	11.8	11.0	9.2
50	0.015	0.007	0.002	2.156	1.480	1.005	2.156	1.480	1.005	9.3	9.6	7.6
100	–0.033	–0.014	0.004	1.590	1.061	0.746	1.590	1.061	0.746	9.4	7.3	6.4
200	–0.022	–0.004	–0.002	1.292	0.856	0.590	1.292	0.856	0.590	9.5	7.2	5.9
Design 2: heterogeneous slopes uncorrelated with $x_{\mathcal{l}it}({\sigma }_{\eta }=0.5,{\rho }_{x\eta }=0)$
	Bias($\times 100$)			SD($\times 100$)			RMSE($\times 100$)			Size(%)
T / N	50	100	200	50	100	200	50	100	200	50	100	200
IPC(ER)
25	–4.294	–2.282	–0.963	10.680	7.490	5.247	11.511	7.830	5.334	16.4	12.3	9.7
50	–4.693	–2.506	–1.143	9.069	6.257	4.322	10.212	6.741	4.471	17.6	11.7	8.7
100	–4.901	–2.611	–1.205	8.114	5.763	3.984	9.479	6.327	4.162	19.2	12.7	8.4
200	–5.371	–2.770	–1.311	7.666	5.487	3.824	9.360	6.147	4.042	20.4	13.1	7.9
BC-IPC(ER)
25	0.163	0.099	0.155	10.455	7.240	4.939	10.457	7.241	4.942	12.6	9.3	6.2
50	–0.202	–0.174	–0.017	9.073	6.150	4.219	9.076	6.153	4.219	12.8	8.6	7.1
100	–0.233	–0.197	–0.049	8.328	5.742	4.004	8.331	5.746	4.005	12.9	8.7	7.3
200	–0.263	–0.171	–0.050	7.933	5.555	3.860	7.937	5.557	3.861	12.8	9.1	6.5
JK-IPC(ER)
25	–0.414	–0.199	–0.077	12.167	7.994	5.314	12.174	7.997	5.315	18.0	11.9	9.0
50	–0.184	–0.204	–0.082	9.763	6.342	4.373	9.765	6.345	4.373	15.1	9.5	8.1
100	–0.100	–0.121	–0.060	8.726	5.865	4.062	8.726	5.867	4.062	14.5	9.7	7.4
200	–0.093	–0.083	–0.027	8.145	5.617	3.874	8.145	5.617	3.874	13.9	9.0	6.4
Design 3: heterogeneous slopes correlated with $x_{\mathcal{l}it}({\sigma }_{\eta }=0.5,{\rho }_{x\eta }=0.5,{\psi }_{\mathcal{l}i}^{(1,2,3,4)})$
	Bias($\times 100$)			SD($\times 100$)			RMSE($\times 100$)			Size(%)
T / N	50	100	200	50	100	200	50	100	200	50	100	200
IPC(ER)
25	1.070	3.570	4.964	11.075	7.527	5.137	11.127	8.331	7.144	15.3	12.8	18.6
50	1.066	3.746	5.217	9.357	6.377	4.164	9.418	7.396	6.675	13.5	13.3	22.6
100	0.882	3.735	5.256	8.228	5.596	3.568	8.275	6.728	6.353	11.5	12.5	23.1
200	0.778	3.876	5.493	7.879	5.248	3.320	7.917	6.524	6.418	11.7	13.5	26.9
BC-IPC(ER)
25	4.567	5.160	5.592	10.452	7.129	4.913	11.406	8.801	7.444	15.3	16.3	22.0
50	5.065	5.823	6.250	9.064	6.149	3.964	10.383	8.469	7.401	16.9	19.6	29.6
100	5.221	5.956	6.356	8.318	5.555	3.565	9.821	8.145	7.287	18.5	20.2	34.4
200	5.540	6.291	6.662	8.050	5.322	3.368	9.772	8.240	7.465	21.4	23.7	40.0
JK-IPC(ER)
25	4.798	5.451	5.789	12.302	8.062	5.462	13.205	9.732	7.959	22.3	20.5	26.0
50	5.937	6.444	6.572	9.981	6.522	4.141	11.613	9.168	7.768	21.7	24.1	33.7
100	6.118	6.548	6.656	8.793	5.734	3.634	10.712	8.704	7.583	23.4	25.1	38.2
200	6.560	6.909	6.981	8.521	5.404	3.398	10.754	8.772	7.764	26.4	27.7	43.3

First, we discuss the results for Design 1 (the model with homogeneous slopes). As can be seen in the table, the IPC estimator (without bias-correction) suffers from a large bias. As a result, the size of the robust Wald test deviates significantly from the nominal level. In contrast, the BC-IPC estimator with the analytical bias correction and the JK-IPC estimator with the jackknife bias correction successfully reduce the bias. The BC-IPC estimator and the JK-IPC estimator perform comparably. However, when the sample size is small (e.g., T = 25 or N = 50), the jackknife method can reduce the bias more effectively than the analytical bias correction. Meanwhile, the jackknife bias correction is likely to lead to greater sampling variability than the analytical bias correction. The standard deviation of the JK-IPC estimator is uniformly larger than the standard deviation of the BC-IPC estimator. As a result, the BC-IPC estimator uniformly outperforms the JK-IPC estimator in terms of the mean squared error criterion.¹⁰ Reflecting this small sample bias, there is a mild size distortion in the robust Wald test for small T or small N, which tends to disappear as N and T increase.

The results for Design 2 (random coefficient models) in Table 3 are qualitatively similar to those for Design 1. However, reflecting the fact that the IPC estimators are less efficient for random coefficient models, there is a moderate small sample bias in the BC-IPC and JK-IPC estimators, which dissipates as N increases. Due to this bias, there is a moderate size distortion in the robust Wald test for small N, but it quickly decreases as the sample size increases.

Next, let us look at the results for Design 3, in which the heterogeneity of the slopes depends on the regressors. The bias of the IPC, BC-IPC, and JK-IPC estimators is large and persistent for any sample size. As a consequence, the size distortions of the robust Wald tests are large and increase as the sample size rises. This implies the importance of the LM test for correlated random coefficients in our robust approach.

Table 4 reports the size and power of the LM test for the dependence of the random coefficients on the regressors. We note that in our alternative model (Design 3), the slope heterogeneity contains a fourth-degree polynomial of the defactored regressors (i.e., ${\psi }_{\mathcal{l}i}^{(1,2,3,4)}$), which is highly nonlinear. In this experiment, we consider ${\text{LM}}_{\text{CRC}}^{(1)}$ and ${\text{LM}}_{\text{CRC}}^{(2)}$ tests, which are defined in equation (20)(20) $${\text{LM}}_{\text{CRC}}^{\left(g\right)}={\widehat{\mathbf{s}}}^{†\prime }{\widehat{\mathbf{\Sigma }}}_{ss}^{-1}{\widehat{\mathbf{s}}}^{†}$$(20) . Recall that we recommend the test using the ${\text{LM}}_{\text{CRC}}^{(2)}$ statistic, because it is expected to have nontrivial power against a sufficiently broad class of alternative models, while the ${\text{LM}}_{\text{CRC}}^{(1)}$ test may not be sufficiently powerful against certain alternatives (see discussions in Section 3). The results reported in Table 4 show that the size of the ${\text{LM}}_{\text{CRC}}^{(1)}$ and ${\text{LM}}_{\text{CRC}}^{(2)}$ tests is very close to the nominal level for Design 1 (the model with homogeneous slopes) and Design 2 (random coefficient models). In the results for Design 3 (correlated random coefficients), we find that the power of the ${\text{LM}}_{\text{CRC}}^{(2)}$ test is high and rises as N and T increase, while the ${\text{LM}}_{\text{CRC}}^{(1)}$ test is much less powerful. Table C.3 in the supplementary materials reports the power of the ${\text{LM}}_{\text{CRC}}^{(1)}$ and ${\text{LM}}_{\text{CRC}}^{(2)}$ tests for Designs 4–7. Therein it is found that the ${\text{LM}}_{\text{CRC}}^{(2)}$ test remains as powerful across Designs 4–7 as in Design 3, whilst the ${\text{LM}}_{\text{CRC}}^{(1)}$ test is substantially less powerful in Design 5 (${\psi }_{\mathcal{l}i}^{(2)}$) and Design 7 (${\psi }_{\mathcal{l}i}^{(4)}$). These findings are in line with the results of our asymptotic local power analysis and they support our claim that the ${\text{LM}}_{\text{CRC}}^{(2)}$ test is powerful against a broad class of alternative models of correlated coefficients.

Table 4 Size and power of the ${\text{LM}}_{\text{CRC}}^{(1)}$ and ${\text{LM}}_{\text{CRC}}^{(2)}$ tests of correlated random coefficient.

	Size (Design 1)			Size (Design 2)			Power (Design 3)
	${\sigma }_{\eta }=0,{\rho }_{x\eta }=0$			${\sigma }_{\eta }=0.5,{\rho }_{x\eta }=0$			${\sigma }_{\eta }=0.5.{\rho }_{x\eta }=0.5$
T,N	50	100	200	50	100	200	50	100	200
${\text{LM}}_{\text{CRC}}^{(1)}$
25	2.6	2.6	3.1	2.7	2.4	2.4	5.9	7.1	9.7
50	3.4	3.3	2.3	2.5	3.0	3.3	8.9	12.0	16.8
100	3.3	2.0	2.6	2.7	2.3	2.3	12.1	16.9	28.8
200	3.4	2.9	2.3	2.7	3.0	2.4	18.8	32.2	51.4
${\text{LM}}_{\text{CRC}}^{(2)}$
25	2.9	3.5	3.9	3.4	3.7	3.6	9.9	18.3	34.2
50	2.9	2.7	3.2	4.4	3.9	3.8	22.7	48.8	77.8
100	2.9	2.8	3.4	4.4	3.2	3.0	45.3	83.2	97.2
200	2.9	3.0	2.9	4.5	3.7	3.7	68.4	97.7	100.0

NOTE: For Designs 1–3, see Table 1. The ${\text{LM}}_{\text{CRC}}^{(g)}$ is the proposed LM test of correlated random effects defined by (20), based on the BC-IPC estimator. The test statistics are referred to the 95% quantile of ${\chi }_g^2$ distribution. All results are based on 2000 replications. The number of factors is estimated by the ER method.

5 An Empirical Application: Feldstein–Horioka Puzzle

In this section, we present an empirical illustration of the proposed method. Specifically, we analyze the Feldstein–Horioka (F–H) puzzle, which has attracted recent attention in the literature.

Feldstein and Horioka (1980) found that in OECD countries, the long-run average of the national savings rate is significantly correlated with the long-run average of the domestic investment rate. This empirical finding seems to conflict with macroeconomic theory, which predicts that the determinants of saving and investment are different and that, assuming perfect international capital mobility, the investment decisions in a country should not be constrained by domestic saving. This is known as the F–H puzzle, and is included as one of the six major puzzles in international macroeconomics by Obstfeld and Rogoff (2000).

Several econometric methods have been applied to investigate the F–H puzzle, including instrumental variable methods, cointegration regressions and the standard fixed effect models.¹¹ Recently, Giannone and Lenza (2010) and Ginama, Hayakawa, and Kanemi (2018) analyzed the F–H puzzle by considering a panel regression model with interactive effects.

In this section, we apply our approach to analyze the F–H puzzle by using the OECD panel dataset covering 1968–1996 with N = 24 and $T=29,$ as used in Ginama, Hayakawa, and Kanemi (2018). Specifically, we estimate the following model:$$\left(\frac{I_{it}}{Y_{it}}\right)={\beta }_i\left(\frac{S_{it}}{Y_{it}}\right)+{\mathit{\lambda }}_i^{0\prime }{\mathbf{f}}_t^0+{\epsilon }_{it},$$where I_it , S_it , and Y_it denote investment, gross savings, and GDP for country i and year t, respectively. This specification is more general than considered in the literature in that the coefficient is allowed to be heterogeneous and interactive effects are included.

We now describe how we proceed when using our approach. Our interest lies in estimating the population mean of β_i . First, we need to determine the number of factors in the error term for the model in which the slope heterogeneity (if any) is ignored. Here, we follow the procedure described in Section 4.1. Having estimated the number of factors in the error term, we conduct the LM test for correlated random coefficients, which is proposed in Section 3. If the null hypothesis is not rejected, we proceed to statistical inference using the robust PHAC standard errors; otherwise we need to consider an alternative approach, as described at the end of Section 3.

Table 5 provides the estimation results. As a comparison, we computed a jackknife bias-corrected CCEP estimator proposed by Westerlund (2018) and the associated nonparametric standard error of (Pesaran 2006, eq. (69)) for random coefficient models.¹² With regard to the number of factors, both ${\widehat{r}}_{GR}$ and ${\widehat{r}}_{ER}$ suggest that the number of factors is one. The LM test for correlated random coefficients, ${\text{LM}}_{\text{CRC}}^{(2)}$, does not reject the null hypothesis of random coefficients or homogeneous slopes, which validates our robust approach. Note that the PHAC standard error is asymptotically valid for the models with homogeneous slopes or random coefficients. Meanwhile, the nonparametric standard error used for the CCEP estimator is developed primarily for the models with random coefficients, and a different standard error may have to be employed for the models with homogeneous slopes; see Theorem 4 in Pesaran (2006). Nevertheless, the BC-IPC and JK-CCEP estimation results show that the savings-GDP ratio is highly significant in the regression of the investment-GDP ratio, indicating the existence of the F–H puzzle.

Table 5 Estimation results for the Feldstein–Horioka puzzle based on OECD data covering 1968–1996, (N = 24, T = 29).

Estimated number of factors
$\widehat{r}$	${\widehat{r}}_{,ER}=1$	${\widehat{r}}_{,GR}=1$
LM test for correlated random coefficient
	test statistic	p-value
${\text{LM}}_{\text{CRC}}^{(2)}$	0.5101	0.7749
Estimation results with $\widehat{r}=1$
	BC-IPC	JK-CCEP
Coef.	0.4238	0.3518
S.E.	0.0891	0.1193

NOTE: BC-IPC is the bias-corrected IPC estimator with the PHAC standard error. JK-CCEP is the jackknife bias-corrected CCEP estimator of Westerlund (2018) with the nonparametric standard error of Pesaran (2006) for random coefficients. For the estimation of the number of factors and computation of BC-IPC, data are demeaned in both dimensions while undemeaned data is used for JK-CCEP.

6 Concluding Remarks

In this article, we have proposed a robust approach against heteroscedasticity, error serial correlation and slope heterogeneity in linear models with interactive effects for large panel data. First, consistency and asymptotic normality of the pooled iterated principal component (IPC) estimator for the models with random coefficients and homogeneous slopes have been established. Then, we have proved the asymptotic validity of the associated Wald test for slope parameter restrictions based on the panel heteroscedasticity and autocorrelation consistent (PHAC) variance matrix estimator for the models with random coefficients and homogeneous slopes, which does not require the Newey-West type time-series parameter truncation. These results asymptotically justify the use of the same pooled IPC estimator and the PHAC standard error for both homogeneous-slope and heterogeneous-slope models. This robust approach can significantly reduce the model selection uncertainty for applied researchers. In addition, we have proposed a Lagrange Multiplier (LM) test for correlated random coefficients with covariates. This test has nontrivial power against correlated random coefficients, but not for random coefficients and homogeneous slopes. The LM test is important because the IPC estimator becomes inconsistent with correlated random coefficients. The finite sample evidence and an empirical application support the usefulness and reliability of our robust approach.

We have examined finite sample performance of the estimators, the tests of linear restrictions on the slope parameters, and the LM tests for correlated random coefficients. We have examined the finite sample performance of the estimators, the tests of parameter restrictions using the PHAC variance estimator, and the LM test for correlated random coefficients. The results show that the size of the proposed robust Wald test with the bias-corrected IPC estimator is sufficiently close to the nominal level in both slope homogeneity and slope heterogeneity, and that the LM test for correlated random coefficients has correct size under both slope homogeneity and random coefficients while exhibiting high power for correlated random coefficients. The finite sample evidence, together with the empirical application for the Feldstein–Horioka puzzle, support the reliability and the usefulness of the proposed robust approach.

It is well recognized that bootstrapping in general can provide more accurate inference; see, for example, Hall (1992). Recently, there has been a rapid growth in the study of bootstrap methods for factor models; see Gonçalves and Perron (2014, 2020), among others. This is an interesting research direction, as applying bootstrap methods to panel data models with interactive effects can lead to more precise inferences.

The proposed approach investigated for the IPC estimator is most likely to be applicable for pooled estimators in other approaches, including the Common Correlated Effects (CCE) estimator of Pesaran (2006), the principal component (PC) estimator investigated by Westerlund and Urbain (2015), and the two-step instrumental variable (2SIV) estimator proposed by Cui et al. (2021), among others. In particular, for these estimators the factors in the regressors and in the error term are extracted, whilst the IPC estimator only exploits the information of the factors in the error term. Studying such applications seems to be an interesting extension.

As emphasized in the article, when the null hypothesis of random coefficient models is rejected in favor of the alternatives, it is preferable to employ estimators which are consistent for the models with correlated random coefficients. In this situation, the mean group estimators based on the CCE, 2SIV and ML approaches proposed by Pesaran (2006), Cui et al. (2021), and Li, Cui, and Lu (2020), among others, seem to be possible choices. Investigating inferential methods for such estimators under correlated random coefficients seems to be an intriguing future research theme.

Acknowledgments

The authors thank Christian Hansen (the editor), the associate editor, and two anonymous referees for their valuable comments and suggestions. We are grateful to Milda Norkutė, Joakim Westerlund and the participants of the 2019 Asia Meeting of the Econometric Society for their helpful comments on the earlier version of this article.

Supplementary Materials

The supplementary material consists of online appendices that discuss the experimental results with proofs of key results and additional discussion, and computational codes that replicate the experimental and empirical results.

Additional information

Funding

This work was supported by JSPS KAKENHI (grant numbers 20H01484, 20H05631, 20K20760, 21H00700, and 21H04397).

Notes

1 In the fixed T panel data analysis, the PHAC estimator is proposed by Arellano (1987) and widely used in the related literature.

2 In this article, we propose to use the PHAC estimator for the $\sqrt{N}$-consistent estimator, which is different from Pesaran (2006) and Reese and Westerlund (2018). See Section 2.

3 One important exception is Ando and Bai (2017), in which shrinkage estimators for sparse heterogeneous slopes in the models with interactive effects are considered.

4 In this article, we do not consider the model with the group-wise homogeneous slope model, which is considered by Su and Ju (2018).

5 See Galvao and Kato (2014) for a similar discussion in misspecified dynamic panel data models.

6 We do not consider cross-sectional and/or time-series structural breaks in ${\mathit{\beta }}_i$ which is beyond the scope of this article.

7 For testing $H_0:{\theta }_p=0$ for $p=1,\dots ,\infty$, a nonparametric regression approach (e.g., sieve or spline regressions) could be employed.

8 Alternatively, we may use the criteria proposed by Bai and Ng (2002).

9 We also considered the case in which ${\varpi }_{\mathcal{l}it}$ is generated from normal distribution. The results are not reported because they were qualitatively very similar.

10 A similar superior performance of analytical bias-correction over jackknife bias-correction is found in Fernández and Weidner (2016).

11 See sec. 1 of Ginama, Hayakawa, and Kanemi (2018) for a brief literature survey.

12 For the IPC estimation, the data are within-transformed and cross-sectionally demeaned, while for the CCEP estimation the data are within transformed only.