Estimation and Inference for Three-Dimensional Panel Data Models

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

85,408 characters · 11 sections · 46 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

{7pt} {7pt}

titlepage\begin{center} { \bf Estimation and Inference for \\ Three-Dimensional Panel Data Models \begingroup \footnote{ \\ An earlier version of this paper has been circulated and submitted to working paper series under the title “Multi-Level Panel Data Models: Estimation and Empirical Analysis", which is the Working Paper \#4/22 available at https://ideas.repec.org/p/msh/ebswps/2022-4.html. We thank the participants of 2024 International Workshop on Econometrics with Its Application and Practice in Finance at Nankai University for constructive suggestions and comments. Gao and Peng would like to acknowledge the Australian Research Council Discovery Projects Program for its financial support under Grant Numbers: DP200102769 & DP210100476. Liu's research is financially supported by the National Natural Science Foundation of China under Grant No. 72203114.{ • $^{\ast}$University of North Texas, United States. \\ $^\dag$Monash University, Australia. \\ $^\sharp$Nankai University, China. } } \addtocounter{footnote}{-1} \endgroup } {\sc Guohua Feng$^{\ast}$, Jiti Gao$^{\dag}$, Fei Liu$^\sharp$ and Bin Peng$^{\dag}$} \today \begin{abstract} Hierarchical panel data models have recently garnered significant attention. This study contributes to the relevant literature by introducing a novel three-dimensional (3D) hierarchical panel data model, which integrates panel regression with three sets of latent factor structures: one set of global factors and two sets of local factors. Instead of aggregating latent factors from various nodes, as seen in the literature of distributed principal component analysis (PCA), we propose an estimation approach capable of recovering the parameters of interest and disentangling latent factors at different levels and across different dimensions. We establish an asymptotic theory and provide a bootstrap procedure to obtain inference for the parameters of interest while accommodating various types of cross-sectional dependence and time series autocorrelation. Finally, we demonstrate the applicability of our framework by examining productivity convergence in manufacturing industries worldwide. {\em Keywords}: Asymptotic Theory, Bias Correction, Dependent Wild Bootstrap, Hierarchical Model {\em JEL classification}: C23, O10, L60 \end{abstract} \end{center}

Introduction

\setcounter{equation}{0}

In the past three decades or so, panel data research has undergone significant advancements, leveraging its capacity to extract richer insights from both cross-sectional and time dimensions. Within this domain, a notable strand of research (Bai,FanLiaoWang) has garnered increasing attention, with a particular focus on addressing cross-sectional dependence (CSD) in large $N$ and $T$ settings by accounting for correlations among observations across different entities within a panel. Specifically, strong CSD is often addressed through a factor structure approach involving the estimation of a low-rank representation, while weak CSD poses inference challenges related to estimating the asymptotic covariance matrix. Additionally, accounting for time series autocorrelation (TSA) is essential when analyzing data with a temporal dimension. To tackle both CSD and TSA jointly, various methods have been proposed in the relevant literature. For instance, goncalves_2011 suggests using the moving block bootstrap method, BAI2020 employ a thresholding method in conjunction with a heteroskedasticity and autocorrelation-consistent (HAC) covariance matrix estimator, and GPY2023 propose a dependent wild bootstrap procedure by extending the time series approach of shao2010 to a two--dimensional panel data framework.

Recent developments in panel data research, especially in contexts where observations exhibit hierarchical structures, have led to the emergence of the relevant literature on hierarchical panel data models (Matyas, KSS2020, CYZ2022, Zhang2023, JLS2023, JLS2024). This evolution reflects the recognition of the necessity to extend beyond traditional panel data models to capture dependencies at multiple levels of aggregation, such as individuals nested within groups or regions nested within countries. Consequently, the literature on hierarchical panel data models builds upon and extends methodologies developed in the cross-sectional dependence literature to accommodate these hierarchical structures, ultimately enabling researchers to more accurately model and analyze complex panel datasets.

Another pertinent literature is regression with incidental parameters, a concept pioneered by Neyman. LANCASTER2000391 commends their seminal paper as a remarkable contribution, widely acknowledged for its influence within the statistical community. To see the challenges posed by incidental parameters in the context of 3 dimensional (3D) panels, we present the following table outlining the potential incidental parameters across different data structures:

table[table omitted — 475 chars of source]

While a 3D panel involves only one additional layer, it significantly increases complexity. To the best of our knowledge, no prior studies have attempted to integrate all potential specifications of incidental parameters into a single regression model. Hence, a unified framework is needed. Fortunately, incidental parameters can be integrated using a factor structure. For example, all four conceivable scenarios of a 2D panel can fit within a structure such as $ \pmb{\alpha}_i^\top \mathbf{g}_t$. This approach offers valuable insights into addressing the challenge of 3D panel data modeling.

With the above challenges in mind, the primary objective of this study is to contribute to the literature on hierarchical panel data models by introducing a panel data regression model with three sets of latent factor structures: one set of global factors and two sets of local factors. Compared with pure factor models, this model is capable of estimating the coefficient of a particular observed independent variable while accounting for both commonalities and individual differences across entities. In addition, in comparison with panel data regression models with only one or two sets of factors, the proposed model and estimation theory is capable to model dependencies that exist at multiple levels of aggregation. This model finds applications in various contexts. For example, in the context of modeling economic convergence, the latent global factors represent economic events affecting the growth of all countries and industries, while the two local factors represent events affecting specific countries or industries, respectively. In the context of modeling international trade, the latent global factors represent factors affecting the trade of all countries, while the two local factors represent events affecting importing or exporting countries, respectively.

To better understand the 3D factor structure, suppose that we have two groups of individuals:

eqnarray[eqnarray omitted — 55 chars of source]

where $L$ and $N$ are two positive integers, and $[S]\coloneqq \{1,\ldots, S\}$ for any positive integer $S$. In the literature of modeling empirical growth (e.g., Rodrik2013), $i$ represents industries, and $j$ denotes countries; in the context of bilateral trade (e.g., CYZ2022), $i$ and $j$ respectively represent importing countries and exporting countries; in the stochastic actor-oriented models (e.g., KS2023), $i$ refers to social groups and $j$ refers to individuals within each group; etc. In addition, we let $y_{ijt}$ denote a set of outcome variables:

eqnarray[eqnarray omitted — 68 chars of source]

which is indexed by both $i$ and $j$, and also varies over each time period $t$. The definition of $y_{ijt}$ should be self-evident in the aforementioned studies, so we omit these details here. Throughout this paper, we assume that each layer of information is influenced by distinct factors, enabling us to extract essential insights through careful data processing. Our approach diverges in that we aim to disentangle factors layer by layer and across different blocks within each layer, rather than consolidating factors from disparate nodes.

For a clearer grasp of the 3D factor structure, we visually represent the three sets of latent factors ($\mathbf{f}_{t}$, $\mathbf{f}_{it}^\circ$, and $\mathbf{f}_{jt}^\bullet$) in Figure (ref), utilizing empirical growth data as a demonstration.

{

figure[figure omitted — 1,009 chars of source]

} Specifically, $\mathbf{f}_{t}$, $\mathbf{f}_{it}^\circ$ and $\mathbf{f}_{jt}^\bullet$ represent $\ell\times 1$, $\ell_i^\circ\times 1$, and $\ell_j^\bullet\times 1$ vectors respectively, denoting the global, country-specific, and industry-specific factors or shocks driving $y_{ijt}$. Here, $\ell$, $\ell_i^\circ$, and $\ell_j^\bullet$ are non-negative fixed integers, subject to the constraint:

eqnarray[eqnarray omitted — 136 chars of source]

where $c$ is a fixed positive constant. This condition implies that $\ell$, $\ell_i^\circ$, and $\ell_j^\bullet$ can potentially be 0. In an extreme scenario, if $y_{ijt}$ behaves entirely randomly like an idiosyncratic error, then $\ell=\ell_i^\circ=\ell_j^\bullet=0$ for all $(i,j)$ pairs. Figure (ref) illustrates that $f_t$ influences every industry and country, while each country-specific (or industry-specific) shock can also impact a specific subset of countries or industries. Mathematically, for all $t\in [T]$, Figure 1 reveals the following mapping:

{

eqnarray[eqnarray omitted — 1,078 chars of source]

where $\pmb{\gamma}_{ij}$, $\pmb{\gamma}_{ij}^\circ$, and $\pmb{\gamma}_{ij}^\bullet$} are factor loadings indicating how each factor influences $y_{ijt}$ at each time period $t$. Figure (ref), along with the restriction (ref) and the mapping (ref), together entail the following observations:

itemize[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • Figure (ref), the restriction (ref), and the mapping (ref) encompass all possible combinations of incidental parameters. • The global shock $\mathbf{f}_t$ (if exists) affects all countries and industries, while the country and industry shocks ($\mathbf{f}_{it}^\circ$ and $\mathbf{f}_{jt}^\bullet$) respectively impact a specific subset of countries or industries. • From the a signal-to-noise ratio point of view, $\mathbf{f}_t$ can be recovered by aggregating all available information, whereas $\mathbf{f}_{it}^\circ$ and $\mathbf{f}_{jt}^\bullet$ can be identified block by block. • The dimensions $i$ and $j$ exhibit symmetry, despite potential differences in the values of $L$ and $N$.

Finally, to accommodate potential endogeneity, we posit that the unobservable factors exert influence on regressors, subject to less restrictive conditions (refer to Remark (ref) and Assumption (ref) for more details). Further elaboration on this aspect will be provided when detailing the model setup and associated conditions in Section (ref). Our objective is twofold: to estimate the coefficients of these endogenous regressors and to discern factors across various levels or blocks of the hierarchy.

Concluding this section, we summarize our contributions as follows. In addition to the introduction of a panel data regression model featuring three sets of latent factor structures, this study also provides additional contributions in the following areas:

itemize[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • We develop an estimation method to recover the parameters of interest and unobservable factors across levels and blocks of the hierarchy, and establish its asymptotic properties. The unified hierarchy structure in a 3D framework covers the total 64 possibilities of incidental parameters which have been the centre of a wide range of applications since Neyman. • The main model and its estimation method discussed in Section 2 below are developed for the case where both cross--sectional dependence and serial correlation are allowed for the error terms. By contrast, the existing literature (see, JLS2024, for example) focuses on the case where the error components are independent and identically distributed (i.i.d.). • We propose a bootstrap procedure to obtain inferences for the parameters of interest, accommodating various types of cross--sectional dependence (CSD) and time series autocorrelation (TSA) within the data generating process of the hierarchy. • To validate our theoretical findings, we conduct extensive simulations and analyze real data examples. In the empirical study, we specifically utilize data from manufacturing industries at the ISIC two-digit level to examine the twin hypotheses of conditional and unconditional convergence for manufacturing industries across countries.

The rest of the paper is organized as follows. Section (ref) presents our model and methodology. An asymptotic theory is established accordingly for each step involved in different estimation procedures. Section (ref) conducts extensive numerical studies to examine the theoretical findings. Specifically, Section (ref) examines the theoretical findings using extensive simulations, while Section (ref) uses a set of data from manufacturing industries at the ISIC two-digit level to examine the twin hypotheses of conditional and unconditional-convergence for manufacturing industries across countries. Section (ref) concludes. In the online supplementary Appendices A and B, Appendix A lists the necessary tables and pictures for Section 4. Appendix (ref) provides two detailed numerical implementations; Appendix (ref) proposes a Jackknife based bias correction method; Appendix (ref) includes the necessary preliminary lemmas; Appendix (ref) presents the detailed proofs; some additional estimation results of the empirical study are reported in Appendix B.

Before proceeding further, it is convenient to introduce some notations that will be repeatedly used throughout the article. Vectors and matrices are always written in bold font. For a matrix $\mathbf{A}$, $\|\mathbf{A}\|$ and $\|\mathbf{A}\|_2$ denote the Frobenius norm and the spectral norm of $\mathbf{A}$, respectively, and $\mathbf{A}^\top$ stands for the transpose of $\mathbf{A}$. Provided that $\mathbf{A}$ has full column rank, let $\mathbf{M}_{\mathbf{A}}=\mathbf{I}- \mathbf{P}_{\mathbf{A}}$ with $\mathbf{P}_{\mathbf{A}} =\mathbf{A}(\mathbf{A}^\top \mathbf{A})^{-1}\mathbf{A}^\top$. For two block wise matrices $\mathbf{A}=\{ \mathbf{a}_{ij}\}$ and $\mathbf{B}=\{ \mathbf{b}_{ij}\}$, provided $\mathbf{a}_{ij}\mathbf{b}_{ij}$ is well defined, we let $\circ$ define the block wise Hadamard product, i.e., $\mathbf{A}\circ \mathbf{B} =\{ \mathbf{a}_{ij}\mathbf{b}_{ij}\}$. For each $\mathbf{a}_{ij}$, let $\sum_{i,j} \mathbf{a}_{ij}=\sum_{i=1}^L\sum_{j=1}^N \mathbf{a}_{ij} $. Also, we let $\operatorname*{\normalfont\textrm{diag}}\{\mathbf{A}, \mathbf{B} \}$ return a diagonal matrix with $\mathbf{A}$ and $\mathbf{B}$ on the main diagonal. For a vector $\mathbf{b}$, let $\|\mathbf{b}\|_1$ define the $L^1$ norm. For two scalars $m$ and $n$, $m\wedge n=\min\{m, n\}$, $m\vee n=\max\{m, n\}$. For two random variables $a$ and $b$, $a\asymp b$ stands for $a=O_P(b)$ and $b=O_P(a)$. $I(\cdot)$ represents the conventional indicator function, and “$\to_P$" and “$\to_D$" stand for convergence in probability and convergence in distribution, respectively. $E^*$ and $\text{Pr}^*$ stand for the expectation and probability induced by the bootstrap procedure.

Model, Assumptions and Estimation Method

\setcounter{equation}{0}

In this section, we will introduce our panel data regression model featuring three sets of latent factor structures, outline a method for estimating the model, and establish its asymptotic properties. Before delving into the specifics of our model, we wish to emphasize two crucial points. Firstly, our main focus is on regressing $y_{ijt}$ on a $d\times 1$ vector $\mathbf{x}_{ijt}$, where $d$ is fixed. Secondly, as discussed in the introduction, we hypothesize that the unobservable factors impact the regressors to address potential endogeneity. With this consideration, our panel data regression model with three sets of latent factor structures is formulated as follows:

eqnarray[eqnarray omitted — 557 chars of source]

where $\pmb{\beta}_{ij}$ and every quantity on the right hand side of (ref) are unknown, and $(i,j,t)\in [L]\times [N]\times [T]$ with $(L,N,T)$ can all be large.

Our goal is to infer $\pmb{\beta}_{ij}$, and to recover the spaces spanned by the three sets of unobservable factors. For the time being, we assume that $\ell$, $\ell_i^\circ$'s and $\ell_j^\bullet$'s are known, and we will address their estimation in Section (ref). The setup connects with the existing works, such as Ando, KSS2020 and JLS2023, by addressing many challenges of the relevant literature within one framework. Meanwhile, it is worth mentioning that the paper by JLS2024 has an almost identical setup as our model (ref), and they offer a set of tests on the heterogeneous coefficients which are driven by the independent and identically distributed (i.i.d.) error components attached to the coefficients. In this paper, we pay particular attention to the case where the random errors of the main model are allowed to be both cross--sectionally dependent and serially correlated.

We now rewrite the main equation (i.e., the upper equation) of (ref) in vector form as follows:

eqnarray[eqnarray omitted — 263 chars of source]

where { $\mathbf{Y}_{ij\centerdot} =(y_{ij1},\ldots, y_{ijT})^\top$, $\mathbf{X}_{ij\centerdot} =(\mathbf{x}_{ij1},\ldots, \mathbf{x}_{ijT})^\top$, $\pmb{\mathcal{E}}_{ij\centerdot} =(\varepsilon_{ij1},\ldots, \varepsilon_{ijT})^\top$, $\mathbf{F} =(\mathbf{f}_1,\ldots, \mathbf{f}_T)^\top$}, $\mathbf{F}_i^\circ =(\mathbf{f}_{i1}^\circ,\ldots, \mathbf{f}_{iT}^\circ)^\top$, and $\mathbf{F}_j^\bullet =(\mathbf{f}_{j1}^\bullet,\ldots, \mathbf{f}_{jT}^\bullet)^\top$.

remarkWe refrain from expressing $\mathbf{x}_{ijt}$ in matrix form at this stage, as we aim to minimize constraints on $\pmb{\phi}_{ij}$, $\pmb{\phi}_{ij}^\circ$ and $\pmb{\phi}_{ij}^\bullet$ as much as possible. These parameters may even exhibit characteristics akin to weak factor loadings, as observed in Yamagata2023. If the focus is about the structure of $\mathbf{x}_{ijt}$ only, we refer the interested reader to JLS2023 for a comprehensive investigation on the factor structure without regressors.

According to (ref) and (ref), it is reasonable to assume that for $\forall (i,j)$

eqnarray[eqnarray omitted — 104 chars of source]

has full column rank. Otherwise, one can always reorganize the unobservable factors to achieve a representation with full column rank.

remarkIf $\mathbf{F}_{ij}^*$ is observable, we can immediately rewrite (ref) as \begin{eqnarray} \mathbf{M}_{\mathbf{F}_{ij}^*}\mathbf{Y}_{ij\centerdot}=\mathbf{M}_{\mathbf{F}_{ij}^*}\mathbf{X}_{ij\centerdot} \pmb{\beta}_{ij} +\mathbf{M}_{\mathbf{F}_{ij}^*}\pmb{\mathcal{E}}_{ij\centerdot} , \end{eqnarray} which enables consistent estimation for $\pmb{\beta}_{ij}$ using the OLS method. For the case with unobservable $\mathbf{F}_{ij}^*$, a very intuitive idea following a typical two dimensional panel data approach (such as Ando) is to consider an objective function in the following form: \begin{eqnarray} (\mathbf{Y}_{ij\centerdot}-\mathbf{X}_{ij\centerdot} \mathbf{b}_{ij} )^\top \mathbf{M}_{\mathbf{C}_{ij}} (\mathbf{Y}_{ij\centerdot}-\mathbf{X}_{ij\centerdot} \mathbf{b}_{ij} ), \end{eqnarray} where each $\mathbf{C}_{ij} = (\mathbf{C}, \mathbf{C}_{i}^\circ, \mathbf{C}_{j}^\bullet)$ is a generic matrix satisfying $\frac{1}{T}\mathbf{C}_{ij}^\top \mathbf{C}_{ij} =\mathbf{I}_{\ell+\ell_i^\circ+\ell_i^\bullet}$. One would hope to obtain the estimates of $\pmb{\beta}_{ij}$ and $\mathbf{F}_{ij}^*$ by combining the quantity of (ref) over $(i,j)$. However, the restriction $\frac{1}{T}\mathbf{C}_{ij}^\top \mathbf{C}_{ij} =\mathbf{I}_{\ell+\ell_i^\circ+\ell_i^\bullet}$ simply cannot be fulfilled practically. To illustrate this, suppose we consider a fixed $i$. Then $(\mathbf{C}, \mathbf{C}_{i}^\circ)$ allow us to generate a set of $\mathbf{C}_{|i}^\bullet =\{\mathbf{C}_{j}^\bullet\}$. Once moving on to a different $i^*$, $(\mathbf{C}, \mathbf{C}_{i^*}^\circ)$ will generate another set of $\mathbf{C}_{|i^*}^\bullet =\{\mathbf{C}_{j}^\bullet\}$. In general, $\mathbf{C}_{|i}^\bullet $ and $\mathbf{C}_{|i^*}^\bullet $ are not necessarily the same unless (1): $(L\vee N)\ll T$; and (2): $\{\mathbf{C}_{i}^\circ\}$ and $\{\mathbf{C}_{j}^\bullet\}$ are generated from two orthogonal spaces from the perspective of vector multiplication. However, it rules out the most common restriction, such as $L \asymp N \asymp T$.

To tackle the problem raised in Remark (ref), we put our thoughts in a nutshell below. For simplicity, suppose that as $T\rightarrow \infty$

eqnarray[eqnarray omitted — 160 chars of source]

which does not lose any generality, as it is only a matter of rotating matrices in view of the factor structure admitting the following expression:

eqnarray[eqnarray omitted — 178 chars of source]

in which $\pmb{\gamma}_{ij}^* = (\pmb{\gamma}_{ij}^\top, \pmb{\gamma}_{ij}^{\circ\top}, \pmb{\gamma}_{ij}^{\bullet\top} )^\top$. Lemme (ref) of the online Appendix A shows the feasibility of (ref) when taking $\max$ over $(i,j)$. By (ref), simple algebra shows that

eqnarray[eqnarray omitted — 273 chars of source]

in which the representation $\mathbf{I}_T-\frac{1}{T}\mathbf{F}\mathbf{F}^\top -\frac{1}{T}\mathbf{F}_i^\circ\mathbf{F}_i^{\circ \top} - \frac{1}{T}\mathbf{F}_j^\bullet\mathbf{F}_j^{\bullet\top}$ automatically separates $\mathbf{F}$, $\mathbf{F}_i^\circ$, and $\mathbf{F}_j^\bullet$, and also serves as a projection matrix asymptotically. Thus, (ref) sheds light on how to recover $\mathbf{F}$, $\mathbf{F}_i^\circ$, and $\mathbf{F}_j^\bullet$ within a single framework. We are now ready to proceed.

Estimation Procedure

We start by introducing a set of new symbols. Let $\mathbf{C}$, $\mathbf{C}_i^\circ$, and $\mathbf{C}_j^\bullet$ be some generic matrices each sharing the same dimensions with $\mathbf{F}$, $\mathbf{F}_i^\circ$, and $\mathbf{F}_j^\bullet$ respectively. We define

eqnarray[eqnarray omitted — 276 chars of source]

Similar to $\mathbf{C}^\circ$ and $\mathbf{C}^\bullet$, we define $\mathbf{F}^\circ = (\mathbf{F}_1^\circ,\ldots, \mathbf{F}_L^\circ)$ and $ \mathbf{F}^\bullet = (\mathbf{F}_1^\bullet,\ldots, \mathbf{F}_N^\bullet ).$ Also, let

\[ \mathbf{b}_{\centerdot\centerdot} = (\mathbf{b}_{11},\ldots,\mathbf{b}_{1N},\ldots \ldots, \mathbf{b}_{LN})^\top \] be a generic matrix with dimensions matching $\pmb{\beta}_{\centerdot\centerdot} = (\pmb{\beta}_{11},\ldots, \pmb{\beta}_{1N},\ldots\ldots,\pmb{\beta}_{LN})^\top$.

Using these symbols, we introduce the following objective function:

eqnarray[eqnarray omitted — 311 chars of source]

where $\mathbf{C}_{ij}^\dag = \mathbf{M}_{(\mathbf{C}, \mathbf{C}_i^\circ)}+\mathbf{M}_{(\mathbf{C}, \mathbf{C}_j^\bullet)}$ and the following conditions hold:

eqnarray[eqnarray omitted — 394 chars of source]

remarkWe stress a few key points below. \begin{enumerate}[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • One may regard either $\mathbf{M}_{(\mathbf{C}, \mathbf{C}_i^\circ)}$ or $\mathbf{M}_{(\mathbf{C}, \mathbf{C}_j^\bullet)}$ as a projection matrix to be used to recover the global factors with slightly over-specified number of factors. The fact is reflected in the following operation: \begin{eqnarray*} \mathbf{F}_{ij}^\dag\mathbf{F}_{ij}^* = (\mathbf{M}_{(\mathbf{F},\mathbf{F}_i^\circ)}+\mathbf{M}_{(\mathbf{F},\mathbf{F}_j^\bullet)})(\mathbf{F},\mathbf{F}_i^\circ,\mathbf{F}_j^\bullet)=(\mathbf{0},\mathbf{M}_{(\mathbf{F},\mathbf{F}_j^\bullet)}\mathbf{F}_i^\circ,\mathbf{M}_{(\mathbf{F},\mathbf{F}_i^\circ)}\mathbf{F}_j^\bullet), \end{eqnarray*} where $\mathbf{F}_{ij}^\dag= \mathbf{M}_{(\mathbf{F},\mathbf{F}_i^\circ)}+\mathbf{M}_{(\mathbf{F},\mathbf{F}_j^\bullet)}$. By Assumption (ref).3, it is then easy to show that \begin{eqnarray*} \max_{i,j}\frac{1}{\sqrt{T}}\|(\mathbf{F}_{ij}^\dag\mathbf{F}_{ij}^*) - (\mathbf{0}, \mathbf{F}_i^\circ, \mathbf{F}_j^\bullet)\|=o_P(1), \end{eqnarray*} so a structure like $\mathbf{F}_{ij}^\dag$ only projects out the global factor, and keeps the local factors in the system as residuals. It reveals one may be able to recover global and local factors sequentially due to different signal strength. • The first three conditions of (ref) are standard in the literature. The last condition essentially requires $\mathbf{C} \perp \mathbf{C}^\circ $ and $\mathbf{C} \perp \mathbf{C}^\bullet $. Rather than labeling it as a constraint, it merely dictates the sequence of calculations when minimizing $Q(\mathbf{b}_{\centerdot\centerdot} ,\mathbf{C}_{\centerdot \centerdot})$. Consider $\mathbf{C}_i^\circ$ as an example. By (ref), we can write \[ \mathbf{C}_{ij}^\dag=\mathbf{M}_{(\mathbf{C}, \mathbf{C}_i^\circ)}+\mathbf{M}_{(\mathbf{C}, \mathbf{C}_j^\bullet)}=2\mathbf{I}_T-2\mathbf{P}_{\mathbf{C}}-\mathbf{P}_{\mathbf{C}_i^\circ}-\mathbf{P}_{\mathbf{C}_j^\bullet}, \] so $\mathbf{C}_{ij}^\dag$ maintains a structure similar to that in (ref). Minimizing $Q(\mathbf{b}_{\centerdot\centerdot} ,\mathbf{C}_{\centerdot \centerdot}) $ with respect to $\mathbf{C}_i^\circ$ is equivalent to maximizing the following quantity: \begin{eqnarray} &&\max_{\mathbf{C}_i^\circ} \operatorname*{\normalfontTr} \left(\sum_{j=1}^N(\mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot } \mathbf{b}_{ij})^\top \mathbf{C}_i^\circ\mathbf{C}_i^{\circ\top} (\mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot } \mathbf{b}_{ij}) \right)\nonumber \\ &=&\max_{\mathbf{C}_i^\circ}\operatorname*{\normalfontTr} \left(\mathbf{C}_i^{\circ\top} \mathbf{M}_{\mathbf{C}} \sum_{j=1}^N(\mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot } \mathbf{b}_{ij})(\mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot } \mathbf{b}_{ij})^\top \mathbf{M}_{\mathbf{C}} \mathbf{C}_i^{\circ}\right) , \end{eqnarray} where the equality holds because of $\mathbf{C}\perp \mathbf{C}_i^\circ$. Equation (ref) implies that in order to obtain the estimates of $\mathbf{F}_{i}^\circ$ and $\mathbf{F}_{j}^\bullet$, one should project out the estimate of $\mathbf{F}$ firstly, which is automatically taken into account in view of (ref) and (ref). Now, how to minimize global and local factors sequentially becomes even clearer. • Additionally, since \[ \lambda_{\min}(\mathbf{C}_{ij}^\dag)\ge \lambda_{\min} (\mathbf{M}_{(\mathbf{C}, \mathbf{C}_i^\circ)})+\lambda_{\min}(\mathbf{M}_{(\mathbf{C}, \mathbf{C}_j^\bullet)})= 0, \] we have $Q(\mathbf{b}_{\centerdot\centerdot} ,\mathbf{C}_{\centerdot \centerdot})\ge 0$ to ensure the non-negativity of the objective function. \end{enumerate}

Having explained the above points, we outline our estimation procedure below.

\hrule

itemize[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • Step 1 Perform the following optimization: \begin{eqnarray} (\widehat{\mathbf{b}}_{\centerdot\centerdot}, \widehat{\mathbf{C}}_{\centerdot \centerdot})=\operatorname*{\arg\!\min} Q(\mathbf{b}_{\centerdot \centerdot},\mathbf{C}_{\centerdot \centerdot}) , \end{eqnarray} where $\widehat{\mathbf{b}}_{\centerdot\centerdot} =(\widehat{\mathbf{b}}_{11},\ldots, \widehat{\mathbf{b}}_{LN})^\top$ and $\widehat{\mathbf{C}}_{\centerdot \centerdot} = (\widehat{\mathbf{C}},\widehat{\mathbf{C}}_1^\circ,\ldots, \widehat{\mathbf{C}}_L^\circ, \widehat{\mathbf{C}}_1^\bullet,\ldots, \widehat{\mathbf{C}}_N^\bullet)^\top$. Here, (ref) can be decomposed as follows: \begin{itemize}[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • $\widehat{\mathbf{b}}_{ij} = (\mathbf{X}_{ij\centerdot}^\top \widehat{\mathbf{C}}_{ij}^\dag \mathbf{X}_{ij\centerdot} )^{-1} \mathbf{X}_{ij\centerdot}^\top \widehat{\mathbf{C}}_{ij}^\dag \mathbf{Y}_{ij\centerdot},$ where $\widehat{\mathbf{C}}_{ij}^\dag= 2\mathbf{I}_T-2\mathbf{P}_{\widehat{\mathbf{C}}}-\mathbf{P}_{\widehat{\mathbf{C}}_i^\circ}-\mathbf{P}_{\widehat{\mathbf{C}}_j^\bullet}$; • $\widehat{\mathbf{C}}\widehat{\mathbf{V}} =\widehat{\pmb{\Sigma}} \widehat{\mathbf{C}}$, where $\widehat{\mathbf{V}}$ includes the largest $\ell$ eigenvalues of $\widehat{\pmb{\Sigma}} $, and \begin{eqnarray*} \widehat{\pmb{\Sigma}} = \frac{1}{LNT}\sum_{i=1}^L\sum_{j=1}^N(\mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot } \widehat{\mathbf{b}}_{ij} )(\mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot } \widehat{\mathbf{b}}_{ij})^\top; \end{eqnarray*} • For $\forall i\in [L]$, $\widehat{\mathbf{C}}_i^\circ\widehat{\mathbf{V}}_i^\circ =\widehat{\pmb{\Sigma}}_i^\circ \widehat{\mathbf{C}}_i^\circ$, where $\widehat{\mathbf{V}}_i^\circ$ includes the largest $\ell_i^\circ$ eigenvalues of $\widehat{\pmb{\Sigma}}_i^\circ$, and \begin{eqnarray*} \widehat{\pmb{\Sigma}}_i^\circ = \frac{1}{NT} \sum_{j=1}^N \mathbf{M}_{\widehat{\mathbf{C}}}(\mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot } \widehat{\mathbf{b}}_{ij})(\mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot } \widehat{\mathbf{b}}_{ij})^\top\mathbf{M}_{\widehat{\mathbf{C}}}; \end{eqnarray*} • For $\forall j\in [N]$, $\widehat{\mathbf{C}}_j^\bullet\widehat{\mathbf{V}}_j^\bullet =\widehat{\pmb{\Sigma}}_j^\bullet \widehat{\mathbf{C}}_j^\bullet$, where $\widehat{\mathbf{V}}_j^\bullet$ includes the largest $\ell_j^\bullet$ eigenvalues of $\widehat{\pmb{\Sigma}}_j^\bullet$, and \begin{eqnarray*} \widehat{\pmb{\Sigma}}_j^\bullet = \frac{1}{LT} \sum_{i=1}^L\mathbf{M}_{\widehat{\mathbf{C}}}(\mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot } \widehat{\mathbf{b}}_{ij})(\mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot } \widehat{\mathbf{b}}_{ij})^\top\mathbf{M}_{\widehat{\mathbf{C}}}. \end{eqnarray*} \end{itemize} • Step 2 Finally, update the estimate of $\pmb{\beta}_{ij}$ by \begin{eqnarray} \widetilde{\mathbf{b}}_{ij} = (\mathbf{X}_{ij\centerdot}^\top \mathbf{M}_{\widehat{\mathbf{C}}_{ij}} \mathbf{X}_{ij\centerdot} )^{-1} \mathbf{X}_{ij\centerdot}^\top \mathbf{M}_{\widehat{\mathbf{C}}_{ij}} \mathbf{Y}_{ij\centerdot}, \end{eqnarray} where $\widehat{\mathbf{C}}_{ij} =(\widehat{\mathbf{C}}, \widehat{\mathbf{C}}_i^\circ, \widehat{\mathbf{C}}_j^\bullet)$.

