Factor Network Autoregressions

Factor Network Autoregressions

\noindentKeywords:\ Networks, factor models, principal components, VAR, tensor decomposition.

\def\spacingset#1 \spacingset{1} \thispagestyle{empty} \spacingset{1.7}

\refstepcounter{lettersection} \oldsection{Introduction}

Network models are key tools for analyzing interconnected systems and are increasingly used in many areas of research. Typical applications involve social networks, where a high number of agents are connected along several dimensions (e.g., zhuetal17), or economic networks, which are usually applied to evaluate spillover effects between agents, geographical regions and economic sectors (e.g., acemogluetal15; dieboldyilmaz14; billioetal12). From an empirical perspective, the analysis of networks calls for the development of novel statistical techniques that are able to account for complex interactions and handle large datasets, possibly recorded over time.

Several approaches have been proposed to integrate network structures into well-established time series models; the network autoregression (NAR) by zhuetal17 or the global vector autoregressive (GVAR) model by pesaranetal04 are two leading examples. These models are generally designed to deal with one type of connection at a time between the agents or “nodes" of a network, e.g., import/export flows. They characterize networks in the form of a single adjacency matrix; i.e., a matrix whose $ij$-th element represents a link between node $i$ and node $j$.

The assumption of a single adjacency matrix, however, appears in general restrictive. In fact, modern economies or social networks are complex systems, in which agents interact through many different channels. For instance, countries are simultaneously linked not only by international trade flows, but also by financial linkages, multinational firms' activities, migration flows, etc. Accordingly, models of multilayer networks, i.e., networks where nodes are linked through multiple types of connections (kivelaet14), have important applications in many different areas, such as international economics and financial (systemic) risk assessment.

To state the problem we consider in this paper, let $y_t$ be an $N$-dimensional vector of stationary time series, e.g., GDP growth rates for different countries. To model the dynamics of $y_t$ when $N$ is large, we could consider a VAR model in which the dynamic relations between the $N$ elements of $y_t$ are mediated by $m$ time-varying networks, denoted as $W_{k,t}$, $k=1,\ldots,m$, e.g., measuring trade flows of different commodities or various financial linkages. All together these networks, which are $N\times N$ matrices, form the $m$ layers of a multilayer network. We could then consider a multilayer NAR:

where $b_k$, $k=1,\ldots, m$, $\varrho$, and $a$ are unknown scalars, and $\zeta_{t}:=(\zeta_{1t},\cdots, \zeta_{Nt})'$ is an $N$-dimensional vector of errors. This model generalizes the NAR by zhuetal17 to the case of multilayer networks. However, in empirical applications, the number $m$ of layers can be very large (e.g., for trade data one may get to hundreds or even thousands of product categories); hence, least-squares estimation of (ref) can become quickly infeasible.

In this paper, we assume that the observed layers $W_{k,t}$, $k=1,\ldots,m$, of the multilayer network have a factor structure, where the common factors are unobserved networks, each represented by an $N\times N$ matrix, capturing the commonality between the layers. We call these latent factors network factors; once retrieved, we can use them in (ref) to model the dynamics of $y_t$ in place of all the $m$ layers, thus reducing the dimensionality of the problem. We call the resulting model a factor network autoregression (FNAR). To estimate the FNAR we first develop a novel, tensor-based, principal component analysis approach, which allows to estimate the network factors. We then use these network factors to estimate an $N$-dimensional FNAR.

Our model exploits two different types of dimension-reduction techniques. First, by assuming, as in zhuetal17, that the network effects are homogeneous across nodes (or across a finite number of groups of nodes), the number of unknown parameters in our FNAR is at most linear in the number of networks considered, instead of being proportional to $N^2$ as is the case for ordinary vector autoregressions. Second, since the FNAR is based on the extraction of just few network factors common across the $m$ network layers, the number of parameters does not grow with the number of layers $m$, but with the number of network factors which we assume to be finite and independent of $N$ and $m$.

In terms of theory, we first provide sufficient conditions such that the FNAR admits a (weakly) stationary and causal solution. Second, we provide sufficient conditions such that our estimators of the network factors and of the FNAR coefficients are consistent and asymptotically normal as the sample size $T$, the number of layers $m$, and the number of nodes $N$, diverge. Crucially, the estimator of the FNAR coefficients we propose is consistent even when the FNAR errors admit cross-sectional correlation as well as serial correlation, due to the presence of node specific factors. Finally, the number of network factors can be determined by an eigenvalue ratio approach, by adapting to the present context those recently proposed by hanetal20 and helitrapani2022.

In an empirical application, we use the FNAR model to study the dynamics of GDP growth rates in a multilayer network of $N=24$ countries with $m=25$ different layers, reflecting international trade flows for different categories of goods and services, and a variety of cross-border financial linkages. In a first step we retrieve 6 network factors common across the $m$ layers, which, once identified, have a clear economic meaning. These are then used to construct the regressors which are employed in a FNAR to study and predict GDP growth rates for all countries. A similar idea was explored by chen2022modeling where, however, the tensor data of trade flows is compressed also along the node (country) dimension, thus making the interpretation of the underlying factors less clear.

The paper is organized as follows. Section (ref) gives a summary of the related literature. Section (ref) illustrates the model. Section (ref) presents the estimation approach. Section (ref) introduces the assumptions and the asymptotic properties of the estimators. Section (ref) evaluates the finite sample performance of our estimators by means of Monte Carlo simulations. Section (ref) presents the empirical application and Section (ref) concludes. All proofs are in the online Supplementary Material, which also contains additional possible estimators, further simulation results and details on the empirical application.

Notation. The order of a tensor is the number of dimensions, also known as modes. The fibers of a given mode are defined by fixing every index but one. We deal only with order-3 tensors. The mode-$q$ matricization, or unfolding, of a tensor $ \mathcal{T}$, denoted as $\text{mat}_{(q)} (\mathcal{T})$ or, equivalently, $ \mathcal{T}_{(q)} $, is a matrix having as columns its mode-$q$ fibers. For a generic order-3 tensor $\mathcal{T}$ of dimensions $d_1 \times d_2 \times d_3$ with generic element $\mathcal T_{ijk}$, the mode-1 matricization is a $d_1 \times (d_2 d_3)$ matrix having as columns the $d_1$-dimensional vectors $\mathcal T_{\cdot jk}$, the mode-2 matricization is a $d_2 \times (d_1 d_3)$ matrix having as columns the $d_2$-dimensional vectors $\mathcal T_{i\cdot k}$, and the mode-3 matricization is a $d_3 \times (d_1 d_2)$ matrix having as columns the $d_3$-dimensional vectors $\mathcal T_{ij\cdot}$.

Finally, we denote the mode-$q$ multiplication of tensor $\mathcal{T}$ by a matrix $X$ of size $p \times d_q$ as $\mathcal{T} \times_q X$, which is a tensor of size $p$ in the $q$-th dimension and the same size as $\mathcal{T}$ in the other dimensions.

\refstepcounter{lettersection} \oldsection{Related literature}

Our modelling approach is especially related to the network autoregression (NAR) model by zhuetal17, the community network autoregression (CNAR) model by chenetal2020cnar and the group network autoregression (GNAR) model by zhu2022simultaneous. Differently from these papers, where just one network is considered, our main contribution is to integrate a large multilayer network into a VAR model by representing multilayer networks as a superposition of common network factors.

Our methodology for extracting the network factors is related to the recent statistical literature on factor analysis of tensor time series building on the tensor Tucker decomposition, which assumes the existence of low-rank tensor factors (chenyangzhang22; hanetal2020iter,hanetal20; helitrapani2022; zhangetal2022; chen2024rank). This approach differs from the canonical polyadic (CP) decomposition which always assumes the factors to be vectors chang2023modelling,hanetal2021.

Differently from all the above cited works, since our aim is to reduce the dimensionality of a multilayer network, our setting is a special case of the Tucker decomposition, where the factors are low-rank only along the layers' mode, while they retain all information along the nodes' two modes, i.e., they are $N\times N$ matrices which we can directly interpret as network factors. Moreover, differently from some of the above cited works, we allow for serial dependence in the idiosyncratic tensor, reflecting the assumption, commonly made in the literature on factor modeling of economic data, that factors account for the cross-sectional (for us, cross-layer) variation of the data (see, e.g., bai2003inferential, in the vector case). In Appendix F.5 we provide a simulation-based comparison between our estimates of the network factor loadings and those obtained by means of the approaches by chenyangzhang22 where no serial idiosyncratic correlation is allowed for.

