EconBase
← Back to paper

Decomposition of Bilateral Trade Flows Using a Three-Dimensional Panel Data Model

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.

74,543 characters · 11 sections · 52 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.

\newtheorem{corollary}{Corollary} \newtheorem{definition}{Definition} \newtheorem{lemma}{Lemma} \newtheorem{proposition}{Proposition} \newtheorem{remark}{Remark} \newtheorem{theorem}{Theorem} \newtheorem{assumption}{Assumption} \newtheorem{example}{Example}

\numberwithin{corollary}{section} \numberwithin{definition}{section} \numberwithin{equation}{section} \numberwithin{lemma}{section} \numberwithin{proposition}{section} \numberwithin{remark}{section} \numberwithin{theorem}{section}

\allowdisplaybreaks[4]

titlepage{ \begin{center} { \bf Decomposition of Bilateral Trade Flows Using a Three-Dimensional Panel Data Model\footnote{ Peng acknowledges the Australian Research Council Discovery Grants Program for its financial support under Grant Number DP210100476. Correspondence: Yufeng Mao, Department of Econometrics and Business Statistics, Monash University, Caulfield East, VIC 3145, Australia. Email: [email removed] }} {\sc Yufeng Mao$^\sharp$, Bin Peng$^\sharp$, Mervyn Silvapulle$^\sharp$, Param Silvapulle$^\sharp$ and Yanrong Yang$^*$} $^\sharp$Monash University and $^*$Australian National University \today \begin{abstract} This study decomposes the bilateral trade flows using a three-dimensional panel data model. Under the scenario that all three dimensions diverge to infinity, we propose an estimation approach to identify the number of global shocks and country-specific shocks sequentially, and establish the asymptotic theories accordingly. From the practical point of view, being able to separate the pervasive and nonpervasive shocks in a multi-dimensional panel data is crucial for a range of applications, such as, international financial linkages, migration flows, etc. In the numerical studies, we first conduct intensive simulations to examine the theoretical findings, and then use the proposed approach to investigate the international trade flows from two major trading groups (APEC and EU) over 1982-2019, and quantify the network of bilateral trade. \end{abstract} \end{center} {\em Keywords}: Three-Dimensional Panel Data, Bilateral Trade, Asymptotic Theory {\em JEL classification}: C23, P45 }

Introduction

All countries of the world are nowadays connected with each other more or less through varieties of bilateral trade. Getting reliable and up-to-date statistics on exports and imports of different countries is thus crucial in order to provide a detailed insight into the most recent trading patterns. Given an increasing interest in understanding such a complex network, we see the rising popularity of multi-dimensional models over the past decade, e.g., MNP2013, BHM2016, andreou2019inference, choi2020canonical, kapetanios2020estimation, just to name a few. Excellent reviews on the applications and theoretical developments of multi-dimensional panel data models can be respectively seen in BEP2016 and Jorg2016 for instance.

Despite a vast amount of research on two dimensional factor models (see BN2008 for an excellent review), it seems that the literature of multi-dimensional models has not even settled on how to effectively distinguish pervasive and nonpervasive economic shocks, where pervasive and nonpervasive shocks refer to those affecting the entire network and those affecting only a part of the network respectively (e.g., Wang2008, ER2017). In this regard, the presentation (ref) of Section (ref) provides a clear visualization using matrix form.

To solve the aforementioned issue, different algorithms have been proposed (e.g., Jorg2016 and references therein), but from the theoretical point of view the progress has not been pushed forward much since Wang2008. We now comment on the relevant literature. ER2017 extend the study of Wang2008 to allow for long run dependence, and both papers numerically rely on some initial estimates on the pervasive and nonpervasive factors. In our view, the requirement on initial estimates is due to the fact that both studies aim to estimate the pervasive and nonpervasive factors in one objective function, which as a consequence leads to a complex minimization problem. Thus, the numerical implementation often becomes complex, and is hard to be justified. In another two works, both CKKK2018 and andreou2019inference propose sequential procedures to identify and estimate pervasive and nonpervasive shocks, in which canonical correlation analysis (CCA) are adopted. However, only two of the three dimensions are allowed to diverge in both studies. Han2019 considers a shrinkage estimation approach to explore the group effects of the factor structure, which can be computationally expensive, as the choice of tuning parameter often plays an important role in practice.

From the practical point of view, being able to separate the pervasive and nonpervasive shocks in a multi-dimensional panel data is crucial for a range of applications. First, as mentioned in the beginning of the paper, accounting for pervasive and nonpervasive shocks reveals a detailed network structure of the international trade. We will come back to it in the empirical study section. A second example is better understanding business-cycle fluctuations across countries and regions (KOW2003). Along this line of research, identifying the common fluctuations across macroeconomic aggregates worldwide has always been one of the priorities (e.g., GHR1997). The emergence of multi-dimensional panel data models provides an excellent framework to facilitate the investigation. Another field which urgently calls for development on multi-dimensional panel data models is associated with migration flows. As well understood, the rate of migration between two countries does not depend solely on their relative attractiveness, but also on the one of alternative destinations (BERTOLI201379). Given the increasing mobility of the entire population, how to better capture the bilateral flows therefore becomes vital now more than ever. Other examples requiring multi-dimensional panel data models can also be found in CKKK2018, kapetanios2020estimation, etc.

Having presented the above challenges and necessities, in this study, we specifically consider a three-dimensional panel data model with unobserved global (pervasive) and country-specific (nonpervasive) factors, which has been exposed in the literature but has not been fully solved to the best of the authors' knowledge. On theory, our contributions are the following three-fold: (1). under the scenario that all three dimensions can diverge to infinity, we propose an estimation approach to identify the number of global shocks and country-specific shocks sequentially; (2). the newly proposed approach is easy to implement, and the asymptotic theories are established accordingly; (3). we further conduct intensive numerical studies to examine the finite sample performance of the newly proposed approach using both simulated and real datasets. In the empirical study, we then apply the approach to decompose the network of bilateral trade using country level data from two major trading groups (APEC and EU) over the period 1982-2019. We find that the country-specific shocks become more volatile in recent years, which may indicate the increasing instability of the inward and outward bilateral trade costs over the past couple of decades. In addition, we show that the trade flows involving China mainland, Germany and the United States show relatively strong sensitivity to global shocks, which reflects the fact that, in general, they are leading export and import countries worldwide. We note that the relationship among Canada, Mexico, and the United States is also highly sensitive to different shocks, which somewhat reflects the fact that all three of them are highly economically related through North American Free Trade Agreement (NAFTA) that eliminates some trade barriers and promotes the trading activities.

The structure of this paper is as follows. Section (ref) presents the model with the estimation approach, and establishes the asymptotic properties accordingly. In Section (ref), we conduct intensive simulations to examine the finite sample performance of the newly proposed approach. Section (ref) provides an empirical study using country level bilateral trade data. Section (ref) concludes. Due to the limit of space, the preliminary lemmas and the proofs are given in the online supplementary appendices.