\hrule

The detailed numerical implementation in provided in Appendix (ref). In the next subsection, we propose the above estimation procedure, and establish the corresponding asymptotic properties.

Main Assumptions

In the following, we first provide some notation and mathematical symbols which will be repeatedly used in Assumptions and theoretical development of the paper. Key assumptions with their justifications will follow.

Notation:

Regressors: $\mathbf{X}_{i \centerdot \centerdot } =( \mathbf{X}_{i1\centerdot} ,\ldots, \mathbf{X}_{iN\centerdot} ) $ and $\mathbf{X}_{\centerdot j\centerdot } =( \mathbf{X}_{1j\centerdot} ,\ldots, \mathbf{X}_{Lj\centerdot} ) $;

Errors: $\pmb{\mathcal{E}}_{\centerdot \centerdot \centerdot}=(\pmb{\mathcal{E}}_{1\centerdot\centerdot} ,\ldots,\pmb{\mathcal{E}}_{L\centerdot\centerdot} )^\top$, $\pmb{\mathcal{E}}_{i\centerdot\centerdot} = ( \pmb{\mathcal{E}}_{i1\centerdot}, \ldots, \pmb{\mathcal{E}}_{iN\centerdot})$, $\pmb{\mathcal{E}}_{\centerdot j \centerdot} =( \pmb{\mathcal{E}}_{1j\centerdot}, \ldots, \pmb{\mathcal{E}}_{Lj\centerdot})$, $\mathbf{V}_{\centerdot\centerdot\centerdot}=(\mathbf{V}_{1\centerdot\centerdot},\cdots, \mathbf{V}_{L\centerdot\centerdot})^\top$, $\mathbf{V}_{i\centerdot\centerdot}=(\mathbf{V}_{i1\centerdot},\ldots,\mathbf{V}_{iN\centerdot})$, $\mathbf{V}_{\centerdot j\centerdot}=(\mathbf{V}_{1j\centerdot},\ldots,\mathbf{V}_{Lj\centerdot})$, and $\mathbf{V}_{ij\centerdot}=(\mathbf{v}_{ij1},\ldots, \mathbf{v}_{ijT})^\top$;