Our paper is also related to a growing literature that uses tensor decomposition methods to estimate the coefficients of high-dimensional time series models. First, wangetal21 apply tensor decomposition to the order-3 tensor whose slices are the matrices of (unknown) VAR coefficients at different lags, which makes this approach an alternative to ours. However, unlike our approach, wangetal21 do not exploit data on observed networks to estimate the VAR. We refer to Appendix G.2.4 for an empirical comparison between this approach and our FNAR. Second, wangetal24 consider an autoregressive model for tensor-valued time series and use a Tucker decomposition to estimate the autoregressive coefficients, while chang2023modelling model matrix-valued time series using a tensor CP decomposition. Both these last two approaches are more distant from ours, since we do not aim to estimate an AR model for a multilayer network itself, but we aim to use a multilayer network to estimate a vector autoregression in the same spirit as zhuetal17.

Finally, we also contribute to two further strands of literature. First, to the literature on factor and factor-augmented models (e.g., stockandwatson02, baiandng2006) by developing a new framework where factors are matrices rather than vectors, and enter a factor-augmented autoregression by multiplying (weighting) the lagged vector of endogenous variables, rather than being included directly as regressors. Second, to two important streams of the literature on network econometrics, that is: (i) works that investigate the properties of observed networks, such as production networks (see, e.g., acemogluetal12), trading networks in financial and interbank markets (denbeeetal21) and social networks (zhuetal17), and (ii) works concerned with the estimation of network links from the data (e.g., dieboldyilmaz14, billioetal12, barigozzibrownlees19). Our approach combines these two streams of research: on the one hand, we use data on a large number of observed economic networks; on the other hand, we estimate unobserved common network factors driving them.

\refstepcounter{lettersection} \oldsection{Model}

We assume to observe $m$ networks, each represented by a matrix $W_{k,t}$, $k=1,\ldots, m$, of dimension $N\times N$. The networks have no self-loops, so the diagonal elements of the matrices are zero, and are weighted and {directed}, so in general the matrices have real entries and are not symmetric. In the present high-dimensional setting, it is convenient to assume that the weights are normalized in such a way that the elements in each row of $W_{k,t}$ sum to $N$. This can be done without loss of generality (see Section (ref) for specific comments on this aspect).

We then assume that each observed network can be written as a linear combination of $r$ network factors, $F_{k,t}$, $k=1,\ldots, r$, each of them being an $N\times N$ matrix and with $r\ll m$, plus a network $\mathcal E_{t\cdot\cdot k}$ idiosyncratic to the $k$-th layer, which is an $N\times N$ matrix. Specifically, we assume

where $u_{kh}$, $h=1,\ldots, r$, are the scalar factor loadings for network $k$.

If the factors are sufficiently “pervasive”, we can think of replacing the multilayer NAR in (ref) with a FNAR model, which is given by:

where $\beta_k$, $k=1,\ldots, r$, $\rho$, and $\alpha$ are unknown scalars, and $\nu_{t}:=(\nu_{1t},\cdots, \nu_{Nt})'$ is the $N$-dimensional vector of FNAR errors. Notice that the network factors are rescaled by $N$ in agreement with our normalization assumption on the observed networks.

Furthermore, in order to allow also for factors common across the $N$ nodes, we assume each element of the FNAR errors $\nu_{it}$, $i=1,\ldots, N$, to have a factor structure:

where $G_{kt}$, $k=1,\ldots,q$, are node factors, $\lambda_{ik}$, $k=1,\ldots, q$, are the scalar factor loadings for node $i$, and $\epsilon_{it}$ is the node idiosyncratic component.

We call the term $\sum_{k=1}^{r} \beta_k N^{-1} F_{k,t-1}y_{t-1}$ the {\it network effect}, where the coefficients capture the strength of dynamic network effects between nodes, exerted through different types of relationships, which are summarized by means of few network factors. We call the term $\rho y_{t-1}$ the {\it momentum effect}, which captures the direct dynamic interaction of a node with itself. Last, we call the term $\alpha$ the {\it nodal effect}. This can be generalized to include a random effect by adding a term $Z_{t-1}\delta$, where $\delta$ is a $K$-dimensional parameter vector, common to all nodes, and $Z_{t}$ is a $N\times K$ matrix of node-specific exogenous variables. Because of the factor structure in the FNAR errors, the network nodes are correlated with each other not only through the network relationships, but also through the common factors which characterize the cross-sectional dependence at a global level.

By comparing the FNAR in (ref) with the multilayer NAR in (ref) and a standard VAR(1) for $y_t$ we see that, while the latter requires estimating $N^2+N$ parameters, the multilayer NAR requires estimating $m+2$ parameters, and the FNAR further reduces the number of parameters to $r+2$. Hence, the FNAR model provides two forms of dimension reduction, both along the layer direction, by means of the network factors in (ref) and along the node direction by means of the network mediated interactions in (ref).

\refstepcounter{lettersection} \oldsection{Estimation}

Estimation of the FNAR model requires two distinct steps. The first one involves estimation of the $r$ network factors. The second step involves fitting the FNAR equation (ref). Furthermore, it is necessary to determine the number of network and node factors.

\noindentNetwork factors. We start by discussing estimation of the network factors. To this aim, it is convenient to collect the $m$ weight matrices or layers of the multinetwork into a weight tensor of order 3, which we denote as $\mathcal{W}_t$ and has size $N \times N \times m$. In our notation the $m$ network matrices are the frontal slices of the tensor. The factor network structure in (ref) is then rewritten as:

where $\mathcal{F}_t $ is a $N \times N \times r$ tensor containing, as frontal slices, the $r$ network factors $F_{k,t}$, $k=1,\ldots, r$, each of dimensions $N\times N$, which are common across all layers, $U$ is a $m \times r$ matrix of factor loadings, with entries $u_{kh}$, $k=1,\ldots, m$, $h=1,\ldots, r$, determining how each layer of the original network loads on the network factors, and $\mathcal{E}_t$ is an $N \times N \times m$ tensor containing, as frontal slices, the {idiosyncratic networks}, $\mathcal E_{t\cdot\cdot k}$, $k=1,\ldots, m$. The elements of $\mathcal{E}_t$ are allowed to be (i) weakly correlated in all three modes, and (ii) autocorrelated over time (see Section (ref) for the specific assumptions). Finally, from (ref) the mode-3 matricization of $\mathcal W_t$, denoted as $\mathcal{W}_{(3)t}$, is a $m\times N^2$ matrix such that:

where $\mathcal{F}_{(3)t}$ is $r\times N^2$ and $\mathcal{E}_{(3)t}$ is $m\times N^2$. This expression resembles a conventional factor model, with the major difference that each of the $r$ factors is no longer a scalar but a vector of size $N^2$, containing up to $N(N-1)$ non-zero elements (as there are no self-interactions in the network).

To estimate $\mathcal{F}_t$, or, equivalently, $\mathcal{F}_{(3)t}$, first, we compute the sample $m\times m$ outer product of $\mathcal{W}_{(3)t}$:

Second, letting $\widehat{V}^{\mathcal{W}} $ be the $m \times r$ matrix whose $j$-th column is the normalized eigenvector corresponding to the $j$-th largest eigenvalue, $ \widehat{\mu}_j^{\mathcal{W}}$, of $ \widehat{\Gamma}^{\mathcal{W}}$, and $\widehat{M}^{\mathcal{W}}$ be the $r \times r$ diagonal matrix with $ \widehat{\mu}_j^{\mathcal{W}}$ as its $j$-th diagonal entry, we estimate the $m\times r$ loadings matrix $U$ as:

Third, we estimate the mode-3 matricization of the network factors $\mathcal{F}_{(3)t}$ as the principal components (PCs) of $\mathcal{W}_{(3)t}$, i.e., by linear projection of $\mathcal W_{(3)t}$ onto $\widehat U$: $ \widehat{\mathcal{F}}_{(3)t} := (\widehat{U}'\widehat{U})^{-1}\widehat{U}'\mathcal{W}_{(3)t}= N ( \widehat{M}^{\mathcal{W}} )^{-1/2}\widehat{V}^{\mathcal{W}'} \mathcal W_{(3)t}, $ which is an $r\times N^2$ matrix; once folded into an order-3 tensor, it gives the estimated $N\times N\times r$ tensor of network factors:

having as layers $\widehat F_{k,t}:=\text{mat}_{(1)}( \widehat{\mathcal{F}}_{\cdot ,\cdot ,k, t} )$, $k=1,\ldots, r$, which are the $N\times N$ matrices of estimated network factors.

Like in ordinary principal component analysis (PCA), the key intuition is that by exploiting the cross-sectional variation we can estimate the space spanned by the factors. In this case, the relevant cross-sectional dimension is the dimension $m$ of the layers in a network; i.e., the third dimension of the weight tensor $\mathcal{W}_{t}$. The rescaling by $N$ in estimating the loadings reflects the fact that no dimension reduction is applied to the first and second modes of $\mathcal W_t$.

Notice that the estimated loadings and factor tensor are such that they satisfy the identifying conditions: $m^{-1}\widehat{U}'\widehat U=m^{-1}N^{-2} \widehat{M}^{\mathcal{W}}$, which is a diagonal matrix by construction, and

so that the $r$ rows of $\widehat{\mathcal{F}}_{(3)t}$ are the classical normalized PCs of the $m$ rows of ${\mathcal{W}}_{(3)t}$.

There are two main differences between our approach and the one proposed by chenyangzhang22 for estimating common tensor factors from tensor times series admitting a Tucker decomposition. First, we estimate factors using only contemporaneous sample second moments, thus allowing the idiosyncratic tensor to be also autocorrelated. The second difference is that we extract factors along a single dimension of the tensor, namely the dimension of the network layers. This implies that the extracted factors still have a network interpretation; i.e., they are common network factors. Indeed, by construction, the estimated factor matrices $\widehat F_{k,t} $, $k=1,\dots, r$, have zeros along the main diagonal, as the observed network matrices, and have a scale fixed by means of (ref).

To preserve the network properties of the factors, we do not standardize nor demean the elements in $\mathcal{W}_{t} $ along the time dimension. Indeed, our variables are all expressed in the same unit of measurement and standardization would eliminate the fundamental interpretation of $\mathcal{W}_{t}$ as a network, as centering would affect the zero diagonal entries. The same strategy is adopted by chenyangzhang22.

Our estimators in (ref) and (ref) are defined consistently with the assumption that the observed networks have rows summing to $N$. This implies that the estimated network factors in (ref) have variance growing with $N$, as shown in (ref), and they must be rescaled before being used in the FNAR defined in (ref) to ensure that the scale of the estimated network coefficients $\beta_k$ does not depend on $N$. Clearly, we could equivalently work with row-normalized observed networks and then no rescaling by $N$ would be needed anywhere, although this would imply that, as the number of nodes $N$ grows, the entries of $\mathcal W_t$ would have to become smaller and smaller.

\noindentFNAR coefficients. To describe the estimation of the FNAR, it is convenient to introduce some further notation. Hereafter, with reference to (ref), let $y := (y_1,\cdots ,y_T)' = (y_1,\cdots ,y_N)$ be the $T\times N$ matrix of observed data (with $y_t$, $t=1,\ldots, T$, being $N$-dimensional, and $y_i$, $i=1,\ldots, N$, being $T$-dimensional), and define also the $r$-dimensional vector $\beta:=(\beta_1,\cdots,\beta_r)'$ and the $(r+2)$-dimensional vector $\theta := (\beta', \rho, \alpha )'$ of the FNAR coefficients. Moreover, with reference to the node factor structure in (ref), we let $\epsilon:=(\epsilon_1,\cdots ,\epsilon_T)'=(\epsilon_1,\cdots ,\epsilon_N)$ be the $T\times N$ matrix of node idiosyncratic components (with $\epsilon_t$, $t=1,\ldots, T$, being $N$-dimensional, and $\epsilon_i$, $i=1,\ldots, N$, being $T$-dimensional), and we let $G_t:=(G_{1t},\cdots, G_{qt})'$ and $\Lambda_i=(\lambda_{i1},\cdots,\lambda_{iq})'$, so that $G:=(G_1,\cdots ,G_T)'$ is the $T\times q$ matrix of node factors, and $\Lambda:=(\Lambda_1,\cdots ,\Lambda_N)'$ is the $N\times q$ matrix of loadings. Finally, define the order-3 tensor $\mathcal X$ of dimensions $T\times N\times (r+2)$ having in each of the first $r$ layers one of the $r$ matrices $ \mathsf F_k:=\left(N^{-1} F_{k,0} y_{0},\cdots, N^{-1} F_{k,T-1} y_{T-1}\right)^\prime $, $k=1,\ldots, r$, each of size $T\times N$, in layer $(r+1)$ the $T\times N$ matrix $y_{(-1)}:=(y_0,\cdots,y_{T-1})^\prime$, and in layer $(r+2)$ a $T\times N$ matrix of ones, denoted as $1_{T\times N}$. It follows that $\mathcal X_{(1)}:=\text{mat}_{(1)}(\mathcal X)$ is a $T\times N(r+2)$ matrix having as $t$-th row $\text{vec}(X_t)^\prime$ and $\mathcal X_{(2)}:=\text{mat}_{(2)}(\mathcal X)$ is a $N\times T(r+2)$ matrix having as $i$-th row $\text{vec}(X_i)^\prime$, where we define $X_t := (\mathsf F_{1,t-1,\cdot}',\cdots, \mathsf F_{r,t-1,\cdot}', y_{t-1}, \iota_N)$ and $X_i:=(\mathsf F_{1,\cdot,i},\cdots, \mathsf F_{r,\cdot,i}, y_{i}, \iota_T)$, with $\iota_N$ and $\iota_T$ being $N$- and $T$-dimensional vectors of ones, respectively, with $\mathsf F_{k,t-1,\cdot}$ and $\mathsf F_{k,\cdot,i}$ being the $t$th row and $i$th column of $\mathsf F_{k}$, respectively.

According to the above notation, the FNAR in (ref), jointly with the node factor structure in its errors given in (ref), can be equivalently rewritten as:

which, by stacking all $T$ or $N$ equations, can also be written as:

Since both the network factors (contained in $\mathcal X$) and the node factors are unobserved, estimation of (ref) is infeasible. To make estimation feasible, we start by substituting $\mathcal X$ with $\widehat{\mathcal X}$, which, in turn, is obtained by replacing the network factors $F_{k,t}$ with their estimates $\widehat F_{k,t}$, $k=1,\ldots, r$, given by (ref).

An infeasible OLS estimator of $\theta$ would then be obtained by applying the Frisch-Waugh theorem to partial out the effect of either $G$ or $\Lambda$:

where $M_G:= I_T-G(G'G)^{-1}G'$ and $M_\Lambda:=I_N-\Lambda(\Lambda'\Lambda)^{-1}\Lambda'$ are the $T\times T$ and $N\times N$ linear projectors onto the spaces orthogonal to the factors and loadings spaces, respectively. Likewise, given $\theta$, the estimators of $\Lambda$ and $G$ would be the usual PC estimators applied to the FNAR residuals $\widehat{\nu} := y- \text{mat}_{(1)}(\widehat{\mathcal X}\times_3\theta')$.

Following this reasoning, two equivalent estimators of $\theta$ are given by:

where, letting $\widehat \nu^\dag:=y-\text{mat}_{(1)}(\widehat{\mathcal X}\times_3\widehat{\theta}^{\dag'})$ and $\widehat \nu^*:=y-\text{mat}_{(1)}(\widehat{\mathcal X}\times_3\widehat{\theta}^{*'})$, we defined (see also Appendix A)

with $\widehat V^{\widehat\nu^\dag}$ being the $T\times q$ matrix of normalized eigenvectors of $N^{-1}\widehat \nu^{\dag}\widehat \nu^{\dag'}$, and $\widehat M^{\widehat\nu^*}$ being the $q\times q$ diagonal matrix of eigenvalues of $T^{-1}\widehat \nu^{*'}\widehat \nu^{*}$ with corresponding normalized eigenvectors given by the columns of the $N\times q$ matrix $\widehat V^{\widehat\nu^*}$. By iterating between (ref) and (ref), we solve such minimization and compute the two final estimators of $\theta$. As explained at the end of Section (ref), the two are asymptotically equivalent and they might differ numerically just because of the iterative approaches used to compute them.

Finally, notice that, obviously, once we have computed $\widehat{\theta}^\dag$ we can also estimate the loadings $\Lambda$ by linear projection of $\widehat{G}^\dag$ onto $\widehat \nu^{\dag}$, and once we have computed $\widehat{\theta}^*$ we can also estimate the node factors $G$ by linear projection of $\widehat{\Lambda}^*$ onto $\widehat \nu^{*}$. These estimates, however, are not needed for estimating $\theta$.

The outlined estimators are robust to the presence of autocorrelated node factors $G_t$, and, as a consequence, the FNAR errors $\nu_t$ are allowed to be both cross-sectionally and serially correlated. The adopted algorithm is similar to the one proposed by chenetal2020cnar, who in turn adapted the approach by bai2009panel to the NAR setting. Under the assumption of no autocorrelation in the node factors, the OLS and GLS estimators are also valid estimators of $\theta$ (see Appendices B and C, respectively).

\noindentNumber of factors. Letting $\widehat{\mu}_j^{\mathcal W}$, $j=1,\ldots, m$, be the $j$-th largest eigenvalue of $\widehat{\Gamma}^{\mathcal W}$ defined in (ref), we estimate the number of network factors, $r$, by means of the eigenvalue ratio criterion:

where $r_{\max}$ is a predefined maximum number of network factors such that $r_{\max}<\min\{m,T,N^2\}$. This is the criterion proposed by helitrapani2022, which generalizes the approach proposed by hanetal20 to the tensor factor model with autocorrelated idiosyncratic components.

Likewise, letting $\widehat{\mu}_j^{\widehat \nu}$, $j=1,\ldots, N$, be the $j$-th largest eigenvalue of the sample covariance matrix of the FNAR residuals, the number of node factors, $q$, is determined by means of the criterion:

where $q_{\max}$ is a predefined maximum number of node factors such that $q_{\max}<\min\{T,N\}$; see ahn2013eigenvalue. In practice, since the FNAR residuals $\widehat{\nu}_t$ depend on the chosen value of $q$, we can adopt an iterative procedure, which starts by over-estimating $q$ in the first stage, identical to the one described by bai2009panel.

Alternative approaches to estimating $r$ and $q$ are possible; see, e.g., the information criterion by baing02, or the randomized test by trapani2018randomized which is based on the divergence rates of the eigenvalues. Importantly, the latter could also be used to test for no factors.

\refstepcounter{lettersection} \oldsection{Theory}

Assumptions

In the following, we allow for the number of layers $m$ to grow to infinity. Therefore, all our assumptions are stated for the infinite sequence of $N\times N$ networks $W_{i,t}:=\text{mat}_{(1)}( \mathcal{W}_{\cdot ,\cdot ,i, t} )$ with $i\in\mathbb N$. Equivalently, we could state the assumptions for $i=1,\ldots, m$ with $m\in\mathbb N$. Moreover, all assumptions are stated contemplating the possibility that also the number of nodes $N$ and the sample size $T$ grow to infinity.

Assumptions (ref)(ref) and (ref)(ref) imply that we consider only pervasive, or strong, factors. In other words, the network factors are loaded by most or all the network layers. Notice that the factor tensor has dimension $N\times N\times r$, hence the rescaling introduced in Assumption (ref)(ref), which accounts for its first two modes having diverging dimensions. Furthermore, under Assumption (ref)(ref), the 2nd order moments of the process $\{N^{-1}\text{vec}(\mathcal F_{(3)t}), t\in\mathbb Z\}$ are finite and independent of time, for all $N\in\mathbb N$.

Assumptions (ref)(ref) and (ref)(ref) imply that, given the factors, we can consistently estimate $\Gamma^{\mathcal F} := \mathbb{E} [ \mathcal{F}_{(3)t} \mathcal{F}_{(3)t}' ]$, as proved in Lemma E.5(i). Therefore, given the loadings, we can also consistently estimate $\Gamma^\chi:= U\mathbb E [ \mathcal{F}_{(3)t} \mathcal{F}_{(3)t}' ] U' $. Assumption (ref)(ref) simply states that the number of factors is finite for all $m,N\in\mathbb N$ and that in order to find such factors we need $m$ to be large enough. Finally, Assumption (ref)(ref) is very mild as it simply requires each element of the network factors to have finite first moment.

This assumption controls the serial and cross-sectional dependence of the entries of the idiosyncratic tensor. In particular, according to Assumption (ref)(ref), the covariance across time and layers is controlled in a standard way; hence, the classical summability conditions hold (see Lemma E.1 and bai2003inferential). The covariances across the elements of the $N^2$-dimensional vector $\mathcal E_{(3)t j\cdot}$ are of order $N^4$ for any given $t=1,\ldots, T$ and $j=1,\ldots, m$, and we require to scale their sum by $N^2$, thus assuming a standard summability condition. Assumptions (ref)(ref), (ref)(ref), and (ref)(ref) imply that, given the idiosyncratic tensor, we can consistently estimate $m^{-1}N^{-2}\Gamma^{\mathcal E}$, for any $m,N\in\mathbb N$, as proved in Lemma E.5(ii).

Uncorrelatedness of the processes $\{\mathcal F_{(3)t}\}$ and $\{\mathcal E_{(3)t}\}$ is a natural assumption, while the moment conditions controlling higher-order dependence are standard and are direct extensions of what typically assumed in the vector case bai2003inferential.

Define the $m\times m$ matrix $ \Gamma^{\mathcal W}:=\mathbb E[\mathcal W_{(3)t}\mathcal W_{(3)t}']$ and let $\Gamma^\chi$ be as previously defined. Denote the eigenvalues of $\Gamma^\chi$ as $\mu_j^\chi$, $j=1,\ldots, r$, in decreasing order. Then, as established by Lemmas E.1(i) and E.1(iv), we have that $\mu_j^\chi\asymp mN^2$ and $\norm{N^{-2}\Gamma^{\mathcal{E}}}$ is finite for all $N\in\mathbb N$. As a consequence, by Assumption (ref) and Weyl's inequality, the matrix $\Gamma^{\mathcal W}=\Gamma^\chi+\Gamma^{\mathcal E}$ is characterized by an eigengap between the $r$-th and the $(r+1)$-th largest eigenvalues which widens as $m\to\infty$. This property allows us to identify the number of network factors $r$ and it is the rationale for the eigenvalue ratio criterion for estimating $r$ defined in (ref). This also means that the network factor model in (ref) is always identified as long as $m\to\infty$.

In general, the network factors and their loadings are not separately identified unless we impose further restrictions. To this end we make the following assumption.

Under Assumption (ref) the columns of the loadings matrix $U$ and the layers of the tensor factor $\mathcal F_t$ are identified up to a sign multiplication. This identification scheme is a classical one adopted for example by bai2009panel in the vector factor model case.

Assumption (ref)(ref) is standard in the vector factor model case; i.e., when $N=1$, where it is satisfied for example by strong-mixing processes bai2003inferential. In Assumption (ref)(ref) we directly assume a cross-sectional CLT, which is a standard approach in the vector factor model case bai2003inferential.

The FNAR errors $\nu_t$ follow a factor model given in (ref), characterized by the following assumption.

This assumption is similar to the usual set of assumptions for the vector factor model bai2003inferential. The comments to these assumptions are analogous to those previously made for the network factors and are therefore omitted. Concerning the idiosyncratic components, we follow chenetal2020cnar and assume zero correlations both in time and across units. Given that we are considering a factor structure for the FNAR errors, this assumption is not very restrictive, as most of the correlations are likely to be already captured by the lagged terms in the FNAR and by the common factors $G_t$. Nevertheless, it is possible to develop the following asymptotic theory by allowing for the usual kind of weak cross- and autocorrelations between the components of $\epsilon_t$.

Assumption (ref) is taken from bai2009panel and is made just to simplify the proof. In principle it could be relaxed to allow for weak dependence between $\{G_t\}$ and $\{\epsilon_{t}\}$ by means of a condition similar to those required in Assumption (ref), which, in turn, derives from bai2003inferential.

Because of Assumptions (ref)(ref), (ref)(ref), (ref)(ref), (ref)(ref), and (ref), the FNAR errors have covariance matrix $V= \Lambda\Gamma^G\Lambda' + S$, with $S=\mathbb{E} \left[ \epsilon_{t} \epsilon_{t}' \right]$, which is positive definite for all $N\in\mathbb N$. This also implies that $V^{-1}$ is finite for all $N\in\mathbb N$. Moreover, $V$ has the usual eigengap property; i.e., its largest $q$ eigenvalues diverge at rate $N$, while the remaining $N-q$ stay bounded for all $N\in\mathbb N$. This implies that the factor model in (ref), and therefore the number of factors $q$, is always identified as $N\to\infty$. This is the rationale for the eigenvalue ratio criterion for estimating $q$ defined in (ref).

Under Assumption (ref) the columns of the loadings matrix $\Lambda$ and the factors $G_t$ are identified up to a sign multiplication.

Turning to the FNAR defined in (ref), we make the following assumption.

This assumption is a generalization to the case of random multivariate AR models of the usual stability conditions for a VAR. As shown below it implies, together with Assumptions (ref)(ref) and (ref), that the FNAR has a stationary solution for all $N\in\mathbb N$. Notice that Assumption (ref)(ref) is stated for the general case in which $\mathbb E[F_{j,t}]\ne 0$, otherwise the condition needed to ensure the existence of the mean is simply $|\rho|<1$.

To get an intuition of this assumption, consider the $m\times (r+2)$ matrix process $\{\mathcal E_{(3)t} (y_{t-1}\otimes X_t)\}$. We are saying that this process is weakly correlated along the time dimension, which is a standard requirement, but it is also weakly correlated across its $m$ rows. The latter requirement is fulfilled by the idiosyncratic terms $\mathcal E_{(3)t}$ via Assumption (ref)(ref), and here is extended to the case in which $\mathcal E_{(3)t}$ is multiplied by $y_{t-1}\otimes X_t$ which is weakly dependent of $\mathcal E_{(3)t}$ because of Assumption (ref).

Notice that we do not rule out autocorrelation in the node factors, as only the node specific idiosyncratic components are required to have no correlation for the above CLTs to hold, and this is ensured by Assumption (ref)(ref). These assumptions are similar to the conditions in bai2009panel.

Stationarity

In order to develop the theory for the FNAR, we first discuss under which conditions equation (ref) admits a stationary causal solution. Given the difficulty of the problem we limit ourselves to consider weakly stationary solutions. This poses two issues. First, the FNAR is defined for an $N$-dimensional vector $y_t$ where we allow $N\to\infty$. Second, the FNAR is an autoregressive model with stochastic time-varying coefficients. Regarding the former issue, we adopt the definition proposed by zhuetal17.

Turning to the second problem, we have the following result.

In general, one might object that when $N\to\infty$ a meaningful concept of stationarity cannot be stated, as no causal solution to the FNAR can exist since Assumption (ref) will break down, see, e.g., the remark by zhou2020network, in a similar context. What we mean by Definition (ref) and Proposition (ref) is that we are implicitly assuming that there exists a space where the causal solution is well-defined even when $N \to\infty$. This is the same approach adopted by zhuetal17. An interesting implication of this definition is that, under our assumptions, we can ensure that any finite linear combination of the elements of $\{y_t\}$ satisfies a finite dimensional FNAR with a causal solution.

Asymptotic properties of network factors and FNAR coefficients

Consistency and asymptotic normality of the estimated network factor loadings are given next.

Theorem (ref) shows that, when applying PCA to a given mode of the tensor $\mathcal W_t$, the dimensions of all other modes contribute to a faster convergence rate, hence allowing for more degrees of freedom. This is an advantage with respect to the vector case, since even for moderately small values of $T$ we can still have good estimates of the loadings matrix and therefore of the network factors. In particular, we see that the estimated loadings vector $\widehat u_i$ has a consistency rate $\min(m,N\sqrt T)$ and is asymptotically normal if $N\sqrt T/m\to 0$. This is the generalization to the multilayer network case (i.e., to order-3 tensors) of the usual vector case, which corresponds to setting $N=1$ (see bai2003inferential).

Next we prove consistency and asymptotic normality of the estimated network factors.

Theorem (ref)(ref) proves consistency of the whole network factor tensor. Theorem (ref)(ref) proves asymptotic normality of any given column of $\widehat {\mathcal F}_{(3)t}$, which is equivalent to asymptotic normality of any of the $N^2$ entries of each of the $r$ layers of the multilayer network factor $\widehat {\mathcal F}_t$. This is the natural generalization to the multinetwork case of the usual vector case, i.e., when $N=1$ (see bai2003inferential).

We then turn to the asymptotic properties of the estimated FNAR coefficients. Hereafter, let

with $J$ as in Theorem (ref)(ref). Then, we analyze the properties of the estimators $\widehat{\theta}^\dag$ and $\widehat{\theta}^*$ by noticing that

where $u_i:= (X_i\bar J-\widehat X_i)\bar J\theta$ and $u_t:= (X_t\bar J-\widehat X_t)\bar J\theta$, and recall that $\theta := (\beta', \rho, \alpha )'$.

The following theorem holds.

Parts (ref) and (ref) extend Theorem 2 in bai2009panel to the FNAR case. The interesting cases are parts (ref) and (ref), where we do not impose homoskedastic idiosyncratic components in the FNAR errors. Notice that the network coefficients, $\beta_j$, $j=1,\ldots, r$, which are the first $r$ elements of $\theta$, are consistently estimated only up to a sign, due to the indeterminacy in the identification of the network factors.

Estimators of the asymptotic variance-covariance matrix of $\widehat \theta^\dag$ and $\widehat\theta^*$ under the assumptions in parts (ref) and (ref) are (see also bai2009panel):

where $ \widehat Z_{i} := M_{\widehat G^\dag}\widehat{X}_i-\frac 1N\sum_{k=1}^N \left(\widehat\Lambda_i^{\dag'}\left(\frac{\widehat{\Lambda}^{\dag'}\widehat\Lambda^\dag} N\right)^{-1}\widehat\Lambda_k^\dag\right)M_{\widehat G^\dag}\widehat X_k , $ and $ \widehat W_t:= M_{\widehat \Lambda^*}\widehat X_t-\frac 1T\sum_{s=1}^T (\widehat G_t^{*'}\widehat G_s^* )M_{\widehat \Lambda^*}\widehat X_s . $

Four important comments about this result follow. First, the proof is based on showing that if ${ \sqrt {NT}}/ m \to 0$ and $N/m\to 0$, as $m,N,T\to\infty$, then the network factors can be treated as observed; i.e., the generated regressors bias is asymptotically negligible (see Proposition C.1). Under these conditions, and if also $T/N\to 0$ and $\sqrt N/T\to 0$, the iterative estimators are $\sqrt{NT}$-consistent.

Second, if the node factors are not autocorrelated we can also compare the iterative estimators with the OLS and the GLS estimators studied in Appendix B and C, respectively. In this case, the GLS is also $\sqrt{NT}$-consistent since, similarly to the iterative estimator, it rescales $X_t$ by the FNAR error covariance matrix $V$, which is $O(N)$ by Assumption (ref). For the same reason the OLS estimator is just $\sqrt T$-consistent since it does not control for the FNAR error covariance, and it would be $\sqrt{NT}$-consistent only if $V$ were a diagonal matrix; i.e., when no node factor is present, as assumed by zhuetal17. The same results on OLS and GLS are obtained by chenetal2020cnar for the case of observed networks.

Third, the two estimators are asymptotically equivalent and thus also equally efficient. To see this notice that $\widehat{\theta}^\dag$ and $\widehat{\theta}^*$ are such that they solve

so the two losses are identical and must have the same minimum. Once we fix the identification constraints as in Assumption (ref), the only difference is then about the implementation of the minimizations. For $\widehat{\theta}^{\dag}$, by replacing $\Lambda_i= (G^\prime G)^{-1} G^\prime(y_i-X_i\theta)=T^{-1}G^\prime(y_i-X_i\theta)$ in (ref), we can solve for $G$ and $\theta$ only. For $\widehat{\theta}^{*}$, by replacing $G_t=(\Lambda^\prime\Lambda)^{-1}\Lambda'(y_t-X_t\theta)$ in (ref), we can solve for $\Lambda$ and $\theta$ only. These two procedures lead to the two solutions given in (ref). However, since in practice the solutions are obtained by iteration, the two estimates might not coincide exactly, although they will be very similar. And, as expected, the estimated standard errors will coincide (see the results in Section (ref) and Appendix F.4).

Fourth, and last, we should view Theorem (ref) as giving the asymptotic distribution of the theoretical estimator minimizing (ref) or (ref). This is the same point of view adopted by bai2009panel. In practice, it might be important to investigate how the initialization of the algorithm affects such convergence. In the simpler case of a panel regression having errors with a factor structure, jiang2021recursive show that any initial estimator could still lead to a consistent iterated estimator, depending on the structure of the regressors which can be quite general. However, the estimator computed in practice might have a slower convergence rate if the initial estimator is not consistent. We do not explore this aspect further here, but we limit ourselves to notice that in our numerical exercises of Sections (ref) and (ref), convergence is always achieved in few steps and the iterated estimator works well even in presence of weak serial correlation of the node factors.

\refstepcounter{lettersection} \oldsection{Monte Carlo simulations}

To evaluate the finite sample performance of the proposed estimators, we generate artificial time series of $y_t$ and $\mathcal{W}_t$, for $t=1, \dots, T$, according to the model equations (ref), (ref), and (ref). We fix $r=1$ and $q=1$; i.e., one network factor and one node-specific factor. We also fix the values of FNAR parameters $\beta=0.5$, $\rho=0.3$, and $\alpha=0.2$. We consider $N \in \{10,20,50,100,200\}$ nodes, $m \in \{20,50,100\}$ layers, and $T \in \{50,100\}$ time periods. Also, for each value of the pair ($N ,m$), we randomly generate the entries of the loading vectors $U$ and $\Lambda$ once (and independently) from $\mathcal{N}(1,1)$, and then we keep them fixed across Monte Carlo (MC) iterations (see chenetal2020cnar). All other quantities are generated at each MC iteration. All results are based on 500 iterations. Full details on the data generating process are in Appendix F.1.

Tables (ref)-(ref) report the {RMSE} and Relative RMSE ({ReRMSE}) of the estimates. Here we report results under case II, which corresponds to idiosyncratic terms $\mathcal E_t$ having serial and cross-layer correlation, and using the iterative estimator $\widehat{\theta}^\dag$ in (ref). Additional results for case I of uncorrelated idiosyncratic terms are in Appendix F.2. As predicted by the theory, the accuracy of estimates for $\beta$, $\rho$ and $\alpha$ improves with both $N$ and $T$, and the RMSE of network factors and loadings decreases when the number of layers $m$ increases. Furthermore, as shown in Appendix F.3, the MC distributions of the estimated network coefficient are all strongly centered around the true value $\beta=0.5$ and become narrower as $T$ and $N$ increase.

Last, in Appendix F.5 we compare our estimates of the loadings $U$ with those obtained using the TOPUP and TIPUP estimation methods proposed by chenyangzhang22. As expected our approach improves over those estimators in presence of serial idiosyncratic correlation.

\refstepcounter{lettersection} \oldsection{Empirical application}

In this section, we present an application of the FNAR for studying cross-country macroeconomic interdependence determined by global trade flows and cross-border financial relationships.

\noindentData. For a sample of $N=24$ countries, we use $m=25$ networks constructed using bilateral import/export flows for different good (layers 1--9) and services (10--19) categories, bilateral financial positions for different types of financial claims (20--23) and cross-border mergers and acquisitions classified by sector of economic activity (24--25). The list of countries and network layers, including details on how the networks are built, are given in Appendix G.1.

\noindentNetwork factors. Due to data limitations in the time series of financial positions, we collect data for the networks at the annual frequency from 2001 to 2019, so the factor analysis is conducted on a sample of length $T_1=19$.\footnote{We have few missing values over this period. In these cases, we use the previous year's value or the closest available year's value.} Although this is a short time span, we recall that in tensor PCA the effective sample size when estimating the loadings space is $N^{2}T_1$ (see Theorem (ref)).

We then extract the common network factors from the 25 layers of the network, and we set $\widehat r=6$ network factors, as in chenyangzhang22. From Figures G.9-G.12 in Appendix G.2 we can interpret the six network factors as follows. The first network factor conveys approximately the average country weights across all layers of the network. In particular, the factor values are very close to the average weights (scaled by a constant), and the loading coefficients are almost the same for all layers. The countries with the largest factor weights for the US are its major economic partners: Canada, UK, Mexico, China, Japan and Germany.

The second factor captures a difference between financial relationships and trade in goods. The factor loadings for financial layers have opposite sign (positive) compared to the loadings for trade-in-goods layers (negative). Recall that factors are identified up to a sign. The largest positive weights are assigned to economies having relatively large financial sectors with global reach: UK, US, and Hong Kong. In the case of the US connections, a large positive weight is assigned to the UK, whose tight economic links with the US are mostly concentrated in the financial sector, and large negative weights are assigned to Canada, Mexico, and China, i.e., the US biggest trade partners.

The third factor distinguishes between equity and debt relationships, being the only factor where equity, on the one hand, and debt, on the other hand, show loadings with opposite signs. The fourth factor is strongly associated with M&A relationships. The fifth factor is mainly driven by agricultural/extractive goods (positive weights, especially for vegetable fuels, oils, fats, and waxes). It also loads on trade in manufacturing goods (negative weights). Positive weights are assigned to countries with strong trade links with the US in non-manufacturing sectors, such as Canada, Saudi Arabia, and Italy, while negative weights are associated with large manufacturing partners, like China. This factor also distinguishes between stocks of portfolio holdings and flows associated with M&A deals and banking. Finally, the sixth factor captures a distinction between goods-sector M&A integration and services-sector integration.

Next, in line with conventional PCA, we evaluate the fraction of variance in network layers explained by each factor, denoted as $v^{(k)}$, $k=1,\ldots, 6$ (computed as in Appendix G.2). We have $v^{(1)}=0.68$, $v^{(2)}=0.07$, $v^{(3)}=0.03$, $v^{(4)}=0.03$, $v^{(5)}=0.02$, and $v^{(6)}=0.02$. Thus, overall the 6 factors explain about 85% of the total variance of $\mathcal W$. However, the importance of different factors varies greatly across countries; see Table G.11 in Appendix G.2.

\noindentFNAR coefficients. The endogenous vector $y_t$, $t=1,\ldots, T$, collects (quarterly) real GDP growth rates for all $N$ considered countries and for the sample 2001Q1-2019Q4; i.e., $T_2=76$. To address heterogeneity of nodal and momentum effects, we split the countries into two groups: (1) advanced economies ($N_1=15$), and (2) emerging economies ($N_2=9$) and the vector $y_t$ is partitioned accordingly as $y_t=(y_t^{(1)'}; y_t^{(2)'})'$. We consider the following FNAR, for $t=1,\ldots, T$,

where $\widetilde F_{j,t} = \widehat F_{j,\tau}$ for $4(\tau-1)+1\le t \le 4\tau$, $\tau=1,\ldots, T_1$. In other words, the network factors $F_{k,t}$, $t=1,\ldots, T_1$, which are computed on a yearly basis, are treated as constant throughout all quarters of a given year. Hereafter, we let $\theta:=(\beta', \rho^{(1)},\rho^{(2)},\alpha^{(1)}, \alpha^{(2)})'$.

By means of the criterion defined in (ref), we find evidence of one common node factor, i.e., $\widehat q=1$. We then estimate the model by GLS as described in Appendix C. Last, we consider the iterative estimators $\widehat\theta^\dag$ or $\widehat\theta^*$ defined in (ref) and we initialize the algorithm by using the GLS estimator and the estimated node loadings, $\widehat \Lambda$, and factor, $\widehat G_t$, computed by PCA on the GLS residuals as described in Appendix A. Since these residuals do not display significant autocorrelation, we are confident that the GLS estimator is $\sqrt{NT}$-consistent and, based on the results of jiang2021recursive we conjecture that Theorem (ref) holds for our iterated estimators. Convergence is reached in 8 or 4 iterations for $\widehat\theta^\dag$ or $\widehat\theta^*$, respectively.

Table (ref) reports the estimated coefficients and their standard errors with significance reported according to the usual $Z$-test. The coefficients on $N^{-1}\widehat F_{1,t-1}y_{t-1}$ and $N^{-1}\widehat F_{5,t-1}y_{t-1}$ are always strongly significant, while there is mixed evidence regarding the coefficients on $N^{-1}\widehat F_{2,t-1}y_{t-1}$, $N^{-1}\widehat F_{4,t-1}y_{t-1}$, and $N^{-1}\widehat F_{6,t-1}y_{t-1}$ which are mildly significant and not for all estimates.

Based on the interpretation of the first network factor, the coefficient on $N^{-1}\widehat F_{1,t-1}y_{t-1}$ captures a general network effect operating through aggregate economic linkages. Given the loadings of factor 5 in Figure G.10, the coefficient on $N^{-1}\widehat F_{5,t-1}y_{t-1}$ indicates that trade in mineral fuels and in animal and vegetable oils (major inputs of chemical industry) has the main impact on GDP growth. Apart from this, trade linkages tend to generate larger spillovers in manufacturing sectors (layers 6-9) than in non-manufacturing sectors (layers 1-3), and financial linkages tend to generate larger spillovers when they take the form of M&A or flows of banking assets (rather than portfolio holdings).

Finally, based on these estimates, we can approximate the network effects associated with the original layers of the network, by appropriately rescaling the estimated network effect coefficients $\beta$. Specifically, given the definition of estimated loadings in (ref) and the properties of tensor multiplication, and letting $\widehat{\mathcal{W}}_{t} := \widehat{\mathcal{F}}_{t} \times_3 \widehat{U}$, for a given estimate $\widehat{\beta}$ we have that:

Thus, $N^2 \widehat{U} (\widehat{M}^{\mathcal{W}})^{-1}\widehat{\beta}$ is the vector of network effects in terms of the row-normalized tensor $\widehat{\mathcal{W}}_{t-1} / N$ and its entries are shown in Figure (ref), when computed using the iterated estimator $\widehat{\theta}^*$. The figure confirms a substantial heterogeneity of effects across layers, reflecting their different loadings on the network factors.

\noindentForecasting. We conclude by studying the performance of our FNAR when producing 1-quarter-ahead forecasts of GDP growth rates based on a recursive window exercise from 2002Q1 to 2019Q4. We consider the following competitors: the tensor-based estimators MLR and SHORR by wangetal21; a NAR estimated either using the $m$ layers of the common component (denoted as TUCKER COMMON and further regularized via LASSO or Ridge due to collinearity of the regressors) or just $r$ network factors (denoted as TUCKER FACTORS) both obtained from a full Tucker factor decomposition, estimated via TOPUP as in chenyangzhang22 (see Appendix H for more details); a multilayer NAR estimated via LASSO or Ridge; and an ordinary VAR. Details on the adopted forecasting scheme and on the implementation of the alternative estimation methods are in Appendix G.2.4.

Table (ref) reports the root mean squared forecast errors (RMSFE) in terms of percentage points of GDP growth, for each country considered. In the last two rows of Table (ref), we report the average (across all countries) RMFSE and relative RMSFE (ReRMSFE) with respect to the FNAR, for all forecasting methods (values larger than one indicate a better performance of the FNAR). For most countries, and on average, our approach delivers the forecasts with smallest RMSFEs. In particular, we outperform the approaches based on a full Tucker decomposition which are our most natural competitors. Indeed, contrary to the case of factors extracted by means of a full Tucker decomposition, our network factors still contain terms which are idiosyncratic to the network nodes, i.e., countries, which are potentially relevant for predicting country specific GDP growth rates.

\refstepcounter{lettersection} \oldsection{Conclusions} In this paper, we have introduced a factor network autoregression (FNAR) for time series characterized by multiple network effects. Estimation is based on two steps. First, we extract few network factors common across the layers of the underlying multilayer network. Second, we estimate a factor-augmented NAR or FNAR where the network effects are determined by the latent network factors. The FNAR errors are allowed to have an underlying factor structure capturing common correlations across the network nodes. We prove consistency and asymptotic normality of the proposed estimators as the number of layers, nodes and time observations diverges to infinity.

The results of an empirical application show that, by accounting for cross-country economic and financial linkages, the model provides a rich description of the dynamics of GDP growth rates and produces accurate forecasts.

We outline three possible extensions of this work, which we leave for further research. First, by adapting the works by chenetal2020cnar and zhu2022simultaneous to the FNAR framework, we could consider a FNAR with momentum and nodal coefficients which are group specific, where the group structure is unknown and the number of groups $K$ can grow with the number of nodes $N$. Second, by generalizing to the tensor setting the three-pass regression filter by kelly2015three, we could improve the performance of our estimator by accounting also for the information contained in the vector of dependent variables $y_t$ when extracting the network factors. Third, by extending to the tensor case an approach similar to the one proposed by wu2020adaptive, we could allow for time-varying network factor loadings under the standard assumption of local stationarity.

\newcounter{lettersection} \setcounter{lettersection}{0}

\let\oldsection\section

\setcounter{figure}{0} \setcounter{table}{0}

This supplemental material contains several appendices to our paper. In Appendix (ref) we describe how to estimate the node factors by Principal Components Analysis (PCA) on the FNAR residuals. In Appendices (ref) and (ref), we introduce the OLS estimator and a GLS-type estimator of the FNAR coefficients. In Appendix (ref), we provide the proofs for Propositions and Theorems. In Appendix (ref), we provide and prove all auxiliary results. In Appendix (ref) we provide additional simulation results. In Appendix (ref) we provide additional information on the data used in our empirical application, as well as additional empirical results. Finally, in Appendix (ref) we consider alternative approaches for the estimation of a multilayer NAR, based on full Tucker decompositions.

\refstepcounter{lettersection} \oldsection{Estimation of node factors}

Once we compute the OLS estimator, let $\widehat\nu:=(\widehat\nu_1\cdots \widehat\nu_T)'$ be the $T\times N$ matrix of residuals of the FNAR such that $\widehat \nu_t:=y_t -\widehat X_t \widehat{\theta}^{\text{\upshape\tiny OLS}}$, $t=1,\ldots, T$. Then, the node factors $G_t$ and their loadings $\Lambda$ can be estimated by PCA in two equivalent ways. Specifically, consider either the $N\times N$ or $T\times T$ sample covariance matrices

Then, letting $\widehat G:=(\widehat G_1,\cdots, \widehat G_T)'$ and $\widetilde G:=(\widetilde G_1,\cdots, \widetilde G_T)'$ be the estimated $T\times q$ matrices of factors, the PC estimators are given by:

where $\widehat M^{\widehat\nu}$ is the $q\times q$ diagonal matrix of eigenvalues of $\widehat {\Gamma}^{\widehat \nu}$ with corresponding normalized eigenvectors being the columns of the $N\times q$ matrix $\widehat V^{\widehat\nu}$, and $\widetilde V^{\widehat\nu}$ is the $T\times q$ matrix of normalized eigenvectors of $\widetilde {\Gamma}^{\widehat \nu}$. It is easy to verify that, regardless of the choice made in (ref), $\widehat G\widehat \Lambda'=\widetilde G\widetilde \Lambda'$.

Notice that the estimated loadings and factors are such that $\widehat \Lambda'\widehat \Lambda /N$ and $\widetilde \Lambda'\widetilde \Lambda /N$ are diagonal and \[ \frac 1T\sum_{t=1}^T \widehat G_t\widehat G_t'= I_q,\quad \frac 1T\sum_{t=1}^T \widetilde G_t\widetilde G_t'= I_q \] so that the estimated factors are orthonormal.

\refstepcounter{lettersection} \oldsection{OLS estimator and its asymptotic properties}

The OLS estimator of $\theta$ is given by:

This is the estimator proposed by zhuetal17 and chenetal2020cnar for the case in which the network is observed and, respectively, when the FNAR errors are either uncorrelated or have a factor structure with serially uncorrelated node factors.

We start by making the following standard assumption.

Because of Assumptions (ref)(ref), (ref)(ref), and (ref), we have that the FNAR errors $\{\nu_t\}$ are not autocorrelated. This is necessary for the CLT in the next part of this assumption to hold. Assumption (ref)(ref) is also found in baiandng2006. In fact, Assumptions (ref)(ref) and (ref)(ref) are made for simplicity and could be proved in a similar way as in chenetal2020cnar.

The OLS estimator in (ref) satisfies:

The asymptotic properties of the terms of (ref) are given in the following Proposition.

The next theorem follows.

From Proposition (ref), we see that if $\sqrt{ T}/ m \to 0$, as $m,N,T\to\infty$, then the network factors can be treated as observed and Theorem (ref) follows. In particular, by virtue of Assumptions (ref)(ref) and (ref)(ref), the OLS estimator is $\sqrt T$-consistent and asymptotically normal. Notice that the requirement $\sqrt T/m\to 0$ for Theorem (ref) to hold is analogous to the one assumed in the vector case by baiandng2006.

An estimator of the asymptotic variance-covariance matrix of the OLS estimator is then given by (see baiandng2006): \[ \widehat{\operatorname*{Avar}}\left[ \sqrt{T} (\widehat{\theta}^{\text{\upshape\tiny OLS}} -\bar J\theta) \right] = \left(\frac{1}{NT} \sum_{t=1}^{T} \widehat{X}_{t}^' \widehat{X}_t \right)^{-1} \left(\frac{1}{N^2T} \sum_{t=1}^{T} \widehat{X}_t^' \widehat \nu_t \widehat \nu_t' \widehat{X}_t \right) \left( \frac{1}{NT} \sum_{t=1}^{T} \widehat{X}_{t}^' \widehat{X}_t \right)^{-1}. \]

\refstepcounter{lettersection} \oldsection{GLS estimator and its asymptotic properties}

A GLS extension of the OLS estimator of the FNAR considered in Appendix (ref) can be computed by means of the following procedure, initially proposed by chenetal2020cnar for the special case where the network is observed.

Once we estimate the node factors and their loadings as described in Appendix (ref), let $\widehat{\epsilon}:=\widehat \nu-\widehat G\widehat \Lambda' = \widehat \nu-\widetilde G\widetilde \Lambda'$ and $\widehat S$ be the diagonal matrix with entries the diagonal entries of $T^{-1}\widehat{\epsilon}'\widehat{\epsilon}$. We can estimate the covariance matrix $V$ of $\nu_t$ as $\widehat V :=\widehat \Lambda\widehat \Lambda'+\widehat S$ and, by applying the Sherman-Morrison-Woodbury formula, its inverse as:

where $\widehat S^{-1}$ is a diagonal matrix and hence easy to compute.

The GLS estimator of $\theta$ is then given by:

To study the asymptotic properties of the GLS estimator, we make two more assumptions. First, we extend Assumptions (ref)(ref) and (ref)(ref) to the following.

Because of Assumptions (ref)(ref), (ref)(ref), and (ref), we have that the FNAR errors $\{\nu_t\}$ are not autocorrelated. This is necessary for the CLT in the next part of this assumption to hold. Assumptions (ref)(ref) and (ref)(ref) follow directly from (ref)(ref) and (ref)(ref) since we know that $V^{-1}$ is finite for all $N\in\mathbb N$. Notice, however, the different role of $N$ in the definition of $\Omega_0$ and $\Omega_1$, indeed, as $N\to\infty$ we have $X_t'VX_t=O_p(N^2)$, but $X_t' V^{-1} X_t=O_p(N)$, since $X_t=O_p(\sqrt N)$, $V=O(N)$ and $V^{-1}=O(1)$.

Second, to study the properties of the GLS estimator (ref) we need to prove consistency of the estimated inverse of the FNAR errors covariance $\widehat V^{-1}$ defined in (ref). This is not an easy task, for at least three reasons: first, the FNAR errors are estimated and not observed; second, the matrix $V$ is $N\times N$ so it is a high-dimensional one; third, to study $\widehat V^{-1}$, we need uniform consistency over all $N^2$ entries of the estimated covariance $\widehat V$. These difficulties are reduced if we assume that all considered random variables are sub-Gaussian, which is a classical assumption in high-dimensional statistics vershynin2018high.

This approach is similar to the one adopted in a vector factor model context by fan2013large. Instead, chenetal2020cnar assume a set of moment conditions on the regressors matrix $X_t$, which amount to bounding up the 8th order cross-cumulants (in addition, they also make use of Hanson-Wright concentration inequalities which are based on the assumption of sub-Gaussianity).

The GLS estimator in (ref) satisfies:

The asymptotic properties of the terms of (ref) are then given in the following Proposition.

The next theorem follows.

For observed network factors, the GLS estimator has a faster rate of convergence than the OLS estimator and it is more efficient; see also chenetal2020cnar. The different rates depend on the different scaling needed for Assumptions (ref)(ref) and (ref)(ref) to hold. Indeed, on the one hand $\mathbb{E} \left[ X_t' V X_t \right]=O(N^2)$, while, on the other hand $\mathbb{E} \left[ X_t' V^{-1} X_t \right]=O(N)$. This is because, by Assumption (ref), $\norm{V}=O(N)$ but $\norm{V^{-1}}=O(1)$.

Now, from Proposition (ref), we see that the conditions ${ \sqrt {NT}}/ m \to 0$ and $N/m\to 0$ allow us to treat the factors as observed. These conditions are stronger than in the OLS case.

An estimator of the asymptotic variance-covariance matrix of the GLS estimator is then given by: \[ \widehat{\operatorname*{Avar}}\left[ \sqrt{NT} (\widehat{\theta}^{\text{\upshape\tiny GLS}} -\bar J\theta) \right] = \left(\frac{1}{NT} \sum_{t=1}^{T} \widehat{X}_{t}^' \widehat{V}^{-1} \widehat{X}_t \right)^{-1}, \] where $ \widehat{V}^{-1}$ is defined in (ref). Alternatively, to address possible residual cross-correlation of the node idiosyncratic components, we could use: \[ \widehat{\operatorname*{Avar}}\left[ \sqrt{NT} (\widehat{\theta}^{\text{\upshape\tiny GLS}} -\bar J\theta) \right] = \left(\frac{1}{NT} \sum_{t=1}^{T} \widehat{X}_{t}^' \widehat{V}^{-1} \widehat{X}_t \right)^{-1} \left(\frac{1}{NT} \sum_{t=1}^{T} \widehat{X}_{t}^' \widehat{V}^{-1} \widehat\nu_t\widehat\nu_t' \widehat{V}^{-1} \widehat{X}_t \right) \left(\frac{1}{NT} \sum_{t=1}^{T} \widehat{X}_{t}^' \widehat{V}^{-1} \widehat{X}_t \right)^{-1}. \]

\refstepcounter{lettersection} \oldsection{Proofs of the main results}

General statements of Assumptions (ref), (ref), (ref), and (ref)

All the following results in Appendices (ref) and (ref) are proved under a more general version of the assumptions in the main text. Namely, we replace Assumptions (ref)(ref), (ref)(ref), (ref)(ref), (ref), (ref)(ref), and (ref) with:

Assumption 2 (General statement). {\it

}

Assumption 3 (General statement). {\it For all $i,j\in\mathbb N$, all $k=1,\ldots,r$, and all $t\in\mathbb Z$, $\mathbb E[\mathcal{F}_{(3)tkj}\mathcal{E}_{(3)tij}]=0$, and, for all $m,N,T\in\mathbb N$ and all $t=1,\ldots, T$,

for some finite $C_{\mathcal{F}\mathcal{E}}$ and $C_{\mathcal{F}\mathcal{E}}^\prime$ independent of $t$, $m, N$, and $T$ and some $\gamma \in [0,2]$.}

Assumption 5 (General statement). {\it

}

Assumption 10 (General statement). {\it For all $m,N,T\in\mathbb N$,

for some finite $\mathfrak K_1$ and $\mathfrak K_2$ independent of $m,N$, and $T$ and some $\gamma\in[0,2]$.}

These assumptions all depend on a generic $\gamma\in[0,2]$, thus are more general than those stated in the main text, which correspond to the case $\gamma=2$. All Theorems and Propositions in the main text and in Appendices (ref) and (ref) are stated in the case $\gamma=2$, which is the least favorable one, meaning it is the case giving the slowest possible convergence rates.

Proof of Proposition (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Proposition (ref)

Proof of Theorem (ref)

Proof of Proposition (ref)

Proof of Theorem (ref)

\refstepcounter{lettersection} \oldsection{Auxiliary lemmata}

We start with some notation. Let $\chi_t:=U\mathcal{F}_{(3)t}$ and recall the definition of the $m\times m$ matrices:

with $j$ largest eigenvalues $\widehat{\mu}_j^{\mathcal{W}}$, ${\mu}_j^{\mathcal{W}}$, ${\mu}_j^{\chi}$, and ${\mu}_j^{\mathcal{E}}$, respectively.