Before proceeding further, it is convenient to introduce some notation: $\| \cdot\|_{F}$ denotes the Euclidean norm of a vector or the Frobenius norm of a matrix; for a matrix $\textbf{A}$, its spectral norm is defined as $\| \textbf{A} \|_2 = \sqrt{\lambda_{\max} \{ \textbf{A}^{\prime} \textbf{A} \} }$, where $\lambda_{\max} \{ \cdot \}$ denotes the maximum eigenvalue; $\textbf{M}_{\textbf{A}} = \textbf{I} - \textbf{P}_{\textbf{A}} $ denotes the orthogonal projection matrix generated by matrix $\textbf{A}$, where $ \textbf{P}_{\textbf{A}} = \textbf{A}(\textbf{A}' \textbf{A})^{-1} \textbf{A}'$ and $\textbf{A}$ is a matrix with full column rank; let $\to_P$ and $\to_D$ denote convergence in probability and in distribution, respectively; we write $a \asymp b$ if $a=O_P(b)$ and $b=O_P(a)$; let $\operatorname*{\normalfont\textrm{diag}}(\textbf{A}, \textbf{B})$ denotes the block-diagonal matrix that takes $\textbf{A}$ and $\textbf{B}$ as the upper left and lower right blocks; $\mbox{vec}(\textbf{A})$ stands for the vectorization operation; $\mathbb{I}(\cdot)$ stands for the indicator function.

Model & Methodology

In this section, we first present the model, then provide the estimation approach, and finally establish the asymptotic theories accordingly.

The Setup

Having presented our motivations in Section (ref), we specifically consider the next model in this study.

eqnarray[eqnarray omitted — 154 chars of source]

where $i=1,\ldots,M$ index the exporters, $j=1,\ldots,N$ index the importers, and $t=1,\ldots,T$ index the time periods. We observe $y_{ijt}$'s only, and $u_{ijt}$'s are the idiosyncratic error terms. $\bm{g}_t$ is an $r_g\times 1$ unobservable global factor, which is regarded as global shocks and may capture the globalisation trends. Some detailed explanation on the globalisation trends can be found in kapetanios2020estimation, and we shall be more specific on this so-called “trend" in the empirical study of Section (ref). $\bm{f}_{E,it}$ and $\bm{f}_{I,jt}$ represent the unobservable $r_{E,i}\times 1$ and $r_{I,j}\times 1$ country-specific factors. Specifically, $\bm{f}_{E,it}$ is referred to as an exporter factor which affects all import partners associated with export country $i$ and $\bm{f}_{I,jt}$ is referred to as an importer factor which affects all export partners associated with import country $j$. The country-specific factors may capture the unobservable multilateral trade resistances (MTRs) that are different for exporters and importers. Loosely speaking, MTRs refer to the barriers which each of exporter and importer face in their trade with all their trading partners. We refer interested readers to anderson2003gravity for a comprehensive discussion on MTR. $\bm{\gamma}_{ij}$, $\bm{\lambda}_{E,ij}$ and $\bm{\lambda}_{I,ij}$ are the corresponding factor loadings. Throughout this paper, we always use the subscript $_g$ to denote the variables associated with the global factors, and use the subscripts $_E$ and $_I$ to denote the variables associated with the exporters and importers respectively.

The model (ref) is in fact not new, and has been mentioned in Jorg2016, CKKK2018, and \citet*[eq. 2]{kapetanios2020estimation} among others for different purposes. In what follows, we propose an easily implemented methodology to recover the structure of the right hand side of (ref), when all three dimensions are allowed to diverge to infinity. Precisely, we first estimate the numbers of global and country-specific factors (i.e., the values of $r_g$, $r_{E,i}$'s and $r_{I,j}$'s), and then establish inferences for global and country-specific shocks.

remarkBefore proceeding further, we comment on an important identification issue. For simplicity, we suppose that $r_g =1$ and $g_t \equiv 1$, and suppose further that \begin{eqnarray*} &&\bm{f}_{E, it} =\bm{f}_E +\bm{\eta}_{E, it}\quadwith\quad E[\bm{\eta}_{E, it} ]=0,\nonumber \\ &&\bm{f}_{I, it} =\bm{f}_I +\bm{\eta}_{I, it}\quadwith\quad E[\bm{\eta}_{I, it} ]=0. \end{eqnarray*} Then, the model (ref) becomes \begin{eqnarray} y_{ijt} = \bm{\gamma}_{ij}^* + \bm{\lambda}_{E, ij}' \bm{\eta}_{E, it} + \bm{\lambda}_{I,ij}' \bm{\eta}_{I,j t} + u_{ijt}, \end{eqnarray} where $ \bm{\gamma}_{ij}^* = \bm{\gamma}_{ij}+\bm{\lambda}_{E,ij}' \bm{f}_{E} + \bm{\lambda}_{I,ij}' \bm{f}_{I} $. It then infers that for a model having a multi-layer factor structure, only one layer can have non-zero mean factors.

Having said Remark (ref), without loss of generality, we assume that

eqnarray[eqnarray omitted — 102 chars of source]

for country-specific factors throughout this study.

As repeatedly pointed out in the literature (e.g., Wang2008, Jorg2016, CKKK2018), investigating (ref) relies on how to utilize the sparse structure of the next presentation.

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

In view of (ref), a few facts emerge:

enumerate• In order to estimate (ref), one needs to identify the number of factors for each $\bm{g}_t$, $\bm{f}_{E,it}$ and $\bm{f}_{I,jt}$. Traditional PCA usually requires a low rank setting. However, having $\bm{g}_t$, $\bm{f}_{E,it}$'s and $\bm{f}_{I,jt}$'s in one column as in (ref) yields a factor with a diverging dimension, which suggests that recovering all factors and loadings in one goal seems to be challenging. Thus, it motivates us to consider a multiple steps approach below. • The country-specific factors associated with exporters and importers are interchangeable, as the sparse structure associated with the corresponding factor loadings depends on how we rank $y_{ijt}$ with respect to $i$ and $j$ only. Thus, we would expect to recover the exporter and importer factors in a parallel manner. • As clearly seen in (ref), $\bm{g}_t$ has an impact on every single $y_{ijt}$, although the magnitude depends on the value of $\bm{\gamma}_{ij}$. However, $\bm{f}_{E,it}$ or $\bm{f}_{I,jt}$ affects only an asymptotically negligible subset of $y_{ijt}$'s due to the sparse structure. From the signal-to-noise ratio point of view, we expect that the global factors are easier to be identified. Intuitively speaking, they can be estimated first if principal component analysis (PCA) is employed. As the country-specific factors contain the second tier of signal, they should be recovered after removing the dominating ones.

In Section (ref) below, we propose a multi-step estimation approach based on the aforementioned points.

The Estimation Approach

We are now ready to present the estimation approach, which is a procedure involving multiple steps. The outline is as follows.

enumerate• Conduct PCA to identify the number of global factors $r_g$, and estimate the global factor, which contains the strongest “signal" as explained under (ref). • Remove the estimated global factor, then simultaneously conduct multiple PCA to estimate the number of country-specific factors $r_{E,i}$' and $r_{I,j}$'s, and recover the country-specific factors, which contain “signals" weaker than the global factor but stronger than the error terms.

First, we write (ref) in matrix form to facilitate the development. Throughout, the subscript $_\bullet$ always stands for including all available sample in the corresponding dimension for notational simplicity.

eqnarray[eqnarray omitted — 135 chars of source]

where the response variables and error terms are defined by

eqnarray[eqnarray omitted — 413 chars of source]

the global factors and loadings are defined by

eqnarray[eqnarray omitted — 185 chars of source]

and the country-specific factors and loadings are defined by

eqnarray[eqnarray omitted — 799 chars of source]

With the above notations in hand, we are ready to present the details of each step with necessary discussions.

Step 1 --- Conduct PCA on $\frac{1}{MNT}\textbf{Y}'\textbf{Y} $ as follows.

eqnarray[eqnarray omitted — 122 chars of source]

in which $\frac{1}{T}\widehat{\textbf{G}}'\widehat{\textbf{G}} =\textbf{I}_{k_{\max}}$, $\textbf{V}_{g} =\operatorname*{\normalfont\textrm{diag}}\{ \widehat{\rho}_{g,1},\ldots, \widehat{\rho}_{g,k_{\max}} \}$ with $ \widehat{\rho}_{g,1}\ge \cdots \ge \widehat{\rho}_{g,k_{\max}}$ being the largest $k_{\max}$ eigenvalues, $k_{\max} \ (> r_g)$ is a user-specified fixed large integer. By (ref), we implement the following two sub-steps.

enumerate• Estimate the number of global factors $r_g$ by \begin{equation} \widehat{r}_g = \operatorname*{\arg\!\min}_{0 \le k \le k_{\max}} \left\{ \frac{ \widehat{\rho}_{g,k+1}}{ \widehat{\rho}_{g,k}} \cdot \mathbb{I} \left( \widehat{\rho}_{g,k}\ge\omega_{MNT} \right) + \mathbb{I} \left( \widehat{\rho}_{g,k}< \omega_{MNT} \right) \right\}, \end{equation} where $\omega_{MNT} = 1/\ln (\max\{M, N, T\})$, and $\widehat{\rho}_{g,0} =1$ is a mock eigenvalue. • Estimate $\textbf{G}$ by letting $\widehat{\textbf{G}}$ include the first $\widehat{r}_g$ columns only, where we have slightly abused the notation $\widehat{\textbf{G}}$. The loading matrix is estimated by $\widehat{\bm{\Gamma}} =\frac{1}{T} \textbf{Y} \widehat{\textbf{G}} $.

Step 2 includes two parallel sections: Part 1 and Part 2.

Part 1 --- For each $j=1,\ldots,N$, conduct PCA:

eqnarray[eqnarray omitted — 259 chars of source]

where $ \textbf{Y}_{I,j} = (\textbf{Y}_{\bullet j1},\ldots, \textbf{Y}_{\bullet jT} )$ with $ \textbf{Y}_{\bullet jt} =(y_{1jt},\ldots, y_{Mjt})'$, $\widehat{\bm{\Gamma}}_{I, j}$ includes the $M$ rows of $\widehat{\bm{\Gamma}}$ corresponding the $j^{th}$ importer, $\frac{1}{T}\widehat{\textbf{F}}_{I,j}' \widehat{\textbf{F}}_{I,j}=\textbf{I}_{k_{\max}}$, and $\textbf{V}_{I,j} =\operatorname*{\normalfont\textrm{diag}}\{\widehat{\rho}_{Ij,1},\ldots, \widehat{\rho}_{Ij,k_{\max}} \}$ with $\widehat{\rho}_{Ij,1}\ge \cdots \ge \widehat{\rho}_{Ij,k_{\max}}$ being the largest $k_{\max}\ (> r_{I,j})$ eigenvalues. By (ref), implement the followings.

enumerate• Estimate $r_{I,j}$ by \begin{equation} \widehat{r}_{I,j} = \operatorname*{\arg\!\min}_{0 \le k \le k_{\max}} \left\{ \frac{ \widehat{\rho}_{Ij,k+1}}{ \widehat{\rho}_{Ij,k}} \cdot \mathbb{I} \left( \widehat{\rho}_{Ij,k}\ge \omega_{MNT} \right) + \mathbb{I} \left( \widehat{\rho}_{Ij,k} < \omega_{MNT} \right) \right\}, \end{equation} where $\widehat{\rho}_{Ij,0} =1$ is a mock eigenvalue. • Estimate $\textbf{F}_{I,j} = ( \bm{f}_{I,j1},\ldots, \bm{f}_{I,jT})'$ by letting $\widehat{\textbf{F}}_{I,j}$ include the first $\widehat{r}_{I,j}$ columns only. The loading matrix $\bm{\Lambda}_{I,\bullet j}$ defined in (ref) is estimated by $\widehat{\bm{\Lambda}}_{I, \bullet j} = \frac{1}{T}(\textbf{Y}_{I, j} - \widehat{\bm{\Gamma}}_{I, j} \widehat{\textbf{G}}' ) \widehat{\textbf{F}}_{I,j} $.