Loadings: $\pmb{\phi}_{ij}^* = (\pmb{\phi}_{ij}^\top, \pmb{\phi}_{ij}^{\circ \top}, \pmb{\phi}_{ij}^{\bullet \top})^\top$, $\pmb{\gamma}_{ij}^* = (\pmb{\gamma}_{ij}^\top, \pmb{\gamma}_{ij}^{\circ \top}, \pmb{\gamma}_{ij}^{\bullet \top})^\top$, $\pmb{\Gamma}_{i\centerdot}^\circ =(\pmb{\gamma}_{i1}^\circ,\ldots, \pmb{\gamma}_{iN}^\circ )^\top$, and $\pmb{\Gamma}_{\centerdot j}^\bullet = (\pmb{\gamma}_{1j}^\bullet,\ldots, \pmb{\gamma}_{Lj}^\bullet )^\top$.

assumption• \begin{itemize}[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • For the error components, we impose the following conditions: \begin{itemize}[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • $\|\pmb{\mathcal{E}}_{\centerdot \centerdot \centerdot}\|_2=O_P(\sqrt{LN}\vee \sqrt{T})$, and $\left(\frac{ \max_i\| \pmb{\mathcal{E}}_{i\centerdot\centerdot}\|_2}{\sqrt{NT}}+\frac{ \max_j\| \pmb{\mathcal{E}}_{\centerdot j\centerdot}\|_2}{\sqrt{LT}} \right)\cdot\log(L\vee N\vee T)=o_P(1)$; • $ |\frac{1}{LNT} \sum_{i,j}\pmb{\mathcal{E}}_{ij\centerdot}^\top \mathbf{F}_{ij}^* \pmb{\gamma}_{ij}^* | =O_P ( \frac{1}{\sqrt{LN \wedge T}} )$, and $\left| \frac{1}{LNT} \sum_{i,j}\pmb{\mathcal{E}}_{ij\centerdot}^\top \mathbf{X}_{ij\centerdot}\pmb{\beta}_{ij} \right| =O_P (\frac{1}{\sqrt{LN \wedge T}})$' • $\|\pmb{\mathbf{V}}_{\centerdot \centerdot \centerdot}\|_2=O_P(\sqrt{LN}\vee \sqrt{T})$, $\|\pmb{\mathbf{V}}_{i \centerdot \centerdot}\|_2=O_P(\sqrt{N}\vee \sqrt{T})$, and $\|\pmb{\mathbf{V}}_{\centerdot j\centerdot}\|_2=O_P(\sqrt{L}\vee \sqrt{T})$, for $\forall (i,j)$. \end{itemize} • For the loadings, we impose the following conditions: \begin{itemize}[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • $\frac{1}{LN}\sum_{i,j} (\pmb{\gamma}_{ij}^*\pmb{\gamma}_{ij}^{*\top}-E[\pmb{\gamma}_{ij}^*\pmb{\gamma}_{ij}^{*\top}] )=o_P(1)$ in which \begin{eqnarray*} 0<\min_{i,j} \lambda_{\min}(E[\pmb{\gamma}_{ij}^*\pmb{\gamma}_{ij}^{*\top}] )\le \max_{i,j} \lambda_{\max}(E[\pmb{\gamma}_{ij}^*\pmb{\gamma}_{ij}^{*\top}] )<\infty; \end{eqnarray*} • $\max_i \|\frac{1}{N} \pmb{\Gamma}_{i\centerdot}^{\circ \top}\pmb{\Gamma}_{i\centerdot}^\circ -\pmb{\Sigma}_{\pmb{\gamma}, i}^\circ \| =o_P (1) $ in which \begin{eqnarray*} 0<\min_i \lambda_{\min}(\pmb{\Sigma}_{\pmb{\gamma}, i}^\circ)\le \max_i \lambda_{\max}(\pmb{\Sigma}_{\pmb{\gamma}, i}^\circ)<\infty; \end{eqnarray*} • $\max_j \|\frac{1}{L} \pmb{\Gamma}_{\centerdot j}^{\bullet \top}\pmb{\Gamma}_{\centerdot j}^\bullet -\pmb{\Sigma}_{\pmb{\gamma}, j}^\bullet \| =o_P (1)$ in which \begin{eqnarray*} 0<\min_j \lambda_{\min}(\pmb{\Sigma}_{\pmb{\gamma}, j}^\bullet)\le \max_j \lambda_{\max}(\pmb{\Sigma}_{\pmb{\gamma}, j}^\bullet)<\infty. \end{eqnarray*} \end{itemize} • For the factors, we impose the following conditions: \begin{itemize}[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • $\log(L\vee N\vee T)\cdot\max_{i,j}\| \frac{1}{T} \mathbf{F}_{ij}^{*\top} \mathbf{F}_{ij}^* -\mathbf{I}_{\ell+\ell_i^\circ+\ell_j^\bullet}\| =o_P (1)$; • $\|\mathbf{F}^\circ\|_2 =O_P(\sqrt{L}\vee \sqrt{T})$ and $\|\mathbf{F}^\bullet\|_2 =O_P(\sqrt{N}\vee \sqrt{T})$. \end{itemize} • For the regressors, we impose the following conditions: \begin{itemize}[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • $\max_i\frac{1}{NT} \| \mathbf{X}_{i\centerdot\centerdot}\|^2 =O_P(1)$; • $\max_j\frac{1}{LT} \| \mathbf{X}_{\centerdot j\centerdot}\|^2 =O_P(1)$. \end{itemize} • Suppose that the following conditions are satisfied with probability approaching 1: \[\inf_{\mathbf{C}_{\centerdot\centerdot}} \lambda_{\min}\left(\frac{1}{T}\mathbf{D}_{\centerdot \centerdot}\right)>0, \] where $\mathbf{D}_{\centerdot \centerdot} = \operatorname*{\normalfont\textrm{diag}}\{\mathbf{D}_{11},\ldots, \mathbf{D}_{LN}\}$, $\mathbf{D}_{ij, 1} = \mathbf{X}_{ij\centerdot}^\top\mathbf{C}_{ij}^\dag \mathbf{X}_{ij\centerdot} $, $\mathbf{D}_{ij, 2} = \pmb{\gamma}_{ij} \otimes (\mathbf{C}_{ij}^\dag \mathbf{X}_{ij\centerdot})$, and $\mathbf{D}_{ij} =\mathbf{D}_{ij, 1} -\mathbf{D}_{ij, 2}^\top (E[\pmb{\gamma}_{ij} \pmb{\gamma}_{ij}^{\top}] \otimes \mathbf{I}_T)^{-1}\mathbf{D}_{ij, 2}$. \end{itemize}

Assumption (ref) imposes a set of regular conditions. All of the conditions involving $\max_i$, $\max_j$ and $\max_{ij}$ can be easily justified in the same fashion as in Lemma (ref) in Appendix A below. However, as we cannot define some type of mixing conditions along with the dimension $i$ or $j$, so we state them as they are via a set of high level assumptions. Assumption (ref).1.(a) specifies the spectral norm orders of random errors, which can be satisfied by various examples discussed in the supplementary Appendix S.2 of Moon2017. In particular, it can be justified by i.i.d. $\varepsilon_{ijt}$ with zero mean and uniformly bounded fourth moment $E[\varepsilon_{ijt}^4]$. Following the arguments in Lemma S.2.1 of Moon2017, we can demonstrate that Assumption (ref).1.(a) also holds for the MA($\infty$) process: $\varepsilon_{ijt} = \sum_{\tau=0}^\infty a_{ij\tau}e_{ij,t-\tau}$, where $e_{ijt}$ are i.i.d. with zero mean and uniformly bounded fourth moment, and the MA coefficients ${a_{ij\tau}}$ satisfy certain uniform boundedness conditions.

Assumption (ref).1.(b) is met by strictly exogenous and i.i.d. or $\alpha$-mixing errors. Assumption (ref).2 regulates the behaviour of factor loadings, which can be justified by identically distributed factor loadings with weak cross--sectional dependence. Assumption (ref).3 imposes conditions on the factors, which can be easily satisfied. For instance, Assumption (ref).3.(a) can be met by $\mathbf{f}_t$, $\mathbf{f}_{it}^\circ$, and $\mathbf{f}_{jt}^\bullet$ that are independently generated from a multivariate standard normal distribution. In Assumption (ref).3.(b), the required orders of spectral norm of $\mathbf{F}^\circ$ and $\mathbf{F}^\bullet$ can be fulfilled by the i.i.d. or MA($\infty$) processes, similar to those discussed for Assumption (ref).1.

As explained in Remark (ref), Assumption (ref).4 actually does not impose too many conditions on variables used to generate $\mathbf{x}_{ijt}$. For instance, one may have weak factor signal through the factor loadings associated with $\mathbf{x}_{ijt}$. If that is the case, one may add certain factor structure to $\mathbf{x}_{ijt}$. In short, we focus on the current estimation approach in Section (ref). We comment on Assumption (ref).5 as follows. This assumption only includes the factor loadings of the global factors, so it is similar to Assumption A of Bai, but with a lightly over specified number of factors that kicks in via $\mathbf{C}_{ij}^\dag= \mathbf{M}_{(\mathbf{C}, \mathbf{C}_i^\circ)}+\mathbf{M}_{(\mathbf{C}, \mathbf{C}_j^\bullet)}$ by the construction. This is indeed achievable as argued in Moon. As we have not imposed too many conditions on the components of $\mathbf{x}_{ijt}$ up to this point, we state Assumption (ref).5 as it is. In a very extreme case in which $\mathbf{x}_{ijt}\equiv \mathbf{v}_{ijt}$, this condition reduces to $\min_{i,j}\lambda_{\min}(\mathbf{V}_{ij\centerdot}^\top\mathbf{V}_{ij\centerdot}/T)$ which can be realized under a set of additional conditions on $\mathbf{v}_{ijt}$ and in connection with Lemma (ref) in Appendix A below.

assumption• \begin{itemize}[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • Suppose that $\{ \pmb{\zeta}_{ijt} : = ( \varepsilon_{ijt}, \mathbf{v}_{ijt}^\top)^\top\}$ are independent of the other variables, and $E[\pmb{\zeta}_{ijt} ]=\mathbf{0}$ and $ E[\varepsilon_{ijt} | \mathbf{v}_{ijt}]=0$, and for $\forall (i,j)$, $ \{\pmb{\zeta}_{ijt} | t\in [T]\}$ are strictly $\alpha$-mixing process with the $\alpha$-mixing coefficient \begin{eqnarray*} \alpha_{ij}(t) = \sup_{A\in \mathcal{F}_{ij,-\infty}^0, B\in \mathcal{F}_{ij,t}^\infty} |P(A)P(B) -P(AB)| \end{eqnarray*} satisfying $ \max_{i,j}\sum_{t=1}^\infty [\alpha_{ij}(t)]^{\kappa_{ij}/(4+\kappa_{ij})} < \infty$ with $\kappa_{ij}>0$ and $\max_{i,j} E\|\pmb{\zeta}_{ij1}\|^{2+\kappa_{ij}} <\infty$, where $\mathcal{F}_{ij,-\infty}^0$ and $\mathcal{F}_{ij,t}^\infty$ are the $\sigma$-algebras generated by $\{\pmb{\zeta}_{ijs}: s \leq 0\}$ and $\pmb{\zeta}_{ijs}: s \geq t\}$, respectively. Finally, there exits $\kappa^*\in (0, \min_{i,j}\kappa_{ij})$ such that $LN/T^{1+\kappa^*/4}\to 0$. • $\max_{i,j} (\|\pmb{\gamma}_{ij}^\circ\|^2+\|\pmb{\gamma}_{ij}^\bullet\|^2)=O_P(\log(LN))$. • The minimum eigenvalue of $\pmb{\Sigma}_{\mathbf{v},ij}=\lim_{T\to\infty}\frac{1}{T}E[\mathbf{V}_{ij\centerdot}^\top \mathbf{V}_{ij\centerdot}]$ is positive and uniformly bounded away from zero for all $i$ and $j$. \end{itemize}

The mixing condition of Assumption (ref).1 is necessary to derive Lemma (ref), which is then further required to derive Lemma (ref) and the asymptotic distribution. The requirements about $\kappa_{ij}$'s impose conditions on the existence of some moments of $\pmb{\zeta}_{ijt}$, which can be easily fulfilled, e.g., $\pmb{\zeta}_{ijt}$ follows a normal distribution. Assumption (ref).2 is standard.

Asymptotic Properties

Under the main assumptions, we are now ready to establish the consistency and asymptotic normality results.

lemmaFor {\normalfont Step 1}, under Assumption (ref), as $(L,N,T)\to (\infty,\infty,\infty)$, we have \begin{itemize}[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • $\frac{1}{\sqrt{LN}}\| \pmb{\beta}_{\centerdot \centerdot} -\widehat{\mathbf{b}}_{\centerdot \centerdot} \|=o_P(1)$; • $\|\mathbf{P}_{\mathbf{F}} -\mathbf{P}_{\widehat{\mathbf{C}}} \|=O_P\left( \frac{1}{\sqrt{LN}}\|\pmb{\beta}_{\centerdot \centerdot}-\widehat{\mathbf{b}}_{\centerdot\centerdot}\| +\frac{1}{\sqrt{L\wedge N\wedge T}}\right)$. \end{itemize}

The establishment of Lemma (ref) does not require an excessive number of conditions for the variables involved in (ref). It shows that we can recover $\pmb{\beta}_{\centerdot \centerdot}$ as a whole. As expected, the global factor $\mathbf{F}$ can be estimated firstly by putting everything together.

As shown in (ref) and (ref), in order to recover $\mathbf{F}_i^\circ$ and $\mathbf{F}_j^\bullet$, we need to divide the data into different blocks as presented in the estimation procedure, so $\widehat{\mathbf{b}}_{ij}$'s are automatically partitioned into blocks accordingly. As a consequence, we can no longer use the overall consistency of $\widehat{\mathbf{b}}_{\centerdot \centerdot}$ derived in Lemma (ref). To proceed, we further impose Assumption (ref), and summarize the results associated with $\mathbf{F}_i^\circ$ and $\mathbf{F}_j^\bullet$ in Lemma (ref) below.

lemmaFor {\normalfont Step 1}, under Assumptions (ref) and (ref), as $(L,N,T)\to (\infty,\infty,\infty)$, \begin{itemize}[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • $\max_{i, j}\| \pmb{\beta}_{ij}-\widehat{\mathbf{b}}_{ij}\|=o_P(1)$; • $\max_i\|\mathbf{P}_{ \mathbf{F}_i^\circ} -\mathbf{P}_{ \widehat{\mathbf{C}}_i^\circ}\|=O_P (\max_{i,j} \|\pmb{\beta}_{ij}-\widetilde{\mathbf{b}}_{ij}\| + \frac{\sqrt{\log(LN)}}{\sqrt{N \wedge T}} +\frac{1}{\sqrt{L}}+ \frac{\max_i\|\pmb{\mathcal{E}}_{i\centerdot\centerdot}\|_2}{\sqrt{N T}} +\frac{\max_i \| \mathbf{F}_i^{\circ\top} \mathbf{F}\|_2 }{T} )$; • $\max_j\|\mathbf{P}_{ \mathbf{F}_j^\bullet} -\mathbf{P}_{ \widehat{\mathbf{C}}_j^\bullet}\|=O_P (\max_{i,j} \|\pmb{\beta}_{ij}-\widehat{\mathbf{b}}_{ij}\| + \frac{\sqrt{\log(LN)}}{\sqrt{L \wedge T}} +\frac{1}{\sqrt{N}} + \frac{\max_j\|\pmb{\mathcal{E}}_{\centerdot j\centerdot}\|_2}{\sqrt{L T}}+\frac{\max_j\| \mathbf{F}_j^{\bullet\top} \mathbf{F}\|_2}{T} )$. \end{itemize}

Lemma (ref) justifies the validity of our approach across all blocks involved in Step 1. With this, we have successfully recovered all key components of model (ref), thus concluding our examination of Step 1 in the estimation procedure.

Now, we proceed to Step 2, where we update the estimator for $\pmb{\beta}_{ij}$. Additionally, we consider two mean group estimators to infer the average values of $\pmb{\beta}_{ij}$ on two cross-sectional dimensions, respectively. For $\forall (i,j)$, define

equation*[equation* omitted — 204 chars of source]

To facilitate the theoretical development, we introduce some additional notation. Let $\mathcal{R}_{ij}^{*}=(\mathcal{R}_{F}, \mathcal{R}_{F,i}^\circ, \mathcal{R}_{F,j}^\bullet)$, where $\mathcal{R}_F = \sum_{ij}(\mathbf{F}_{i}^{\circ}\pmb{\gamma}_{ij}^\circ+\mathbf{F}_{j}^{\bullet}\pmb{\gamma}_{ij}^\bullet)\pmb{\gamma}_{ij}^\top (\pmb{\Gamma}_{\centerdot \centerdot }^\top \pmb{\Gamma}_{\centerdot \centerdot })^{-1}$, $\mathcal{R}_{F,i}^\circ=\sum_{j=1}^N(\mathbf{F}\pmb{\gamma}_{ij}+\mathbf{F}_{j}^{\bullet}\pmb{\gamma}_{ij}^\bullet)\pmb{\gamma}_{ij}^{\circ\top} (\pmb{\Gamma}_{i\centerdot }^{\circ\top} \pmb{\Gamma}_{i\centerdot }^\circ)^{-1}$ and $\mathcal{R}_{F,j}^\bullet = \sum_{i=1}^L(\mathbf{F}\pmb{\gamma}_{ij}+\mathbf{F}_{i}^{\circ}\pmb{\gamma}_{ij}^\circ)\pmb{\gamma}_{ij}^{\bullet\top} (\pmb{\Gamma}_{\centerdot j }^{\bullet\top} \pmb{\Gamma}_{\centerdot j }^\bullet)^{-1}$.

theoremLet Assumptions (ref) and (ref) hold. As $(L,N,T)\to (\infty,\infty,\infty)$, \begin{itemize} • for $\forall (i,j)$, if $\sqrt{T}/(L\wedge N) \to 0$, \begin{eqnarray*} \sqrt{T} (\widetilde{\mathbf{b}}_{ij} - \pmb{\beta}_{ij} ) \to_D N(\mathbf{0}, \pmb{\Sigma}_{\mathbf{v},ij}^{-1}\pmb{\Sigma}_{\mathbf{v},\varepsilon,ij}\pmb{\Sigma}_{\mathbf{v},ij}^{-1}), \end{eqnarray*} where $\pmb{\Sigma}_{\mathbf{v},\varepsilon,ij} =\lim_{T\to\infty}\frac{1}{T}E[\mathbf{V}_{ij\centerdot}^\top\pmb{\mathcal{E}}_{ij\centerdot } \pmb{\mathcal{E}}_{ij\centerdot }^\top \mathbf{V}_{ij\centerdot}]$. • for $\forall i$, if $\sqrt{T}/(L\wedge N) \to 0$ and $\sqrt{N}/(L\wedge T) \to 0$, \begin{eqnarray*} \sqrt{NT} (\widetilde{\mathbf{b}}_{i\centerdot} - \pmb{\beta}_{i\centerdot} -\mathcal{B}_{i\centerdot} ) \to_D N(\mathbf{0}, \pmb{\Sigma}_{\mathbf{b},i\centerdot}), \end{eqnarray*} where $\mathcal{B}_{i\centerdot}=\frac{1}{NT}\sum_{j=1}^N\pmb{\Sigma}_{\mathbf{v},ij}^{-1} \pmb{\phi}_{ij}^{*\top} \mathcal{R}_{ij}^{*\top} \mathbf{M}_{\widehat{\mathbf{C}}_{ij}} \mathcal{R}_{ij}^* \pmb{\gamma}_{ij}^*$ is the bias term with the probability order of $O_P(\frac{1}{L\wedge N \wedge T})$, and $\pmb{\Sigma}_{\mathbf{b},i\centerdot}=\lim_{N,T\to\infty} \frac{1}{NT}\sum_{j_1,j_2=1}^N \pmb{\Sigma}_{\mathbf{v},ij_1}^{-1} E[\mathbf{V}_{ij_1\centerdot}^\top\pmb{\mathcal{E}}_{ij_1\centerdot } \pmb{\mathcal{E}}_{ij_2\centerdot }^\top \mathbf{V}_{ij_2\centerdot}]\pmb{\Sigma}_{\mathbf{v},ij_2}^{-1}$. • for $\forall j$, if $\sqrt{T}/(L\wedge N) \to 0$ and $\sqrt{L}/(N\wedge T) \to 0$, \begin{eqnarray*} \sqrt{LT} (\widetilde{\mathbf{b}}_{\centerdot j} - \pmb{\beta}_{\centerdot j}-\mathcal{B}_{\centerdot j} ) \to_D N(\mathbf{0}, \pmb{\Sigma}_{\mathbf{b},\centerdot j}), \end{eqnarray*} where $\mathcal{B}_{\centerdot j}=\frac{1}{LT}\sum_{i=1}^L\pmb{\Sigma}_{\mathbf{v},ij}^{-1} \pmb{\phi}_{ij}^{*\top} \mathcal{R}_{ij}^{*\top} \mathbf{M}_{\widehat{\mathbf{C}}_{ij}} \mathcal{R}_{ij}^* \pmb{\gamma}_{ij}^*$ is the bias term with the probability order of $O_P(\frac{1}{L\wedge N \wedge T})$, and $\pmb{\Sigma}_{\mathbf{b},\centerdot j}=\lim_{L,T\to\infty} \frac{1}{LT}\sum_{i_1,i_2=1}^L \pmb{\Sigma}_{\mathbf{v},i_1j}^{-1} E[\mathbf{V}_{i_1j\centerdot}^\top\pmb{\mathcal{E}}_{i_2j\centerdot } \pmb{\mathcal{E}}_{i_2j\centerdot }^\top \mathbf{V}_{i_2j\centerdot}]\pmb{\Sigma}_{\mathbf{v},i_2j}^{-1}$. \end{itemize}

Theorem (ref).(1) provides the asymptotic distribution for each $\widetilde{\mathbf{b}}_{ij}$. Meanwhile, Theorem (ref).(2-3) establish the asymptotic normality for the mean group estimators, which incorporate bias terms to achieve the optimal rates of convergence ($\sqrt{NT}$ and $\sqrt{LT}$, respectively). The presence of bias terms arises due to the hierarchical factor structure, which can be addressed using conventional methods such as the analytical approach or the Jackknife procedure Chen2021. In Appendix (ref) below, we propose a Jackknife method to remove the bias terms present in the theorem.

To conduct practical inference, it is essential to deal with the time series correlation involved in the term $\pmb{\Sigma}_{\mathbf{v},\varepsilon,ij}$. Therefore, we now focus on Theorem (ref).(1) and then develop a new dependent wild bootstrap procedure by considerably extending GPY2023 for the standard time--series and cross-sectional panel data setting to the following setting where the double--index $(i,j)$ involved can be both cross--sectional. For the second and third parts of Theorem (ref), the corresponding bootstrap procedure involves a two--step approach in each case. The first step is a bias correction procedure as discussed in Appendix (ref) below, and the second step develops a bootstrap procedure for each case in almost the same way as proposed below.

Bootstrap Procedure: \hrule

itemize[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • Let $\pmb{\xi}=(\xi_{1},\ldots, \xi_{T})^{\top}$ be an $m$-dependent time series for each bootstrap draw, and let the components of $\pmb{\xi}$ satisfy that \begin{eqnarray*} E[\xi_{t}]=0,\quad E|\xi_{t}|^{2} =1,\quad E|\xi_{t}|^{2+\kappa_{ij}/2 }<\infty,\quad E[\xi_{t}\xi_{s}]=a\left( \frac{t-s}{m} \right) , \end{eqnarray*} where $m\to\infty$, $\kappa_{ij}$ is a constant and defined in Assumption (ref), and $a(\cdot)$ is a symmetric kernel defined on $[-1,1]$ satisfying that $a(0)=1$ and $K_{a}(x)=\int_{\mathbb{R}} a(u)e^{-iux}du\ge0$ for $x\in\mathbb{R}$. • For $\forall (i,j)$, construct a new set of dependent variables by $\mathbf{Y}_{ij\centerdot}^{*} = \mathbf{X}_{ij\centerdot} \widehat{\mathbf{b}}_{ij}+ \widehat{\pmb{\mathcal{E}}}_{ij\centerdot}\circ\pmb{\xi}$, where $\widehat{\pmb{\mathcal{E}}}_{ij\centerdot} = \mathbf{Y}_{ij\centerdot} -\mathbf{X}_{ij\centerdot} \widehat{\mathbf{b}}_{ij}$. Accordingly, the new estimate of $\pmb{\beta}_{ij}$ is obtained using \begin{eqnarray*} \widetilde{\mathbf{b}}_{ij}^*=\left(\mathbf{X}_{ij\centerdot}^\top \mathbf{M}_{\widehat{\mathbf{C}}_{ij}} \mathbf{X}_{ij\centerdot}\right)^{-1} \mathbf{X}_{ij\centerdot}^\top \mathbf{M}_{\widehat{\mathbf{C}}_{ij}} \mathbf{Y}_{ij\centerdot}^*, \end{eqnarray*} • We repeat the above procedure $\mathcal{L}$ times.

\hrule

theoremLet Assumptions (ref) and (ref) hold. Assume further that $\frac{m}{\sqrt{Th}} \to 0$ and $mT^{2/\underline{\kappa}_{ij}-1}\rightarrow 0$, where $\underline{\kappa}_{ij}=\min(\kappa^\ast_{ij},4)$, $\kappa^\ast_{ij}=1/ (\frac{1}{2+\kappa_{ij}}+\frac{\kappa_{ij}}{2(4+\kappa_{ij})} )$, and $\kappa_{ij}$ is defined in Assumption (ref). As $(L,N,T)\to (\infty,\infty,\infty)$, for $\forall (i,j)$ \begin{eqnarray*} \sup_{w\in \mathbb{R}}\left| \normalfont Pr^* (\sqrt{T}(\widetilde{\mathbf{b}}_{ij}^* -\widetilde{\mathbf{b}}_{ij})\le w) -\Pr (\sqrt{T}(\widetilde{\mathbf{b}}_{ij}-\pmb{\beta}_{ij})\le w)\right| =o_P(1). \end{eqnarray*}

Theorem (ref) provides a procedure to recover the asymptotic distribution for each given $(i,j)$. It is worth mentioning that the DWB method is initially introduced in shao2010 in the context of time series analysis, wherein a comprehensive comparison between the DWB and some existing bootstrap methods can be found.

Selection of the Numbers of Factors

To conclude our theoretical investigation, we provide a procedure to select the numbers of factors. For notational simplicity, we define $\pmb{\ell}^\circ =(\ell_1^\circ,\ldots, \ell_L^\circ)^\top$ and $\pmb{\ell}^\bullet =(\ell_1^\bullet,\ldots, \ell_N^\bullet)^\top.$ Thus, our goal is to estimate $\ell$, $\pmb{\ell}^\circ$, and $\pmb{\ell}^\bullet$.

First, we emphasize that we can still get consistent estimation of $\widehat{\mathbf{b}}_{\centerdot\centerdot}$ by prescribing a $d_{\max}$ to be the numbers of global, country and industry factors. As long as $d_{\max}\ge (\ell \vee (\max_i \ell_i^\circ) \vee (\max_i\ell_j^\bullet))$, the estimation procedure of Section (ref), and the results developed in Section (ref) still hold true. In connection with the estimation procedure of Section (ref), we provide the following steps to estimate $\ell$ firstly, and then estimate the elements of $\pmb{\ell}^\circ$ and $\pmb{\ell}^\bullet$ sequentially.

\hrule

itemize[leftmargin=*, itemsep=0.5pt, parsep=0.5pt, topsep=0.6pt] • Step 1* Prescribe a $d_{\max}$ to be the numbers of global, country and industry factors, and implement the estimation procedure of Section (ref) to get $\widetilde{\mathbf{b}}_{\centerdot \centerdot} =(\widetilde{\mathbf{b}}_{11},\ldots, \widetilde{\mathbf{b}}_{1N},\ldots, \widetilde{\mathbf{b}}_{LN})^\top$. • Step 2* Using $\widetilde{\mathbf{b}}_{\centerdot \centerdot}$, we calculate $\widehat{\pmb{\Sigma}}$ as in Section (ref), and estimate the number of global factors as follows: \begin{eqnarray} \widehat{\ell} =\operatorname*{\arg\!\min}_{0\le s \le d_{\max} } \Big\{\frac{\widehat{\lambda}_{s+1}}{\widehat{\lambda}_s }\cdot I(\widehat{\lambda}_s \ge \omega)+I(\widehat{\lambda}_s < \omega) \Big\}, \end{eqnarray} where $\omega =[\log (L\vee N\vee T)]^{-1}$, $\widehat{\lambda}_0 =1$ is a mock eigenvalue, and $\widehat{\lambda}_s$ stands for the $s^{th}$ largest eigenvalue of $\widehat{\pmb{\Sigma}}$. Record the eigenvectors corresponding to the largest $\widehat{\ell}$ eigenvalues in $\widehat{\mathbf{C}}$. • Step 3* Having obtained $\widetilde{\mathbf{b}}_{\centerdot \centerdot}$ and $\widehat{\mathbf{C}}$, we sequentially estimate $\ell_i^\circ$'s and $\ell_j^\bullet$'s as follows: \begin{eqnarray} \widehat{\ell}_i^\circ &=&\operatorname*{\arg\!\min}_{0\le s \le d_{\max} } \Big\{ \frac{\widehat{\lambda}_{i,s+1}^\circ}{\widehat{\lambda}_{i,s}^\circ} \cdot I(\widehat{\lambda}_{i, s}^\circ \ge \omega)+ I(\widehat{\lambda}_{i, s}^\circ < \omega)\Big\},\nonumber \\ \widehat{\ell}_j^\bullet &=&\operatorname*{\arg\!\min}_{0\le s \le d_{\max} } \Big\{ \frac{\widehat{\lambda}_{j,s+1}^\bullet}{\widehat{\lambda}_{j,s}^\bullet} \cdot I(\widehat{\lambda}_{j, s}^\bullet \ge \omega)+ I(\widehat{\lambda}_{j, s}^\bullet< \omega)\Big\}, \end{eqnarray} where $\widehat{\lambda}_{i,0}^\circ=\widehat{\lambda}_{j,0}^\bullet=1$ are mock eigenvalues, and $\widehat{\lambda}_{i,s}^\circ$ and $\widehat{\lambda}_{j,s}^\bullet$ stand for the $s$-th largest eigenvalues of $\widehat{\pmb{\Sigma}}_i^\circ $ and $\widehat{\pmb{\Sigma}}_j^\bullet$ respectively.

\hrule

The estimators in (ref) and (ref) can be considered as extensions of LY12. However, as pointed out by LY12, it remains unresolved to bound the ratio associated with the eigenvalues which converge to 0 from below. To bypass this unresolved issue, we introduce a data driven tuning parameter $\omega$. The idea is that although it is challenging to study a ratio with a denominator converging to 0, we can discard this ratio and construct a V-shape curve by employing the indicator function in (ref) and (ref), respectively.

theoremUnder Assumptions (ref) and (ref), as $(L,N,T)\to (\infty,\infty,\infty)$, $\Pr(\widehat{\ell} =\ell, \widehat{\pmb{\ell}}^\circ= \pmb{\ell}^\circ, \widehat{\pmb{\ell}}^\bullet= \pmb{\ell}^\bullet)\to 1,$ where $\widehat{\pmb{\ell}}^\circ =(\widehat{\ell}_1^\circ,\ldots, \widehat{\ell}_L^\circ)^\top $ and $\widehat{\pmb{\ell}}^\bullet=(\widehat{\ell}_1^\bullet,\ldots, \widehat{\ell}_L^\bullet)^\top $.

Up to this point, we have successfully recovered all the unknown quantities as mentioned in the beginning of Section (ref). In the next section, we examine the above theoretical results using extensive simulations.

remarkTo conclude this section, we note that practically it is not guaranteed that each country (or industry) has the same number of individuals, so it brings certain unbalanceness to the data set. Formally, $\forall i\in [L]$, we may have $j\in [N_i]$, which is exactly the same as our empirical study of Section (ref). After carefully checking the proofs, we can see that as long as $\min_i N_i\to \infty$ and $\frac{T}{\min_i N_i^2}\to 0$, all the aforementioned results are still valid with certain modifications.

Simulation Studies

\setcounter{equation}{0}

In this section, we conduct extensive numerical studies to examine the theoretical findings. Specifically, we provide simulation studies, and the main data generating process (DGP) follows (ref). In what follows, we provide details about each component involved.

To capture heterogeneity of the coefficients, we generate them by $\pmb{\beta}_{ij} =(0.5+i/L, 0.5+j/N)^\top$ for $(i,j)\in [L]\times[N]$, which implies $d=2$. For factors and loadings, we first generate the numbers of factors. Specifically, we let $\ell =2$, and let $\ell_i^\circ$ and $ \ell_j^\bullet$ be randomly selected from the set $\{ 0,1,2\}$ with equal probability attached to each possible outcome. The design is to show that we allow $\ell_i^\circ$ and $ \ell_j^\bullet$ to be 0, which is meaningful for practical analysis. Therefore, for each generated dataset, only $\ell$ is fixed. The values of $\pmb{\ell}^\circ$ and $\pmb{\ell}^\bullet$ vary at each simulation replication.

Sequentially, we generate factors and loadings, and let them be independently drawn from some normal distributions: $\mathbf{f}_t\sim N(\mathbf{0},\mathbf{I}_\ell)$, $\mathbf{f}_{it}^\circ\sim N(\mathbf{0},2\cdot\mathbf{I}_{\ell_i^\circ})$ if $\ell_i^\circ>0$, $\mathbf{f}_{jt}^\bullet\sim N(\mathbf{0}, 2\cdot\mathbf{I}_{\ell_j^\bullet})$ if $\ell_j^\bullet>0$, $\pmb{\gamma}_{ij}\sim N(\mathbf{1},\mathbf{I}_\ell)$, $\pmb{\gamma}_{ij}^\circ\sim N(\mathbf{0},\mathbf{I}_{\ell_i^\circ})$ if $\ell_i^\circ>0$, $\pmb{\gamma}_{ij}^\bullet\sim N(-\mathbf{1},\mathbf{I}_{\ell_j^\bullet})$ if $\ell_j^\bullet>0$, $\pmb{\phi}_{ij, s}\sim N(\mathbf{1},\mathbf{I}_\ell)$, $\pmb{\phi}_{ij, s}^\circ\sim N(\mathbf{0},\mathbf{I}_{\ell_i^\circ})$ if $\ell_i^\circ>0$, $\pmb{\phi}_{ij, s}^\bullet\sim N(-\mathbf{1},\mathbf{I}_{\ell_j^\bullet})$ if $\ell_j^\bullet>0$, where $s\in [d]$, $\pmb{\phi}_{ij, s}$, $\pmb{\phi}_{ij, s}^\circ$, and $\pmb{\phi}_{ij, s}^\bullet$ stand for the $s^{th}$ columns of $\pmb{\phi}_{ij}$, $\pmb{\phi}_{ij}^\circ$, and $\pmb{\phi}_{ij}^\bullet$ respectively. We note that normal distributions are not really necessary here, as the choice of some other distributions doesn't affect the finite-sample performance. To present the main steps for the simulation and save the space, we do not further explore alternative options.

We generate residuals as follows:

eqnarray*[eqnarray* omitted — 367 chars of source]

where $\rho_\varepsilon=\rho_v =0.1$, $\pmb{\Sigma}_\varepsilon=\pmb{\Sigma}_v = \{ 0.2^{\|(i_1,j_1)-(i_2,j_2)\|}\}_{LN\times LN}$ with $(i,j)\in [L]\times[N]$, and

eqnarray*[eqnarray* omitted — 289 chars of source]

in which $\mathbf{v}_{ijt,s}$ is the $s^{th}$ element of $\mathbf{v}_{ijt}$, and $s\in [d]$. Based on the above DGP, we have weak serial correlation and cross-sectional dependence among $(\varepsilon_{ijt}, \mathbf{v}_{ijt}^\top)^\top$.

Finally, we assemble $y_{ijt}$'s and $\mathbf{x}_{ijt}$'s as in (ref). For each dataset, we carefully calculate all the unknown quantities following the numerical implementation of Appendix (ref).

To evaluate the performance of our proposed method, we repeat the above procedure $R=1000$ times, and define several criteria. For notational simplicity, the subindex $_k$ always stands for the corresponding quantities obtained at the $k^{th}$ replication in what follows. That said, we use root mean square errors (RMSE) to measure our estimates about $\pmb{\beta}_{\centerdot\centerdot}$, $\mathbf{F}$, $\mathbf{F}^\circ$, and $\mathbf{F}^\bullet$:

eqnarray[eqnarray omitted — 782 chars of source]

where we let $\mathbf{P}_{\widehat{\mathbf{C}}_k} =\mathbf{0}$ if $\widehat{\ell}=0$, and similarly, we let $\mathbf{P}_{\widehat{\mathbf{C}}_{i,k}^\circ}=\mathbf{0}$, $\mathbf{P}_{\mathbf{F}_{i,k}^\circ}=\mathbf{0}$, $\mathbf{P}_{\widehat{\mathbf{C}}_{j,k}^\bullet}=\mathbf{0}$, and $\mathbf{P}_{\mathbf{F}_{j,k}^\bullet}=\mathbf{0}$ if $\widehat{\ell}_{i,k}^\circ =0$, $\ell_{i,k}^\circ =0$, $\widehat{\ell}_{j,k}^\bullet =0$, and $\ell_{j,k}^\bullet =0$ respectively. For the purpose of comparison, we estimate $\pmb{\beta}_{\centerdot\centerdot}$, $\mathbf{F}$, $\mathbf{F}^\circ$ and $\mathbf{F}^\bullet$ by assuming that $\ell$, $\pmb{\ell}^\circ$, and $\pmb{\ell}^\bullet$ are known, and calculate $\text{RMSE}_{\pmb{\beta}}^*$, $\text{RMSE}_{\mathbf{F}}^*$, $\text{RMSE}_{\mathbf{F}^{\circ}}^*$, and $\text{RMSE}_{\mathbf{F}^{\bullet}}^*$ as in (ref).

To evaluate the estimation about the numbers of factors, we define the following rates:

eqnarray[eqnarray omitted — 1,044 chars of source]

Apparently, $\text{Rate}_\ell$, $\text{Rate}_{\pmb{\ell}^\circ}$, and $\text{Rate}_{\pmb{\ell}^\bullet}$ measure the percentages of correctly identifying the numbers of global, country and industry factors; $\text{Rate}_{\ell}^-$, $\text{Rate}_{\pmb{\ell}^\circ}^-$, and $\text{Rate}_{\pmb{\ell}^\bullet}^-$ measure the percentages of under estimation; $\text{Rate}_{\ell}^+$, $\text{Rate}_{\pmb{\ell}^\circ}^+$, and $\text{Rate}_{\pmb{\ell}^\bullet}^+$ measure the percentages of over estimation.

Finally, we evaluate the coverage rates (CR) of the bootstrap procedure. For $s\in [d]$,

eqnarray[eqnarray omitted — 168 chars of source]

where $\pmb{\beta}_{ij,s}$ stands for the $s^{th}$ element of $\pmb{\beta}_{ij}$, $\text{CI}_{ij,s,k}$ stands for the 95% confidence interval obtained from the bootstrap draws associated with $\widetilde{\mathbf{b}}_{ij,s,k}^*-\widetilde{\mathbf{b}}_{ij,s,k}$, and $\widetilde{\mathbf{b}}_{ij,s,k}$ and $\widetilde{\mathbf{b}}_{ij,s,k}^*$ respectively stand for the $s^{th}$ elements of $\widetilde{\mathbf{b}}_{ij}$ and $\widetilde{\mathbf{b}}_{ij}^*$ obtained at the $k^{th}$ simulation replication.

We summarize results in Tables (ref), (ref), and (ref). As $i$ and $j$ are symmetric in a sense, we alter $N$ and $T$ only in the following tables, while keeping $L=60$ for simplicity. We now discuss each table in details. In Table (ref), as expected, as either $N$ or $T$ increases, we have improved probability of correctly identifying $\ell$, $\pmb{\ell}^\circ$, and $\pmb{\ell}^\bullet$. Even for the case $(L, N, T)=(60,60,60)$, we are still able to identify different numbers of factors with probabilities greater than 70%. In Table (ref), as expected, $\text{RMSE}_{\pmb{\beta}}^*$, $\text{RMSE}_{\mathbf{F}}^*$, $\text{RMSE}_{\mathbf{F}^{\circ}}^*$, and $\text{RMSE}_{\mathbf{F}^{\bullet}}^*$ outperform $\text{RMSE}_{\pmb{\beta}} $, $\text{RMSE}_{\mathbf{F}}$, $\text{RMSE}_{\mathbf{F}^{\circ}} $, and $\text{RMSE}_{\mathbf{F}^{\bullet}} $ respectively. However, the differences are not very significant, and gradually vanish as the sample size increases. This is not surprising, as we have better probabilities to identify different types of factors with large sample size. As a result, it is less likely to have misspecified factor structures as the sample size increases. Table (ref) reports the coverage rates of the bootstrap procedure. Overall, the numbers are very close to the nominal rate 95%, and even for $T=60$. To conclude, as explained in Remark (ref), we emphasize that the above results show our method works well for the case with $L\asymp N\asymp T$.

{

table[table omitted — 1,016 chars of source]

table[table omitted — 952 chars of source]

}

{

table[table omitted — 299 chars of source]

}

An Empirical Study

\setcounter{equation}{0}

The analysis of economic and productivity convergence stands as a compelling area of inquiry within economics, offering an ideal application for our hierarchical panel data model as outlined in Section (ref). Specifically, starting with the seminal study by Baumol, numerous studies have been devoted to testing whether income or productivity of poorer economies are converging to those of richer economies. As Durlauf2003 puts it, “Few issues in empirical growth economics have received as much attention as the question of whether countries exhibit convergence". A main technique employed by these studies is “cross-country growth regressions", where aggregate- or industry-level cross-country data are used to regress the average growth rates of per capita income (or labour productivity) over a long period on the initial level of income per capita (or labour productivity) and some additional control variables\footnote{In this study we follow Rodrik2013 and define labour productivity of an industry as the industry's real value added divided by its number of employees. As is well known, the value added of an industry, also referred to as gross domestic product (GDP)-by-industry, is the contribution of a private industry or government sector to overall GDP (The U.S. Bureau of Economic Analysis, 2006, available at \url{https://www.bea.gov/help/faq/184}). The definitions of labour productivity and value added, therefore, imply that an industry's labour productivity can be considered as the industry's “GDP per capita", which in turn implies that neoclassical growth models predict not only conditional convergence in GDP per capita among economies but also industry-level convergence in labour productivity among economies.}. A negative and significant coefficient on the initial conditions is taken to be evident in favour of $\beta $-convergence. For excellent surveys of cross-country convergence studies, see Durlauf2003 for example. Despite the increasing availability of disaggregated data at industry level, the hierarchical structure of these data, to the best of our knowledge, has rarely been explored.

Data

We begin our empirical analysis by introducing the data. For labor productivity (or real value added per employee), we follow Rodrik2013 and use the dataset of UNIDO's INDSTAT2, which provides data on value added (in nominal U.S dollars) and employment for 23 manufacturing industries at the ISIC two-digit level for a large number of countries. With this dataset, real value added can be computed by deflating the nominal value added by the US producer price index, and labor productivity can then be obtained by further dividing real value added by employment (i.e., number of employees). Growth in labor productivity is then measured as percentage change in labor productivity.

Our control variables include a wide range of factors that have been found to be important for assessing convergence. These include human capital (as measured by school enrollment), investment price, trade openness and terms of trade, institutions (measured by civil liberties), and government consumption share. We refer interested readers to Martin for detailed definitions of a wide range of variables. A summary of the dataset is given in Table (ref) and Table (ref). It should be noted that since geographical factors are usually time-invariant, they will be captured by the factor structures and thus are not included in the control variables. When measuring dependent and independent variables, we follow Salimans, and treat them differently. Specifically, the dependent variable is measured as a five-year moving average of economic growth, while all explanatory variables are measured at the beginning of each five year period. Due to data availability, our sample starts at 1963 and ends at 2018. For the same reason, the number of countries varies across industries.

Estimation Results

We now investigate conditional convergence for the manufacturing industries, which can be done by estimating equation (ref) where all components (i.e., initial productivity, control variables and three types of factors) are included. The factor selection is achieved using the method that is described in Section (ref).

For the global factor, we identify one factor that affects the growth in labour productivity of each country in every manufacturing industry and it accounts for 22.50% of the variations in $y_{ijt}$. Panels A and B of Table (ref) present the estimated numbers of industry and country factors. As evident from the table, the number of factors varies among industries and countries, ranging from 1 to 10. These factors collectively explain 56.65% of $y_{ijt}$'s variations. Overall, the factor structure can account for approximately 80% of the variations in $y_{ijt}$, indicating a superior goodness of fit for our model.

Figures (ref) and (ref) provide the boxplots of country-specific estimated coefficients for each explanatory variable across 23 industries. As evident in Figure (ref), the estimated partial effects from the initial productivity are generally negative with an average value of -0.157. This suggests that when country characteristics are controlled, initial labour productivity is generally negatively related to the subsequent rate of growth in labour productivity. In other words, conditional convergence in labour productivity exists for the manufacturing industry. This finding is in line with that of Rodrik2013 who, by applying a fixed effects panel data model to the UNIDO's INDSTAT dataset, also finds conditional convergence in labour productivity for the total manufacturing industry. It is also consistent with the income convergence literature (Martin) that finds that once country characteristics are controlled for, the coefficient on initial income becomes negative and statistically significant. To show the disparity in conditional convergence of productivity among countries, we present the country-specific estimated partial effects for initial productivity within the food and beverages industry in Table (ref) as an example. Despite the consistent negative estimates across all countries, there exists significant divergence in their magnitudes. For instance, coefficients for some countries such as Albania, Croatia, and Malawi exceed -0.25, whereas the effects are substantially lower in others like Pakistan and Peru, or even insignificant as in Jordan and Syria. This finding demonstrates the importance of employing heterogeneous panel data in the examination of productivity convergence. Moreover, such country-specific variations are also present in the estimates for other industries as well as the estimated partial effects for control variables.

In order to confirm our results regarding conditional convergence in labour productivity, we conduct two robustness checks. First, we re-estimate the model for the following two subperiods: 1973-2018 and 1983-2018. The results are provided in Figure (ref). It is evident that the estimated partial effects for initial productivity generally remain negative, confirming the conditional convergence of productivity. Secondly, we follow Rodrik2013 and exclude OCED countries from our sample of countries. The results, which are also presented in Figure (ref), show that our findings of conditional convergence is very robust to the exclusion of OECD countries.

To summarize this section, we have examined the hypothesis of the conditional convergence of productivity for manufacturing industries across countries. The empirical results presented in this section suggest that there is strong and consistent evidence of convergence once factors that affect steady-state levels of labour productivity are controlled for. As comparison, we also examine the unconditional convergence of productivity by excluding the control variables. We provide the estimation results and discussion on such issue in Appendix A.

Conclusion

The hierarchy of large datasets has garnered considerable attention recently (e.g., CYZ2022, JLS2023, JLS2024, Zhang2023, along with numerous real data examples collected in Matyas). In this study, we contribute to the literature by proposing a panel data regression model with three sets of latent factor structures, which, to the best of our knowledge, has not been extensively explored in the existing literature. Rather than consolidating factors from various nodes, as seen in the literature of distributed PCA (Fanetal2019), we propose an estimation method to recover the parameters of interest by peeling off factors layer by layer and across different blocks within certain layers. We establish estimation theory and asymptotic properties accordingly. Additionally, we present a bootstrap procedure to infer the parameters of interest, allowing for different types of cross-sectional dependence (CSD) and time-series autocorrelation (TSA) in the data generating process (DGP). Finally, we validate our theoretical findings using extensive simulated and real data examples. In our empirical study, we utilize data from manufacturing industries at the ISIC two-digit level to examine the twin hypotheses of conditional and unconditional convergence for manufacturing industries across countries.

\setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{section}{0} \setcounter{table}{0} \setcounter{figure}{0} \setcounter{remark}{0} \setcounter{corollary}{0} \setcounter{assumption}{0}

{