Part 2 --- For each $i=1,\ldots,M$, conduct PCA:

eqnarray[eqnarray omitted — 261 chars of source]

where $ \textbf{Y}_{E,i} = (\textbf{Y}_{i \bullet 1},\ldots, \textbf{Y}_{i\bullet T} )$ with $\textbf{Y}_{i\bullet t} =(y_{i1t},\ldots, y_{iNt})'$, $\widehat{\bm{\Gamma}}_{E, i}$ includes the $N$ rows of $\widehat{\bm{\Gamma}}$ corresponding to the $i^{th}$ exporter, $\frac{1}{T}\widehat{\textbf{F}}_{E,i}' \widehat{\textbf{F}}_{E,i}=\textbf{I}_{k_{\max}}$, and $\textbf{V}_{E,i} =\operatorname*{\normalfont\textrm{diag}}\{\widehat{\rho}_{Ei,1},\ldots, \widehat{\rho}_{Ei,k_{\max}} \}$ with $\widehat{\rho}_{Ei,1}\ge \cdots \ge \widehat{\rho}_{Ei,k_{\max}}$ being the largest $k_{\max}\ (> r_{E,i})$ eigenvalues. By (ref), we conduct the followings.

enumerate• Estimate $r_{E,i}$ by \begin{equation} \widehat{r}_{E,i} = \operatorname*{\arg\!\min}_{0 \le k \le k_{\max}} \left\{ \frac{\widehat{\rho}_{Ei,k+1}}{\widehat{\rho}_{Ei,k}} \cdot \mathbb{I} \left( \widehat{\rho}_{Ei,k}\ge \omega_{MNT} \right) + \mathbb{I} \left( \widehat{\rho}_{Ei,k} < \omega_{MNT} \right) \right\}, \end{equation} where $\widehat{\rho}_{Ei,0} =1$ is a mock eigenvalue. • Estimate $\textbf{F}_{E,i} = (\bm{f}_{E,i1},\ldots, \bm{f}_{E,iT})'$ by letting $\widehat{\textbf{F}}_{E,i}$ include the first $\widehat{r}_{E,i}$ columns only. The loading matrix $\bm{\Lambda}_{E,i\bullet} =(\bm{\lambda}_{E,i1}, \ldots, \bm{\lambda}_{E,iN})'$ is estimated by $\widehat{\bm{\Lambda}}_{E,i \bullet} = \frac{1}{T}(\textbf{Y}_{E, i} - \widehat{\bm{\Gamma}}_{E, i} \widehat{\textbf{G}}' ) \widehat{\textbf{F}}_{E,i} $.
remarkWe make a few comments on the estimation approach. (1). The use of eigenvalue ratio in (ref), (ref) and (ref) is in the same spirit of LamYao2012 and AhnHorenstein2013. (2). The threshold $\omega_{MNT} $ is to bypass a technical challenge raised in LamYao2012, and the mock eigenvalues $\widehat{\rho}_{g,0}$, $\widehat{\rho}_{Ij,0}$'s and $\widehat{\rho}_{Ei,0}$'s are designed to capture the cases where there are no global factors, or some of the country-specific factors do not exist. From the dimension reduction point of view, it is crucial to have a procedure which accounts for zero factors under the three dimensional panel data framework. (3). $k_{\max}$ is a user-defined fixed integer. Practically, one can adopt any reasonable large value which suits the empirical study (e.g., FLM13, PX2019).

Consistency

In this subsection, we show that the number of factors can be identified consistently in each step with necessary conditions. The asymptotic distributions are established in the next subsection.

To facilitate the development, we impose the following conditions.

assumption\begin{enumerate} • As $T\to \infty$, $\normalfont \frac{1}{T}\textbf{G}' \textbf{G}\to_P \bm{\Sigma}_{\textbf{G}}$, where $\normalfont \bm{\Sigma}_{\textbf{G}}$ is a deterministic positive definite matrix. Also, $\max_{t\ge 1}E \| \bm{g}_t \|_{F}^4 < \infty$. • Suppose that (ref) holds. Moreover, $\max_{i\ge 1,t\ge 1}E\| \bm{f}_{E, it} \|_F^4 < \infty$ and $\|\normalfont\textbf{F}_{E}\|_2 = O_P(\sqrt{T} \vee \sqrt{M})$. Also, $\max_{j\ge 1,t\ge 1}E\| \bm{f}_{I, jt} \|_F^4 < \infty$ and $\|\normalfont\textbf{F}_{I}\|_2 = O_P(\sqrt{T} \vee \sqrt{N})$. \end{enumerate}
assumption\begin{enumerate} • As $(M, N) \to (\infty,\infty)$, $\frac{1}{MN}\bm{\Gamma}' \bm{\Gamma} \to_P \bm{\Sigma}_{\bm{\Gamma}}$, where $\bm{\Sigma}_{\bm{\Gamma}}$ is a deterministic positive definite matrix. Also, $\max_{i\ge1, j\ge 1}E \| \bm{\gamma}_{ij} \|_{F}^4 < \infty$. • Suppose that $\max_{i\ge1,j\ge 1}E\| \bm{\lambda}_{E, ij}\|_F^4<\infty$ and $\max_{i\ge1,j\ge 1}\| \bm{\lambda}_{E, ij}\|_F = O_P(\sqrt{\ln (MN)})$. Also, $\max_{i\ge1,j\ge 1}E\| \bm{\lambda}_{I, ij}\|_F^4<\infty$ and $\max_{i\ge1,j\ge 1}\| \bm{\lambda}_{I, ij}\|_F = O_P(\sqrt{\ln (MN)})$. \end{enumerate}
assumption\begin{enumerate} • Let $\{u_{ijt}\ |\ i\ge 1,j\ge 1, t\ge 1\}$ be independent of the other variables. Let $\mathcal{F}_{-\infty}^0$ and $\mathcal{F}_\tau^\infty$ denote the $\sigma$-algebras generated by $\normalfont\{\textbf{U}_{\bullet \bullet t} \ | \ t \le 0\}$ and $\normalfont\{\textbf{U}_{\bullet \bullet t} \ | \ t \ge \tau\}$ respectively, where $\normalfont \textbf{U}_{\bullet \bullet t} = (u_{11t},\ldots, u_{M1t},\ldots,u_{1Nt},\ldots, u_{MNt})'$. Define the mixing coefficient $\alpha(\tau) = \sup_{A\in \mathcal{F}_{-\infty}^0, B\in \mathcal{F}_\tau^\infty} \left|\Pr(A)\Pr(B) -\Pr(AB) \right|$. \begin{enumerate} • Let $\normalfont\{\textbf{U}_{\bullet \bullet t} \ | \ t\ge 1\}$ be strictly stationary and $\alpha$-mixing such that for some $\nu>0$, $ \max_{i\ge 1,j\ge 1} E|u_{ijt}|^{4+\nu} <\infty$, and the mixing coefficient satisfies $ \sum_{t=1}^\infty [\alpha(t)]^{\nu/(2+\nu)}$ $< \infty$. • $E[u_{ijt}]=0$, $ \max_{i\ge 1, j\ge 1}\sigma_{ij}^2<\infty$ and $ \sum_{(i,j) \neq (m,n)} |\sigma_{ijmn}|=O(MN)$, where $\sigma_{ij}^2 = E[u_{ijt}^2]$ and $\sigma_{ijmn} =E[u_{ijt}u_{mnt}]$ for $t\ge 1$. In addition, suppose that\\ $\sum_{i,m=1}^M \sum_{j,n=1}^N \sum_{t,s=1}^T |E[ u_{ijt}u_{mns}] |=O(MNT)$. \end{enumerate} • Suppose that $r_g<\infty$, $ \max_{i\ge 1}r_{E,i}<\infty$, and $ \max_{j\ge 1}r_{I,j}<\infty$. \end{enumerate}

Assumption (ref) imposes restrictions on the global and country-specific factors, which are not more restrictive than Assumption 1.i of CKKK2018. The conditions on the spectral norm of $\textbf{F}_E$ and $\textbf{F}_I$ are widely adopted in the literature (e.g., \citealp*[Assumption A.1.iii]{LiQianSu} and LuSu). Extensive discussions with examples on this type of assumption can be found in Moon.

Assumption (ref) puts restrictions on the loadings associated with the global and country-specific factors. The bounds on $\max_{i\ge1,j\ge 1}\| \bm{\lambda}_{E, ij}\|_F$ and $\max_{i\ge1,j\ge 1}\| \bm{\lambda}_{I, ij}\|_F$ are fairly standard. See Assumption A7 of CHL2012 for example.

Assumption (ref).1 assumes that the error terms $u_{ijt}$'s follow stationary time series process over $t$, and simultaneously allow for weak cross-sectional dependence over $i$ and $j$. Assumption (ref).2 requires $r_g$, $r_{E,i}$'s and $r_{I,j}$'s to be bounded, which nests $r_g=0$, $r_{E,i}=0$ and $r_{I,j}=0$ as special cases.

Under these conditions, we present the first theorem of this paper below.

theoremUnder Assumptions (ref)-(ref), as $(M,N,T) \to (\infty,\infty,\infty)$, \begin{enumerate} • in {\normalfont Step 1.1}, $\Pr(\widehat{r}_g = r_g) \to 1$; • in {\normalfont Step 1.2}, $\normalfont\frac{1}{\sqrt{T}} \| \widehat{\textbf{G}} - \textbf{G} \textbf{H} \|_{F} = O_P \left( \frac{\sqrt{ \ln (MN)}}{\min \{\sqrt{M}, \sqrt{N}, \sqrt{T}\}} \right)$, where $\normalfont\textbf{H} = \frac{1}{MN}\bm{\Gamma}' \bm{\Gamma}\cdot \frac{1}{T}\textbf{G}' \widehat{\textbf{G}} \cdot(\textbf{V}_{g}^{\dag})^{-1}$, and $\normalfont\textbf{V}_{g}^{\dag}$ is the $r_g\times r_g$ leading principal submatrix of $\normalfont \textbf{V}_g$. \end{enumerate}

Theorem (ref).1 shows that $r_g$ can be estimated consistently, while Theorem (ref).2 indicates that we can only recover $\textbf{G}$ up to a rotation matrix. From the signal-to-noise ratio point of view, only the space spanned by the global factors can be recovered in Step 1.

remarkIt is noteworthy that when establishing Theorem (ref), no harsh conditions are imposed between the global factor structure and the country-specific ones. In this sense, although the rate of Theorem (ref) is slow, we show that the global factors can be identified from the data first with minimum cost. In the traditional literature, the fact has barely been mentioned. To the best of the authors' knowledge, the only exception is Remark 4 of Han2019. In Appendix A of the online supplementary file, we provide a sharper rate for the estimation of the global factor when more structures are adopted. The details are summarized in Lemma (ref).

Having presented the results associated with the global factors, we investigate the country-specific ones, and further impose the following conditions.

assumption\begin{enumerate} • For $i=1,\ldots,M$ and $j=1,\ldots,N$, suppose that the following conditions hold: \begin{enumerate} • $\normalfont \frac{1}{T} \| \textbf{G}'\textbf{F}_{E, i} \|_F = O_P(T^{a_{E,i}})$ and $\normalfont\frac{1}{T} \| \textbf{G}'\textbf{F}_{I,j} \|_F = O_P(T^{a_{I,j}})$, where $\normalfont\textbf{F}_{E, i}$ and $\normalfont\textbf{F}_{I,j}$ are defined in (ref), $\max_{i\ge 1}a_{E,i}< 0$, and $\max_{j\ge 1}a_{I,j}< 0$; • $\max_{i\ge 1}\|\normalfont\frac{1}{T}\textbf{F}_{E,i}' \textbf{F}_{E, i} - \bm{\Sigma}_{\textbf{F}_{E,i}} \|_F =o_P(1)$ and $\max_{j\ge 1}\|\normalfont\frac{1}{T}\textbf{F}_{I, j}' \textbf{F}_{I, j} - \bm{\Sigma}_{\textbf{F}_{I,j}} \|_F =o_P(1)$, where $\bm{\Sigma}_{\textbf{F}_{E,i}}$ and $\bm{\Sigma}_{\textbf{F}_{I,j}}$ are deterministic positive definite matrices; • $\normalfont \max_{i\ge 1}\|\frac{1}{N}\bm{\Lambda}_{E, i\bullet }' \bm{\Lambda}_{E, i\bullet } - \bm{\Sigma}_{\bm{\Lambda}_{E,i\bullet }}\|_F =o_P(1)$ and $\normalfont \max_{j\ge 1}\|\frac{1}{M}\bm{\Lambda}_{I,\bullet j}' \bm{\Lambda}_{I,\bullet j} - \bm{\Sigma}_{\bm{\Lambda}_{I,\bullet j}}\|_F =o_P(1)$, where $\normalfont\bm{\Lambda}_{E,i\bullet } = (\bm{\lambda}_{E,i1}, \ldots, \bm{\lambda}_{E,iN})'$, $\normalfont\bm{\Lambda}_{I,\bullet j}$ is defined under (ref), and $\bm{\Sigma}_{\bm{\Lambda}_{E,i\bullet }}$ and $\bm{\Sigma}_{\bm{\Lambda}_{I,\bullet j}}$ are deterministic positive definite matrices. \end{enumerate} • Suppose that $\max_{j\ge 1} \sum_{i \ne m} \sigma_{ijmj}=O(M)$, and $\max_{i\ge 1} \sum_{j\ne n} \sigma_{ijin}=O(N)$, where $\sigma_{ijmn}$ is defined in Assumption (ref). \end{enumerate}

Assumption (ref).1.(a) requires certain orthogonality between the global factors and country-specific factors. Specifically, the values of $a_{E,i}$ and $a_{I,j}$ measure the degree of orthogonality between the global and country-specific factors. If $a_{E,i}=a_{I,j} =-\infty$, this condition essentially reduces to Assumption A of Ando, where they show the necessity of orthogonality in order to identify the common and group-specific factors under a two-dimensional panel data framework. Similar discussions on orthogonality can also be seen in andreou2019inference. Assumptions (ref).1.(b) and (ref).1.(c) impose more conditions on the blocks of factors and loadings associated with exporters and importers, which are fairly standard. Assumption (ref).2 further regulates the weak cross-sectional dependence of the error terms.

With Assumption (ref) in hand, the country-specific factor structures can be successfully recovered in Step 2. The details are summarized in the next theorem.

theoremUnder Assumptions (ref)-(ref), as $(M,N,T) \to (\infty,\infty,\infty)$, \begin{enumerate} • For $j=1,\ldots,N$, \begin{enumerate} • in {\normalfont Part 1.1 of Step 2}, $\Pr(\widehat{r}_{I,j}= r_{I,j}) \to 1$; • in {\normalfont Part 1.2 of Step 2}, $\normalfont\frac{1}{\sqrt{T}} \| \widehat{\textbf{F}}_{I,j} - \textbf{F}_{I,j} \textbf{H}_{I,j} \|_{F} = O_P \left( \frac{\sqrt{ \ln (MN)}}{\min \{\sqrt{M}, \sqrt{N}, \sqrt{T}\}} + T^{a_{I,j}} \right)$, where $\normalfont\textbf{H}_{I,j} = \frac{1}{M}\bm{\Lambda}_{I, \bullet j}' \bm{\Lambda}_{I, \bullet j}\cdot \frac{1}{T}\textbf{F}_{I,j}' \widehat{\textbf{F}}_{I,j} \cdot(\textbf{V}_{I,j}^{\dag})^{-1}$, and $\normalfont\textbf{V}_{I,j}^{\dag}$ is the $r_{I,j}\times r_{I,j}$ leading principal submatrix of $\normalfont \textbf{V}_{I,j}$. \end{enumerate} • For $i=1,\ldots,M$, \begin{enumerate} • in {\normalfont Part 2.1 of Step 2}, $\Pr(\widehat{r}_{E,i}= r_{E,i}) \to 1$; • in {\normalfont \textbf{Part 2.2} of \textbf{Step 2}}, $\normalfont\frac{1}{\sqrt{T}} \| \widehat{\textbf{F}}_{E,i} - \textbf{F}_{E,i} \textbf{H}_{E,i} \|_{F} = O_P \left( \frac{\sqrt{ \ln (MN)}}{\min \{\sqrt{M}, \sqrt{N}, \sqrt{T}\}} +T^{a_{E,i}} \right)$, where $\normalfont\textbf{H}_{E,i} = \frac{1}{N}\bm{\Lambda}_{E, i\bullet}' \bm{\Lambda}_{E,i \bullet }\cdot \frac{1}{T}\textbf{F}_{E,i}' \widehat{\textbf{F}}_{E,i} \cdot(\textbf{V}_{E,i}^{\dag})^{-1}$, and $\normalfont\textbf{V}_{E,i}^{\dag}$ is the $r_{E,i}\times r_{E,i}$ leading principal submatrix of $\normalfont \textbf{V}_{E,i}$. \end{enumerate} \end{enumerate}

Theorem (ref) shows that $r_{E,i}$ and $r_{I,j}$ can be estimated consistently. Moreover, $\widehat{\textbf{F}}_{I,j}$ and $\widehat{\textbf{F}}_{E,i}$ respectively recover $\textbf{F}_{I,j}$ and $\textbf{F}_{E,i}$ up to rotation matrices.

Till now, we conclude that we have successfully recovered the network presented by (ref). To establish inferences for the estimation approach, we study the asymptotic distributions associated with Step 1 and Step 2 in the next subsection.

Asymptotic Distribution

In order to establish the asymptotic distributions, the following assumptions are necessary to facilitate the development.

assumption\begin{enumerate} • Let $\normalfont \frac{1}{\sqrt{MN}} \| \bm{\Gamma}'\bm{\Lambda}_{E} \|_F = O_P(1)$ and $\normalfont\frac{1}{\sqrt{MN}} \| \bm{\Gamma}'\bm{\Lambda}_{I} \|_F = O_P(1)$. • $\normalfont \frac{1}{T} \textbf{G}'\textbf{G} = \textbf{I}_{r_g}$ and $\normalfont \bm{\Gamma}' \bm{\Gamma}$ is a diagonal matrix with distinct entries. • Suppose that $\normalfont \frac{1}{\sqrt{MN}} \sum_{i=1}^{M} \sum_{j=1}^{N} \bm{\gamma}_{ij} \upsilon_{ijt} \to_{D} N(\textbf{0}, \bm{\Phi}_{t})$ for $t=1,\ldots, T$, where $\upsilon_{ijt} = \bm{\lambda}_{E, ij}' \bm{f}_{E, it} + \bm{\lambda}_{I,ij}' \bm{f}_{I,j t} + u_{ijt}$. \end{enumerate}
assumption\begin{enumerate} • Suppose that $\normalfont \frac{1}{T} \| \textbf{F}_{E, i}'\textbf{F}_{I, j} \|_F = O_P(T^{b_{EI, ij}})$, where $\max_{i\ge 1,j\ge1} b_{EI,ij} < 0$. • \begin{enumerate} • $\normalfont \frac{1}{T} \textbf{F}_{I,j}'\textbf{F}_{I,j} = \textbf{I}_{r_{I,j}}$ and $\normalfont \bm{\Lambda}_{I,\bullet j}' \bm{\Lambda}_{I,\bullet j}$ is a diagonal matrix with distinct entries; • $\normalfont \frac{1}{T} \textbf{F}_{E,i}'\textbf{F}_{E,i} = \textbf{I}_{r_{E,i}}$ and $\normalfont \bm{\Lambda}_{E, i \bullet}' \bm{\Lambda}_{E, i \bullet}$ is a diagonal matrix with distinct entries. \end{enumerate} • \begin{enumerate} • $\normalfont \frac{1}{\sqrt{M}} \sum_{i=1}^{M} \bm{\lambda}_{I,ij} (\bm{\lambda}_{E, ij}' \bm{f}_{E, it} + u_{ijt}) \to_{D} N(\textbf{0}, \bm{\Omega}_{I,jt})$ for each pair of $(j,t)$; • $\normalfont \frac{1}{\sqrt{N}} \sum_{j=1}^{N} \bm{\lambda}_{E,ij} (\bm{\lambda}_{I, ij}' \bm{f}_{I, jt} + u_{ijt}) \to_{D} N(\textbf{0}, \bm{\Omega}_{E,it})$ for each pair of $(i,t)$. \end{enumerate} \end{enumerate}

Assumption (ref).1 requires certain orthogonality between global factor loadings and country-specific factor loadings, which is not unusual in the literature. For instance, LamYao2012 explain the rational behind such a setting at length. Assumption (ref).2 further imposes conditions for the purpose of identification, which has been extensively discussed in BN2013 and FanLiaoWang. In view of Remark (ref), Assumption (ref).3 is fairly standard. We further explain Assumption (ref).3 together with Assumption (ref).3 below.

Similar to Assumption (ref).1, Assumption (ref).1 requires certain orthogonality but focusing on the export factors and importer factors, while Assumption (ref).2 is for the purpose of identification. Assumption (ref).3 is somewhat interesting. Take

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

as an example, which says the asymptotic distribution associated with the $i$-th exporter factor at time $t$ is not only driven by the error component, but also is driven by its entire importer network. The same argument applies to the importer factor. In this way, the networks of export and import are entangled with each other. Mathematically, it requires country-specific shocks to have mean 0, which is ensured by (ref). See Assumption 1.ii of CKKK2018 and Assumption 1.a of Han2019 for similar settings.

To close our theoretical investigation, we summarize the asymptotic distributions associated with the global and country-specific factors in the next theorem.

theoremUnder Assumptions (ref)-(ref), Let $(M,N,T) \to (\infty,\infty,\infty)$. \begin{enumerate} • If $\sqrt{MN} ( \frac{1}{T} +\Delta_{g, MNT}^*) \to 0$, then $\normalfont \sqrt{MN} (\widehat{\bm{g}}_t - \bm{g}_t) \to_{D} N (\textbf{0}, \bm{\Sigma}_{\bm{\Gamma}}^{-1} \bm{\Phi}_{t} \bm{\Sigma}_{\bm{\Gamma}}^{-1} )$ for each $t $. \end{enumerate} In addition, let Assumption (ref) also hold. \begin{enumerate} • If $\sqrt{M} ( \frac{1}{\sqrt{T}} + \Delta_{Ij, MNT}^{*}) \to 0$, then $\normalfont \sqrt{M} (\widehat{\bm{f}}_{I,jt} - \bm{f}_{I,jt}) \to_{D} N (\textbf{0}, \bm{\Sigma}_{\bm{\Lambda}_{I,\bullet j}}^{-1} \bm{\Omega}_{I, jt} \bm{\Sigma}_{\bm{\Lambda}_{I,\bullet j}}^{-1} ) $ for each $(j,t)$; • If $\sqrt{N} (\frac{1}{\sqrt{T}} + \Delta_{Ei, MNT}^{*}) \to 0$, then $\normalfont \sqrt{N} (\widehat{\bm{f}}_{E,it} - \bm{f}_{E,it}) \to_{D} N (\textbf{0}, \bm{\Sigma}_{\bm{\Lambda}_{E, i \bullet}}^{-1} \bm{\Omega}_{E, it} \bm{\Sigma}_{\bm{\Lambda}_{E, i \bullet}}^{-1} )$ for each $(i,t)$. \end{enumerate} In the above, $\Delta_{g, MNT}^{*} $, $\Delta_{Ij, MNT}^{*} $ and $\Delta_{Ei, MNT}^{*} $ are defined as follows. \begin{eqnarray*} \Delta_{g, MNT}^{*} &=& \frac{T^{\max_{i}a_{E,i}}}{\sqrt{N}} + \frac{T^{\max_{j}a_{I,j}}}{\sqrt{M}} + \frac{\sqrt{\ln(MN)} \cdot \left( T^{\max_{i}a_{E,i}} + T^{\max_{j}a_{I,j}} \right) }{\sqrt{T}} \\ && + \ln(MN) \cdot ( T^{2\max_{i}a_{E,i}} + T^{2\max_{j}a_{I,j}}); \\ \Delta_{Ij, MNT}^{*} &=& \frac{\ln(MN) \cdot T^{\max_{j}a_{I,j}}}{ \min\{\sqrt{M}, \sqrt{N}, \sqrt{T}\} } + \sqrt{\ln(MN)} \cdot ( T^{\max_{i}a_{E,i}} + T^{\max_i b_{EI, ij}} ) + T^{a_{I,j}}; \\ \Delta_{Ei, MNT}^{*} &=& \frac{\ln(MN) \cdot T^{\max_{i}a_{E,i}}}{ \min\{\sqrt{M}, \sqrt{N}, \sqrt{T}\} } + \sqrt{\ln(MN)} \cdot ( T^{\max_{j}a_{I,j}} + T^{\max_j b_{EI, ij}} ) + T^{a_{E,i}}. \end{eqnarray*}

The condition $\frac{\sqrt{MN}}{T}\to 0$ in the first result of Theorem (ref) is equivalent to $\frac{\sqrt{N}}{T}\to 0$ in Theorem 1 of BN2013 in which a two dimension model is considered. The condition $\sqrt{MN} \cdot \Delta_{g, MNT}^*\to 0$ requires the orthogonality between the global and country-specific factor structures are strong enough in order to achieve the optimal rate $\sqrt{MN}$. If we adopt the orthogonality as in Ando and andreou2019inference, then this condition will completely vanish.

In order to achieve asymptotic normality for the country-specific factors, slightly stronger restrictions (such as $\frac{M}{T}\to 0$ and $\frac{N}{T}\to 0$) are imposed in the body of this theorem on top of Assumption (ref), which is due to the fact that we need to account for the estimation bias caused by Step 1 of the estimation approach. It is noteworthy that $\frac{M}{T}\to 0$ and $\frac{N}{T}\to 0$ imply $\frac{MN}{T^2}\to 0$, which has been discussed above. Therefore, we claim the newly imposed conditions are reasonable, and are only slightly stronger than those used in traditional two dimensional analysis. The conditions $\sqrt{M} \cdot \Delta_{Ij, MNT}^*\to 0$ and $\sqrt{N} \cdot \Delta_{Ei, MNT}^*\to 0$ require the orthogonality between the global and country-specific factor structures are strong enough in order to achieve the optimal rates $\sqrt{M}$ and $\sqrt{N}$. Again, if orthogonality is adopted, these conditions will disappear automatically.

Simulation

In this section, we examine the finite sample performance of the methodology proposed in Section (ref). Specifically, the data generating process (DGP) is as follows.

eqnarray[eqnarray omitted — 140 chars of source]

where $i=1,\ldots,M$, $j=1,\ldots,N$, and $t=1,\ldots,T$. The global factors, country-specific factors and idiosyncratic errors are generated by the following AR(1) processes

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

where $i.i.d.$ stands for independent and identically distributed. The factor loadings are generated as: $\bm{\gamma}_{ij} \sim i.i.d. \ N(\textbf{0}, \textbf{I}_{r_g})$, $\bm{\lambda}_{E,ij} \sim i.i.d. \ N(\textbf{0}, \textbf{I}_{r_{E,i}})$, and $\bm{\lambda}_{I,ij} \sim i.i.d. N(\textbf{0}, \textbf{I}_{r_{I,j}}).$

We consider the following two cases.

enumerate• Let $\phi_g = \phi_{E,i} = \phi_{I,j} = \phi_{u} = 0$, $r_{g} = 3$, $r_{E,i} = 2$ for $i=1,\ldots,M$, and $r_{I,j} = 1$ for $j=1,\ldots,N$; • Let $\phi_g = \phi_{E,i} = \phi_{I,j} = \phi_{u} = 0.5$, and the rest values are the same as those in DGP 1.

For each DGP, we conduct the estimation approach of Section (ref) by letting $M,N,T \in \{20,40,60,80\}$, and implement 1000 replications for each given sample size.

To measure the performance of the proposed estimation approach, we define a few criteria below. First, we measure the detection on different factors, and start from the global factor structure.

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

where $\widehat{r}_g^\ell$ defines the estimated of $r_g$ at the $\ell^{th}$ replication. It is clear that $P_{g,c}$, $P_{g,u}$ and $P_{g,o}$ define the probabilities of correctly, under and over select the number of global factors. For the export factors, we define

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

where $\widehat{r}_{E,i}^\ell $ stands for the estimated $\widehat{r}_{E,i}$ at the $\ell^{th}$ replication. Also, it is obvious that $P_{E,c}$, $P_{E,u}$ and $P_{E,o}$ define the probabilities of correctly, under and over select the number of export factors. Similarly, we can define $P_{I,c}$, $P_{I,u}$ and $P_{I,o}$ for the import factors. The details are omitted for the sake of conciseness.

Second, we measure the estimation on different factors. Recall that we have defined $\textbf{G}$, $\textbf{F}_{E,i}$ and $\textbf{F}_{I,j}$ under (ref), and then further define

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

In the above formulas, we let $\widehat{\textbf{G}}^\ell$ and $\textbf{G}^\ell $ include the estimated and true global factors from the $\ell^{th}$ replication. Similarly, we define $\widehat{\textbf{F}}_{E,i}^\ell$ and $\textbf{F}_{E,i}^\ell$ for the exporter factors, and define $\widehat{\textbf{F}}_{I,j}^\ell$ and $\textbf{F}_{I,j}^\ell$ for the importer factors.

We summarize the simulation results in Table (ref) to Table (ref). Note that due to the limit of space, the results of some combinations of $(M,N,T)$ are dropped in all tables. In Table (ref), it is clear that as the sample size goes up, the values of $P_{g,c}$, $P_{E,c}$ and $P_{I,c}$ converge to 1. When the sample size is relatively small, it seems that we tend to under select the number of factors. Once all $M$, $N$, $T$ are greater than and equal to 40, the selection on the factors is quite accurate. In Table (ref), we consider a DGP with more time series correlation, and the pattern is almost identical to those presented in Table (ref). Table (ref) reports the results of $\mbox{RMSE}_{\textbf{G}} $, $\mbox{RMSE}_{\textbf{E}}$ and $\mbox{RMSE}_{\textbf{I}} $. It is not surprising that all values of $\mbox{RMSE}$ converge to 0, as the sample size goes up. Moreover, the values of $\mbox{RMSE}_{\textbf{E}} $ and $\mbox{RMSE}_{\textbf{I}} $ are larger than $\mbox{RMSE}_{\textbf{G}}$ in general, which should be expected. The reason is that Step 2 includes the estimation bias associated with Step 1, although the bias is negligible in the asymptotic sense under certain restrictions. It is noteworthy that the values of $\mbox{RMSE}_{\textbf{E}}$ are larger than those of $\mbox{RMSE}_{\textbf{I}} $, which is due to the fact that more unobservable factors are included for the exporters.

Having justified the validity of the proposed estimation approach through simulations, we are now ready to move on to the empirical study in the next section.

Empirical Study

In this section, we use the proposed methodology to investigate the international trade flows.

The Data

We use monthly bilateral export volumes of commodity goods among 23 countries/region over the period of 1982-2019. The export flows data are collected from the Direction of Trade Statistics (DOTS) of International Monetary Fund (IMF) available at \url{https://www.imf.org/external/index.htm}. We use the FOB (free on board) value of exports of goods denominated in U.S. dollars and restrict the sample to 506 country-pairs of 23 countries/regions from two major trading groups over a 456-month period from January, 1982 to December, 2019.

itemize• Asia-Pacific Economic Cooperation (APEC): Australia (AUS), China Mainland (CHN), Hong Kong (HKG), Indonesia (IDN), Japan (JPN), Korea (KOR), Malaysia (MYS), New Zealand (NZL), Singapore (SGP), Thailand (THA), Canada (CAN), Mexico (MEX), United States (USA) • European Union (EU): Denmark (DNK), Finland (FIN), France (FRA), Germany (DEU), Ireland (IRL), Italy (ITA), Netherlands (NLD), Spain (ESP), Sweden (SWE), United Kingdom (GBR)

Canada, Mexico and United States are also the members of North American Free Trade Agreement (NAFTA). As they are already included in APEC, we no longer specifically mention NAFTA in this study. It is worth pointing out that a similar dataset is considered in chen2019modeling to investigate the patterns in the dynamic network of international trade. The difference between their study and our paper lies on the setting of factor structure. While we consider multiple layers of the factor structure, their study focuses on one layer only with a different presentation. As a consequence, the two models and the corresponding estimation approaches are not directly comparable.

In what follows, the combination of an export country/region and one of its import partner is referred to as a country pair. For example, the export flow from the United States to Australia and the export flow from Australia to the United States are the bilateral export flows for two different country pairs.

Estimation Results

We first report the estimated numbers of global and country-specific factors. Specifically, only one global factor is identified from the sample. The estimated numbers of exporter factors and importer factors are summarized in Table (ref). As shown in the table, majorities have only 1 or 2 factors with the importer factors of IDN being the only exception.

Figure (ref) shows the estimated global factor which has a clear upward trend. First, let's explain why such a behaviour can be captured under the proposed framework. Note that Assumption (ref) requires $\frac{1}{T}\textbf{G}'\textbf{G}\to_P \bm{\Sigma}_{\textbf{G}}$ only. As a special case, it may possess a form like

eqnarray[eqnarray omitted — 92 chars of source]

where $\tau_t=t/T$, and $g(\cdot)$ can be functions such as $g(w)=w$, $g(w)=w^2$, etc. Therefore, the upward trending is obviously included. Detailed discussions on trending behaviour like (ref) can be seen in YGP2020. As explained in Wang2008 and Jorg2016, the global factor may be interpreted as global shocks on the entire network of international trade, e.g., the Global Financial Crisis. Our finding is somewhat consistent with their arguments. For example, there is a sudden and severe drop around 2009 which captures the so-called “great trade collapse", a consequence of the 2008 financial crisis, occurred between the third quarter of 2008 and the second quarter of 2009. We refer interested readers to BJY2012 for more details on great trade collapse. In addition, we note that the global factor becomes more volatile over the sample period, which may indicate the increasing vulnerability of countries to shocks on trade due to globalization over the past couple of decades.

Figures (ref) - (ref) show the estimated exporter factors and importer factors. Specifically, Figure (ref) and Figure (ref) present he exporter factors associated with the countries of APEC and EU respectively. Figure (ref) and Figure (ref) show the importer factors associated with the countries of APEC and EU respectively. The exporter factors can be interpreted as country-specific shocks of export countries which affect the trade volumes from the exporters to the import partners. Similarly, the importer factors can be interpreted as country-specific shocks of import countries which affect the trade volumes from the importers to the export partners. As mentioned in Section (ref), the exporter and importer factors may capture the unobservable outward and inward multilateral trade resistances (MTRs) for different exporters and importers respectively, which can be seen as measures of outward and inward bilateral trade costs for different exporters and importers. The detailed discussions on the connection between multilateral resistances and country-specific factors can be found in kapetanios2020estimation, where the exporter and importer factors are always referred to as source and destination country factors. For almost all country-specific factors, we can observe the increase of the volatility, especially from the beginning of the 21st century, indicating the increasing instability of the inward and outward bilateral trade costs for most of the countries in our sample. Under the assumption of bilateral trade costs symmetry, it follows that the inward and outward multilateral resistances are the same for the same country anderson2003gravity. By comparing the estimated exporter and importer factors for the same country, it can be seen that this symmetry in the multilateral resistances is partially supported by the data. For example, the exporter and importer factors for USA share the similar trend.

Figure (ref) presents the heat map of the global factor loadings for different country pairs. Since the global factor loadings are positive for all country pairs, we rescale them to $[0,1]$ for better presentation. The global factor loading can be interpreted as the responses of the trade volumes for different country pairs to the global shocks. The colour of each cell reflects the sensitivity of the trade volume between two countries to the global shocks. For example, in Figure (ref), the darkest cell corresponding to the export flow from CAN to USA indicates that the export volume from CAN to USA is the most sensitive relationship among all bilateral export flows in the sample. Also, the country pairs like CHN and HKG, CHN and USA, MEX and USA also show strong sensitivity to the global shocks. The relationship among USA, MEX and CAN partially can be explained by the fact that all three of them are the members of NAFTA, which eliminates some trade barriers among the three parties and promotes the trading activities. The similar patterns can also be observed among countries from EU and Asia respectively. Overall, by comparing the values in different rows and columns of the plot, it can be seen that the trade flows involving USA, CHN and DEU show relatively strong sensitivity to global shocks, which indicates that, in general, they are leading export and import countries worldwide.

Figure (ref) and Figure (ref) present the heat maps of the exporter factor loadings and the importer factor loadings, respectively. Each column in the plots represents a country-specific factor loading corresponding to an exporter or importer factor. Similar to the global factor loading, the exporter and importer factor loadings are rescaled to have values between $-1$ and $1$. The exporter factor loadings corresponding to different export countries measure the responses of their import partners to the shocks on those export countries. The importer factor loadings can be interpreted in the same manner. As shown in Figure (ref), the export flow from CAN to USA is relatively sensitive to the exporter shocks of CAN. This is also the case for the export flow from JPN to USA which is shown to be sensitive to the exporter shocks of JPN. On the other hand, Figure (ref) shows that the export flows from both CAN and JPN to USA are also sensitive to the country-specific importer shocks of USA. The country-specific factor loadings for other countries can be interpreted similarly.

Conclusion

In this study, we specifically consider a three-dimensional panel data model, which has been exposed in the literature but has not been fully solved to the best of the authors' knowledge. On theory, our contributions are the following three-fold: (1). under the scenario that all three dimensions can diverge to infinity, we propose an estimation approach to identify the number of global shocks and country-specific shocks sequentially; (2). the newly proposed approach is easy to implement, and the asymptotic theories are established accordingly; (3). we further conduct intensive numerical studies to examine the finite sample performance of the newly proposed approach using both simulated and real datasets. In the empirical study, we then apply the approach to decompose the network of bilateral trade using country level data from two major trading groups (APEC and EU) over the period 1982-2019. We find that the country-specific shocks become more volatile in recent years, which may indicate the increasing instability of the inward and outward bilateral trade costs over the past couple of decades. In addition, we show that the trade flows involving China mainland, Germany and the United States show relatively strong sensitivity to global shocks, which reflects the fact that, in general, they are leading export and import countries worldwide. We note that the relationship among Canada, Mexico, and the United States is also highly sensitive to different shocks, which somewhat reflects the fact that all three of them are highly economically related through NAFTA that eliminates some trade barriers and promotes the trading activities.

table[table omitted — 3,129 chars of source]
table[table omitted — 3,127 chars of source]
table[table omitted — 2,690 chars of source]
table[table omitted — 938 chars of source]
figure[figure omitted — 320 chars of source]
figure[figure omitted — 352 chars of source]
figure[figure omitted — 341 chars of source]
figure[figure omitted — 352 chars of source]
figure[figure omitted — 341 chars of source]
figure[figure omitted — 255 chars of source]
figure[figure omitted — 263 chars of source]
figure[figure omitted — 263 chars of source]
center[center omitted — 368 chars of source]

This file includes two appendices. Appendix A presents the preliminary lemmas, and the proofs of the main results. We relegate the secondary results and the associated proofs to Appendix B. Specifically, Appendix (ref) scratches the outline of the proofs, while Appendix (ref) presents the preliminary lemmas, which facilitate the development of the main results. Appendix (ref) summaries the proofs for each step. In Appendix B, Appendix (ref) states the secondary lemmas, while Appendix (ref) includes all the corresponding proofs.

\setcounter{page}{1}