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.
146,244 characters · 14 sections · 99 citation commands
IV Estimation of Heterogeneous Spatial Dynamic Panel Models with Interactive Effects
\thispagestyle{fancy}
\baselineskip=15.0pt
Economic outcomes are shaped by complex dependencies that span both temporal dynamics and spatial interactions. Temporal dependencies arise from phenomena such as habit formation, adjustment costs, and economic slack, wherein past behavior influences current outcomes (e.g., Hamermesh1995, JappelliPistaferri2017). Spatial interactions, on the other hand, reflect the influence of peers, spatial networks, and spillover mechanisms (e.g., Case1991, Manski1993, BramulleEtal2009), where the behavior of one unit is affected by those of others. These dependencies are further complicated by the frequently pervasive impact of aggregate shocks, such as technological advancements, global market fluctuations, and economy-wide regulatory changes (e.g., SarafidisWansbeek2021). Together, these factors highlight the considerable challenges associated with modeling economic behavior.
Earlier contributions in the econometric panel data literature addressed temporal dynamics, spatial interactions and aggregate shocks largely in a fragmented manner. For instance, in the context of large-$T$ panels, a considerable body of work focused on dynamic models with additive fixed effects, with limited consideration of spatial interactions or aggregate shocks (e.g., HahnKuersteiner2002, AlvarezArellano2003, and Hayakawa2015). Over the past decade, progress has been made with the development of dynamic panels incorporating interactive effects to account for aggregate shocks. Examples include ChudikPesaran2015, MoonWeidner2017, NorkuteEtal2021, DeVosEveraert2021 and JuodisSarafidis2022. Parallel advancements have also been made in spatial dynamic panel data models with additive fixed effects, as explored by YuEtAl2008, Korniotis2010, LeeYu2014 among others. In these models, spatial interdependence is captured through a pre-specified $N \times N$ adjacency matrix $\mathbf{W}$, which encodes the structure of the interactions among individual units.\footnote{Comprehensive overviews of spatial panel data models with additive effects can be found in Elhorst2014 and LeeYu2015.}
Advances in econometrics have since sought to bridge these three strands by developing spatial dynamic panel data models that integrate temporal dependencies, spatial interactions, and interactive error components, as in ShiLee2017, BaiLi2021, CuiSarafidisYamagata2023, and HigginsMartellosio2023.\footnote{In HigginsMartellosio2023, it is assumed that $\mathbf{W}$ is observed only partially.} While these methodological advancements represent significant progress, a common assumption underlying much of this literature is slope-parameter homogeneity. That is, the magnitude of the relationships between dependent and independent variables is assumed to be identical across all cross-sectional units. Unfortunately, such assumption can be unduly restrictive in contexts characterized by substantial heterogeneity, which is a common feature of many economic systems. For instance, when cross-sectional heterogeneity in coefficients is captured by a random-coefficient model, it is well-established that dynamic pooled estimators fail to consistently estimate the population average, even for large $T$; see e.g., RobertsonSymons1992 and PesaranSmith1995.
To address this limitation, this paper develops a Mean Group Instrumental Variables (MGIV) estimator for spatial dynamic panel data models with interactive effects and heterogeneous slope coefficients, under large $N$ and $T$ asymptotics. Valid instruments are constructed by projecting out common factors from exogenous covariates using Principal Components Analysis (PCA), following the methodology of Bai2003. Consequently, the individual-specific IV estimates are consistent. Our MGIV estimator then combines these IV estimates and averages them to obtain consistent estimates of population-level effects.
The present extension to incorporate heterogeneous slope coefficients represents a major step forward in the spatial econometrics literature. As noted by LesageChih2016, “space-time panel data samples covering longer time spans allow us to produce parameter estimates for all $N$ spatial units, an exciting point of departure for future work. Allowing for heterogeneous coefficients for each spatial unit holds a natural appeal when contrasted with conventional static spatial panel models.”
In line with the aforementioned spatial literature, our analysis focuses on the estimation of heterogeneous slopes and spillover effects conditional on a pre-specified $\mathbf{W}$, which is treated as fixed and known. Consequently, the proposed method can be particularly appealing for panel datasets involving geographical entities, such as countries or administrative regions, where measures of spatial and economic proximity are readily available and may naturally inform the specification of $\mathbf{W}$. Furthermore, such datasets often include a substantial number of time-series observations per unit, as well as a large number of cross-sectional units, aligning well with the large $N,T$ asymptotics considered in this paper.\footnote{A complementary strand of the economics literature examines settings where $\mathbf{W}$ is latent and estimated from the data using high-dimensional methods, as in DePaula2017, LamSouza2020, and DePaulaEtal2024. This approach offers appealing generality by avoiding the need for a priori specification of $\mathbf{W}$. However, to date, this literature also uniformly assumes slope parameter homogeneity and typically excludes temporal dynamics and aggregate shocks, both of which are explicitly addressed in the present paper. Extending our framework to accommodate the estimation of an unknown $\mathbf{W}$ remains an avenue for future research.}
Our MGIV estimator is shown to be consistent and asymptotically normally distributed as $N,T \to \infty$. Importantly, MGIV is asymptotically unbiased. This property addresses the incidental parameters problem, and thereby standard inference procedures remain valid without the need for bias correction. Additionally, the estimator is linear and computationally efficient, making it particularly well-suited for empirical applications involving large datasets. Recently, our MGIV estimator was implemented in Stata by KripfganzSarafidis2025 via the spxtivdfreg command, providing researchers and practitioners with a readily available and user-friendly tool for empirical analysis.
The MGIV estimator developed in this paper extends the approach of NorkuteEtal2021 to incorporate spatial interactions, addressing issues related to the identification of heterogeneous spatial parameters and the development of asymptotic theory for large $N,T$ settings with spatially interdependent observations. Notably, our approach allows for more flexible expansion rates of $N$ and $T$ compared to NorkuteEtal2021, permitting $T$ to grow faster than $N$ (but slower than $N^{2}$), proportionally to $N$, or slower than $N$. This broader flexibility in the relative growth of $N$ and $T$ enhances the applicability of the estimator in diverse settings where data dimensions vary significantly.
Recently, ChenEtAl2022 proposed a related approach to our MGIV estimator, but their focus is limited to static panels without temporal dynamics.\footnote{Extending the approach of ChenEtAl2022 to dynamic models is non-trivial, as it requires the construction of factor proxies from suitable cross-sectional averages of the dependent variable, the number of which grows with $T$. As noted on page 56 of their paper, the finite-sample bias arising from the correlation between cross-sectional averages and idiosyncratic errors is likely to be exacerbated in the presence of spatial lags. This issue becomes especially pronounced when the spatial weighting matrix remains relatively dense as the sample size increases, which often occurs when $\mathbf{W}$ is specified based on geographic proximity.} Additionally, in their framework interactive effects are captured using cross-sectional averages \`a la Pesaran2006, which relies on the so-called rank condition (e.g., non-zero mean factor loadings). In contrast, our method remains asymptotically valid even if the rank condition is violated.\footnote{Although DeVosEtal2024 proposed a method for evaluating the rank condition for CCE estimators, it is currently applicable only to static panels and does not directly extend to settings with temporal dynamics.} Another related study is that of AquaroEtal2021, which assumes a purely idiosyncratic error structure without accounting for additive or interactive effects.
We illustrate the practical relevance of our method by estimating a regional spatial growth model, focusing on the magnitude of regional growth spillovers in Europe. Our approach explicitly accounts for heterogeneity across regions, reflecting differences in industrial structure, urbanization, and geographic characteristics. For instance, growth drivers in industrial Bavaria (Germany) likely differ from agricultural Thessaly (Greece), as do spillover effects between urban Greater London and rural Lapland (Finland). The results provide evidence of conditional convergence in regional growth dynamics. Spillovers play a dominant role, particularly for investment rates, where roughly four-fifths of the total impact on GDP per capita growth is attributable to neighboring regions, emphasizing the importance of inter-regional linkages. Human capital and R&D spillovers further reinforce the role of knowledge diffusion and innovation in fostering growth across regions.
Throughout, for an $m\times m$ matrix $\mathbf{A}=(a_{ij})_{1\leq i,j\leq m}$, we denote its trace by $\mathrm{tr} (\mathbf{A})=\sum_{i=1}^m a_{ii}$. For an $m\times n$ matrix $\mathbf{B}=(b_{ij})_{1\leq i\leq m,1\leq j\leq n}$, denote its column sum norm by $\|\mathbf{B}\|_1=\max_{1\le j\le n}\sum_{i=1}^m |b_{ij}|$, its Frobenius norm by $\|\mathbf{B}\|=\sqrt{\mathrm{tr} (\mathbf{B}^{\prime }\mathbf{B})}$, and its row sum norm by $\|\mathbf{B}\|_\infty=\max_{1\le i\le m}\sum_{j=1}^n|b_{ij}|$. Also define $\mathbf{P}_{\mathbf{B}}=\mathbf{ B}(\mathbf{B}^{\prime }\mathbf{B})^{-1}\mathbf{B}^\prime$ and $\mathbf{M}_{\mathbf{B}}=\mathbf{I}_m-\mathbf{P}_{\mathbf{B}}$, where $\mathbf{I}_m$ is the $m\times m$ identity matrix. Denote by $C$ a generic positive constant which need not be the same at each appearance, and $\delta_{N\!T}^2=\min\{N,T\}$. We use $N,T \rightarrow \infty$ to denote that $N$ and $T$ pass to infinity jointly, and $\operatorname*{plim}$ to denote the probability limit.
We consider the following spatial dynamic panel data model with $N$ cross-sectional units and $T$ time periods:
$y_{it}$ denotes the outcome of interest for individual unit $i$ at time $t$, while $\sum_{j=1}^N w_{ij}y_{jt}$ is the so-called “spatial-lag”, a weighted sum of neighbor outcomes, where $w_{ij}$ denotes the weight assigned to neighbor $j$ in relation to $i$. These weights capture the connectedness structure among individuals and are specified within the $N \times N$ adjacency matrix $\mathbf{W}=[w_{ij}]$.\footnote{As discussed earlier, this paper aligns with standard practice in the spatial literature by assuming that $\mathbf{W}$ is fixed and known. Accordingly, the focus is placed on estimating the magnitude of heterogeneous spillover effects, conditional on a predetermined network structure.} The vector $\mathbf{x}_{it}=\left(x_{1it},\dots, x_{kit}\right)^{\prime}$ of dimension $k \times 1$, contains observed characteristics for individual $i$ at time $t$. Additionally, in the composite error term $\mathbf{f}_t$ and $\boldsymbol{g}_t$ denote vectors of latent factors with dimensions $r_1 \times 1$ and $r_2 \times 1$, respectively, influencing $y_{it}$. These factors are associated with factor loadings $\boldsymbol{\lambda}_i$ and $\boldsymbol{\phi}_i$, which vary across individuals. The inclusion of latent factor components reflects the premise that individual agents inhabit a common economic environment and, as such, are subject to aggregate, economy-wide or “global” shocks that affect the entire population, albeit with different intensities. Examples of such shocks include technological disruptions, natural disasters, financial crises, pandemics, geopolitical conflicts, global market fluctuations and regulatory changes (e.g., Bai2009 and SarafidisWansbeek2012). Note that the number of parameters associated with factors and factor loadings increases with $T$ and $N$, respectively, as $N,T \to \infty$, leading to the presence of incidental parameters. Lastly, $\varepsilon_{it}$ is a purely idiosyncratic error term.
There are three primary sources of endogeneity in model (ref). First, $\sum_{j=1}^N w_{ij}y_{jt}$ is endogenous by construction. This term essentially represents the formal specification of an equilibrium outcome of a spatial interaction process, wherein the value of the dependent variable for one individual unit is simultaneously determined alongside that of its neighbours. This reciprocal interdependence underscores the networked nature of these relationships, as highlighted by Elhorst2021.
Second, the lagged dependent variable, $y_{i,t-1}$ is also endogenous due to the presence of incidental parameters (e.g., Nickell1981 and PhillipsSull2007).
Third, endogeneity also arises from the potential dependence of the covariates on the latent factors. To account for this, and consistent with the frameworks of Pesaran2006 and NorkuteEtal2021 among many others, we assume
That is, $\mathbf{x}_{it}$ is influenced by $\mathbf{f}_{t}$, with the associated loadings represented by $\boldsymbol{\Gamma}_{i}$, a $k \times r_{1}$ matrix. Importantly, for the sake of generality, we allow the latent factors governing $y_{it}$ and $\mathbf{x}_{it}$ to differ. This distinction justifies the inclusion of the additional term $\boldsymbol{\phi}_i^{\prime }\boldsymbol{g}_t$ in model (ref). Moreover, the loadings $\boldsymbol{\Gamma}_{i}$ are permitted to exhibit correlation with both $\boldsymbol{\lambda}_i$ and $\boldsymbol{\phi}_i$, further accommodating potential interdependencies that may arise in the model. Finally, $\mathbf{v}_{it}$ denotes the idiosyncratic error component for $\mathbf{x}_{it}$.
The individual-specific structural parameters reflect distinct mechanisms: $\rho_{i}$ captures habit formation and adjustment costs, facilitating an important distinction between short- and long-run responses. $\boldsymbol{\beta}_{i}$ reflects the direct effects of an individual's own characteristics, and $\psi_{i}$ encapsulates the influence of neighbours' outcomes, also known as spillover effects (e.g., KelejianPiras2017 and JingEtAl2018).\footnote{The main results of this paper naturally extend to models incorporating “contextual effects” (e.g., Manski1993), also referred to as the spatial Durbin model (Elhorst2014). This is further discussed in Remark (ref).}
Stacking the equations in ((ref)) over $t$ yields
where $\mathbf{y}_{i}=(y_{i1},\ldots,y_{iT})^{\prime }$, $\mathbf{y}_{j}=(y_{j1},\ldots,y_{jT})^{\prime}$, $\mathbf{y} _{i,-1}=(y_{i0},\ldots,y_{i,T-1})^{\prime }$, $\mathbf{X}_{i}=(\mathbf{x}_{i1},\cdots, \mathbf{x}_{iT})^{\prime }$, $\mathbf{F}=(\mathbf{f}_1,\cdots,\mathbf{f} _T)^{\prime }$, $\mathbf{G}=(\boldsymbol{g}_1,\cdots,\boldsymbol{g} _T)^{\prime }$, $\boldsymbol{\varepsilon}_{i}=(\varepsilon_{i1},\cdots, \varepsilon_{iT})^{\prime }$ and $\mathbf{V}_{i}=(\mathbf{v}_{i1},\ldots, \mathbf{v}_{iT})^{\prime }$.
Define $\boldsymbol{\theta} _i=(\psi_i,\rho_i,\boldsymbol{\beta}_i^{\prime })^{\prime }$, $\mathbf{C} _{i}=(\sum_{j=1}^{N} w_{ij} \mathbf{y}_{j},\mathbf{y}_{i,-1},\mathbf{X}_{i})$, and $\mathbf{u}_{i}=\mathbf{F}\boldsymbol{\lambda}_i+\mathbf{G}\boldsymbol{ \phi}_i+\boldsymbol{\varepsilon}_{i}$. Then, the first equation in ((ref)) can be reformulated as
We use the method of Instrumental Variables to estimate $\boldsymbol{\theta}_{i}$. To this end, define the “defactoring” matrices that project out $\mathbf{F}$ and $\mathbf{F}_{-1}=(\mathbf{f}_0,\ldots,\mathbf{f}_{T-1})^\prime$ as $\mathbf{M}_{{\mathbf{F}}}=\mathbf{I}_T - \mathbf{ F}(\mathbf{F}^{\prime }\mathbf{F})^{-1}\mathbf{F}^\prime$ and $\mathbf{M}_{{\mathbf{F}_{-1}}}=\mathbf{I}_T - \mathbf{ F}_{-1}(\mathbf{F}_{-1}^{\prime }\mathbf{F}_{-1})^{-1}\mathbf{F}_{-1}^\prime$. Letting $\mathbf{X}_{i,-1}=(\mathbf{x}_{i0},\ldots,\mathbf{x}_{i,T-1})^\prime$, further define
whose elements are instruments for the elements of $\mathbf{C}_i$.\footnote{Loosely speaking, the term $\sum_{j=1}^Nw_{ij}\mathbf{M}_{\mathbf{F}}\mathbf{X}_{j}$ instruments $\sum_{j=1}^{N} w_{ij} \mathbf{y}_{j}$, the term $\mathbf{M}_{\mathbf{F}} \mathbf{M}_{\mathbf{{F}}_{-1}}\mathbf{X}_{i,-1}$ instruments $\mathbf{y}_{i,-1}$, and ${\mathbf{M}_{\mathbf{{F}}}}\mathbf{X}_{i}$ instruments $\mathbf{X}_{i}$.} It's straightforward to see that, due to the de-factorisation, we have
where $\mathbf{V}_{i,-1}=(\mathbf{v}_{i0},\ldots,\mathbf{v}_{i,T-1})^\prime$. Therefore, these instruments are exogenous.
Since $\mathbf{F}$ and $\mathbf{F}_{-1}$ are not observed, they are estimated using PCA on $\mathbf{X}$ and $\mathbf{X}_{-1}$. Assuming $T^{-1}\mathbf{F}^\prime\mathbf{F}= \mathbf{I}_{r_{1}}$ and $T^{-1}\mathbf{F}_{-1}^\prime\mathbf{F}_{-1}= \mathbf{I}_{r_{1}}$, the estimates $\widehat{\mathbf{F}}$ and $\widehat{\mathbf{F}}_{-1}$ are obtained as $\sqrt{T}$ times the eigenvectors corresponding to the $r_1$ largest eigenvalues of the $T \times T$ matrices $(NT)^{-1}\sum_{i=1}^N \mathbf{X}_{i}\mathbf{X} _{i}^{\prime }$ and $(NT)^{-1}\sum_{i=1}^N \mathbf{X}_{i,-1}\mathbf{X}_{i,-1}^{\prime }$, respectively.\footnote{For simplicity and without loss of generality, $r_1$ is treated as known. However, in practice it can be consistently estimated using established methods in the literature, such as the information criterion approach of BaiNg2002 or the eigenvalue methods of AhnHorenstein2013.} Note that, since the factors are extracted from observed covariates, no estimation error arises in $\widehat{\mathbf{F}}$ that is associated with estimation of the slope coefficients. Furthermore, the factor loadings $\boldsymbol{\Gamma}_i$ can be estimated as $\widehat{\boldsymbol{\Gamma}}_i=T^{-1}\widehat{ \mathbf{F}}^{\prime }\mathbf{X}_{i}$.
Feasible instruments for $\boldsymbol{\theta}_{i}$ are constructed as
and the resulting IV estimator of $\boldsymbol{\theta}_i$ is given by
where
In this section, we examine the limiting properties of the individual-specific IV estimator, $\hat{\boldsymbol{\theta}_i}$, and the MGIV estimator. To facilitate the asymptotic analysis, we first introduce a set of assumptions.
Assumption (ref) is in line with existing spatial literature, see e.g. LeeYu2014, CuiSarafidisYamagata2023. Cross-sectional and time-series homoskedasticity is imposed to simplify the asymptotic analysis of the variance-covariance estimator in panels where both $N$ and $T$ are large. In contrast, NorkuteEtal2021 allow for cross-sectional/time-series heteroskedasticity by leveraging the results in Hansen2007. However, Hansen2007 assumes independence across cross-sectional units, a condition that is violated in the present setup. Although we do not formally derive theoretical results under heteroskedasticity, the finite-sample performance of a robust variance-covariance estimator is thoroughly investigated in Section (ref).
Assumption (ref) ensures that $\mathbf{x}_{it}$ is strictly exogenous with respect to $\varepsilon_{it}$, as e.g., in Pesaran2006 and Bai2009. Therefore, the defactored regressors are valid instruments. Moreover, this assumption accommodates cross-sectional and time series heteroskedasticity, as well as autocorrelation in $\mathbf{v}_{it}$. Unlike $\varepsilon_{it}$, here it is crucial to explicitly allow for this broader structure since, conditional on $\mathbf{F}$, the dynamics in $\mathbf{X}_{i}$ are driven by $\mathbf{V}_{i}$. Additionally, in contrast to NorkuteEtal2021, $\mathbf{v}_{it}$ is allowed to exhibit weak cross-sectional correlation, aligning with the assumption of weak dependence in the process of $y$.
Assumptions (ref) and (ref) align with standard conditions in the PCA literature; see, for example Bai2003. Assumption (ref) permits both interdependence between $\mathbf{f}_{t}$ and $\mathbf{g}_{t}$, as well as within-group correlations in each term. Similarly, Assumption (ref) allows for non-zero correlations not only between $\boldsymbol{\Gamma}_{i}$, $\boldsymbol{\varphi}_{i}$ and $\boldsymbol{\lambda}_{i}$, but also within each of these components. These assumptions are particularly relevant in settings where $y_{it}$ and $\mathbf{x}_{it}$ may be jointly influenced by common shocks.
Assumption (ref) is standard in the spatial literature, as outlined in KelejianPrucha2001. Specifically, Assumption (ref)(ref) serves as a normalization condition, ensuring that no individual is treated as its own neighbor. Assumptions (ref)(ref)-(ref)(ref) ensure the absence of a dominant unit, i.e., a unit that becomes asymptotically correlated with all others. Such scenarios are instead accommodated by the inclusion of latent factors. Assumption (ref)(ref) pertains to the space of the autoregressive and spatial parameters, and are discussed in detail by KelejianPrucha2010. Importantly, these assumptions are invariant to the ordering of the data, which can be arbitrary as long as Assumption (ref) is satisfied. Additionally, it is noteworthy that the spatial weighting matrix $\mathbf{W}$ is not required to be row-normalized.
Assumption (ref) ensures IV-based identification, as e.g. in Wooldridge2002. For instance, Assumption (ref)(ref) implies that the covariates and instruments are not perfectly collinear, and that their correlation is sufficiently strong to avoid degeneracy.
Finally, Assumption (ref) is commonly employed in the random coefficients literature (e.g., Pesaran2006). Assumption (ref)(ref) imposes upper bounds on the random coefficients, ensuring the stability of the process.
The following theorem demonstrates the asymptotic properties of the individual-specific IV estimator $\widehat{\boldsymbol{\theta}}_i$.
Note that the large-$N$ requirement is indispensable because the validity of the instruments used by $\widehat{\boldsymbol{\theta}}_{i}$ necessitates consistent estimation (up to rotation) of the $T \times r_{1}$ matrix $\mathbf{F}$.
The relative expansion rate of $N$ and $T$ employed in Theorem (ref) is more general than that in NorkuteEtal2021, which imposes $T/N\rightarrow c$, where $0 < c < \infty$. Specifically, the relative expansion rate adopted here allows for greater flexibility, permitting $T$ to grow faster than $N$ (but slower than $N^{2}$), proportionally to $N$, or slower than $N$.\footnote{It is straightforward to verify that any case where $T=cN$ also satisfies $T=o\left(N^{2} \right)$ but the converse does not necessarily hold.}
Assuming the cross-sectionally heterogeneous coefficients $\boldsymbol{\theta}_i$ follow the random-coefficient model, as in Assumption (ref), it is known that the dynamic pooled estimator of the population average $\boldsymbol{\theta}=\mathbb{E}(\boldsymbol{\theta}_i)$ will be inconsistent (see RobertsonSymons1992, PesaranSmith1995 and ChenEtAl2022). In the present case, the same holds true even in the absence of temporal dynamics in model (ref) because the spatial lag variable, $\sum_{j=1}^{N} w_{ij} y_{jt}$, is endogenous by construction. For this reason, we develop a Mean Group IV estimator of $\boldsymbol{\theta}$, which combines the individual-specific IV estimates and averages them to obtain consistent estimates of population-level effects.
Specifically, once $\widehat{\boldsymbol\theta}_i$ in ((ref)) are obtained, the Mean Group IV (MGIV) estimator of $\boldsymbol{\theta}$ is obtained as
Theorem (ref) below establishes the asymptotic properties of $\widehat{\boldsymbol\theta}_{MG}$.
We investigate the finite sample behaviour of the proposed approach by means of Monte Carlo experiments. We shall focus on the mean, bias, RMSE, empirical size and power of the t-test.
{0pt} {0pt} {0pt} {0pt} We consider the following heterogeneous, spatial dynamic panel data model:
$i=1,...,N$, $t=-49,...,T$, where
with ${\zeta_{st}\sim i.i.d.N(0,1)}$ for $s=1,...,r_{y}$. We set $\varrho_{fs}=0.5$ $\forall s$, $k=2$ and $r_{y}=3$.
The spatial weighting matrix, $\mathbf{W}=[w_{ij}]$ is an invertible rook matrix of circular form (KapoorKelejianPrucha2007), such that its $i$th row, $1<i<N$, has non-zero entries in positions $i-1$ and $i+1$, whereas the non-zero entries in rows $1$ and $N$ are in positions $(1,2)$, $(1,N)$, and $(N,1)$, $(N,N-1)$, respectively. This matrix is row normalized so that all of its nonzero elements are equal to $1/2$.
The idiosyncratic error, $\varepsilon_{it}$, is non-normal and heteroskedastic across both $i$ and $t$, such that $\varepsilon_{it}=\varsigma_{\varepsilon}\sigma_{it}(\epsilon_{it}-1)/\sqrt{2}$, $\epsilon_{it}\sim i.i.d.\chi_{1} ^{2},$ with $\sigma_{it}^{2}=\eta_{i}\phi_{t}$, $\eta_{i}\sim i.i.d.\chi _{2}^{2}/2$, and $\phi_{t}=t/T$ for $t=0,1,...,T$ and unity otherwise.
{0pt} {0pt} {0pt} {0pt} The stochastic process for the covariates is given by
for $\ell=1,2$. We set $r_{x}=2$. Thus, the first two factors in $u_{it}$, {$f_{1t},f_{2t}$}, also drive the DGP for $x_{\ell it}$, $\ell=1,2$. However, $f_{3t}$ does not enter into the DGP of the covariates directly.\footnote{Observe that, using notation of earlier sections, $f_{3t}=g_{t}$ in Eq. (ref).}
{0pt} {0pt} {0pt} {0pt} The idiosyncratic errors in the covariates are serially correlated, such that
for $\ell=1,2$. We set $\varrho_{\upsilon,\ell}=\varrho_{\upsilon}=0.5$ for all $\ell$.
{0pt} {0pt} {0pt} {0pt} All individual-specific effects and factor loadings are generated as correlated and mean-zero random variables. In particular, the individual-specific effects are drawn as
where $\omega_{\ell i}\sim i.i.d.N(0,\left( 1-\varrho \right)^{2})$, for $\ell=1,2$. We set $\varrho_{\mu,\ell}=0.5$ for $\ell=1,2$.
{0pt} {0pt} {0pt} {0pt} The factor loadings in $u_{it}$ are generated as $\varphi_{si}\sim i.i.d.N(0,1)$ for $s=1,...,r_{y}(=3)$, and the factor loadings in $x_{1it}$ and $x_{2it}$ are drawn as
respectively, for $s=1,...,r_{x}=2$. The process in Eq. ((ref)) allows the factor loadings to $f_{1t}$ and $f_{2t}$ in $x_{1it}$ to be correlated with the factor loadings corresponding to the factor that does not enter into the DGP of the covariates, i.e., $f_{3t}$. On the other hand, Eq. ((ref)) ensures that the factor loadings to $f_{1t}$ and $f_{2t}$ in $x_{2it}$ are correlated with the factor loadings corresponding to the same factors in $u_{it}$, $f_{1t}$ and $f_{2t}$. We consider $\varrho_{\gamma,11}=\varrho_{\gamma,12} \in \left \{0\text{, }0.5\right \}$, whilst $\varrho_{\gamma,21}=\varrho_{\gamma,22}=0.5$.
It is straightforward to see that the average variance of $\varepsilon_{it}$ depends only on $\varsigma_{\varepsilon}^{2}$. Let $\pi_{u}$ denote the proportion of the average variance of $u_{it}$ that is due to $\varepsilon_{it}$. That is, we define $\pi_{u}:=\varsigma_{\varepsilon}^{2}/\left( r_{y} +\varsigma_{\varepsilon}^{2}\right)$. Thus, for example, $\pi_{u}=3/4$ means that the variance of the idiosyncratic error accounts for 75% of the total variance in $u$. In this case most of the variation in the total error is due to the idiosyncratic component and the factor structure has relatively minor significance.
{0pt} {0pt} {0pt} {0pt} Solving in terms of $\varsigma_{\varepsilon}^{2}$ yields
We set $\varsigma_{\varepsilon}^{2}$ such that $\pi_{u}\in \left \{1/4\text{, }3/4\right \}$.\footnote{These values of $\pi_u$ are motivated by the results in Sargent_Sims_1977, in which they find that two common factors explain $86\%$ of the variation in unemployment rate and $26\%$ of the variation in residential construction.}
The slope coefficients are generated as $\rho _{i}=\rho +\eta _{\rho, i}$, $\psi_{1i}=\beta _{1}+\eta _{\psi, i}$, $\beta _{1i}=\beta _{1}+\eta _{\beta _{1},i}$ and $\beta _{2i}=\beta _{2}+\eta _{\beta _{2},i}$. We set $\rho=0.4$, $\psi=0.25$, and $\beta_{1}=3$ and $\beta_{2}=1$, following Bai2009. In addition, we specify $\eta _{\rho, i}\sim i.i.d.~U\left[ -c_{\rho},+c_{\rho}\right] $, $\eta _{\psi, i}\sim i.i.d.~U\left[ -c_{\psi},+c_{\psi}\right] $ and
where $\xi _{\beta_{\ell},i}$ is the standardised squared idiosyncratic error in $x_{\ell it}$, computed as
with $\overline{v_{\ell i}^{2}}=T^{-1}\sum_{t=1}^{T}v_{\ell it}^{2}$, $\overline{ v_{\ell }^{2}}=N^{-1}\sum_{i=1}^{N}\overline{v_{\ell i}^{2}}$, for $\ell =1,2$. We set $c_{\rho}=0.2$, $c_{\psi}=0.15$, $\rho_{\beta}=0.4$ for $\ell=1,2$.
{0pt} {0pt} {0pt} {0pt} We define the signal-to-noise ratio (SNR) conditional on the factor structure, the individual-specific effects and the spatial lag, as follows:
where $\mathcal{L}$ is the information set that contains the factor structure, the individual-specific effects and the spatial lag\footnote{The reason for conditioning on these variables is that they influence both the composite error of $y_{it}$, as well as the covariates.}, whereas $\overline{\text{var}}\left( \varepsilon_{it}\right) $ is the overall average of $E\left( \varepsilon _{it}^{2}\right) $ over $i$ and $t$. Solving for $\varsigma_{\upsilon}^{2}$ yields
We set $SNR=4$, following JuodisSarafidis2018ER and CuiSarafidisYamagata2023.
We consider two sets of instruments, namely
and
Both matrices are of dimension $T\times 4K$. The difference between the two matrices is that $\mathbf{\widehat{Z}}^{1}_{i}$ uses an extra (second) lag of $\mathbf{M}_{\widehat{\mathbf{F}}}\mathbf{X}_{i}$, which helps identification of the autoregressive parameter, whereas $\mathbf{\widehat{Z}}^{2}_{i}$ replaces the second lag of $\mathbf{M}_{\widehat{\mathbf{F}}}\mathbf{X}_{i}$ with the first lag of $\sum_{j=1}^Nw_{ij}\mathbf{M}_{\widehat{\mathbf{F}}}\mathbf{X}_{j}$.
The mean group estimator of $\boldsymbol{\theta}$ is defined as
where $\boldsymbol{\hat{\theta}}_{i}$ is given by
where
with $\ell \in \{1,2\}$ and $\mathbf{C}_{i}=(\mathbf{y}_{i,-1},\mathbf{X}_{i},\mathbf{Y}\mathbf{w}_i)$.
The variance-covariance matrix of the MG estimator is given by
As a benchmark estimator, we consider the pooled two-stage IV (2SIV) estimator developed by CuiSarafidisYamagata2023, which is designed for spatial dynamic panel data models with homogeneous parameters. This estimator is defined as follows:
where \[\widetilde{\mathbf{A}}=\frac{1}{NT}\sum_{i=1}^N \widehat{\mathbf{Z}}_{i}'\mathbf{M}_{\widehat{\mathbf{H}}}\mathbf{C}_{i}\,,\widetilde{\mathbf{B}}=\frac{1}{NT}\sum_{i=1}^N \widehat{\mathbf{Z}}_{i}'\mathbf{M}_{\widehat{\mathbf{H}}}\widehat{\mathbf{Z}}_{i}\,,\widetilde{\mathbf{c}}_y=\frac{1}{NT}\sum_{i=1}^N \widehat{\mathbf{Z}}_{i}'\mathbf{M}_{\widehat{\mathbf{H}}}\mathbf{y}_{i}\,,\] and $\mathbf{M}_{\widehat{\mathbf{H}}}=\mathbf{I}-\widehat{\mathbf{H}}(\widehat{\mathbf{H}}'\widehat{\mathbf{H}})^{-1}\widehat{\mathbf{H}}'$, with $\widehat{\mathbf{H}}$ defined as $\sqrt{T}$ times the eigenvectors corresponding to the $r_{y}$ largest eigenvalues of the $T \times T$ matrix $(NT)^{-1}\sum_{i=1}^N \widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}'$, where $\widehat{\mathbf{u}}_{i}$ is the residual of the first-stage homogeneous IV estimator.
In terms of the sample size, we consider three cases. Case I specifies $N=100\tau$ and $T=25\tau$ for $\tau=1,2,4$. This implies that while $N$ and $T$ increase by multiples of 2, the ratio $N$ over $T$ remains equal to $4$ in all circumstances. Case II specifies $T=100\tau$ with $N=25\tau$ for $\tau=1,2,4$. Therefore, $N/T=0.25$, as both $N$ and $T$ grow. Finally, Case III sets $N=T=50\tau $, $\tau=1,2,4$. These choices allow us to consider different combinations of $(N,T)$ in relatively small and large sample sizes. Note that Case I implies that the rate at which $N,T\to \infty$ violates the conditions outlined in Theorem (ref) for MGIV. Studying the performance of the estimator under these circumstances is valuable, as in many applications $N$ can be significantly larger than $T$. This is exemplified by our empirical application discussed in the following section.
All results are obtained based on 2,000 replications, and all tests are conducted at the 5% significance level. For the power of the “t-test”, we specify $H_0:\rho=\rho^{0}+0.1$ (or $H_0:\psi=\psi^{0}+0.1$, and $H_0:\beta_{\ell}=\beta_{\ell}^{0}+0.1$ for $\ell=1,2$) against two sided alternatives, where ${\rho^0,\psi^0,\beta_1^0,\beta_2^0}$ denote the true parameter values.
Tables (ref)--(ref) report results for $\rho=0.4$, $\psi=0.25$ and $\beta_{2}=1$ in terms of the Mean, RMSE, ARB (Absolute Relative Bias), as well as empirical size (nominal level is $5\%$) and size-corrected power, which is computed based on the 2.5% and 97.5% quantiles of the empirical distribution of the t-ratio under the null hypothesis.\footnote{Results for $\beta_{1}=3$ are qualitatively similar to those for $\beta_{2}=1$ and therefore they are reported in the Appendix.} ARB is defined as $ARB \equiv \left(|\widehat{\theta}_{\ell}-\theta_{\ell}|/\theta_{\ell}\right)100$, where $\theta_{\ell}$ denotes the $\ell$th entry of $\boldsymbol{\theta}=\left(\psi, \rho, \boldsymbol{\beta}' \right)'$. “IVMG” denotes the Mean Group IV estimator of $\boldsymbol{\theta}$ as defined in (ref) with the matrix of instruments given by $\mathbf{\widehat{Z}}^{1}_{i}$.\footnote{The results for the IVMG and 2SIV estimators based on $\mathbf{\widehat{Z}}^{2}_{i}$ are reported in the Appendix.} In each of the tables, Panel A corresponds to $\pi_{u}=3/4$ and Panel B to $\pi_{u}=1/4$.
In regards to Table (ref), IVMG appears to have very little bias under all circumstances. Furthermore, its ARB values decrease steadily with larger sample sizes. Empirical size is very close to nominal one; even in cases where $T$ is rather small, there are only mild size distortions. Moreover, size approaches $5\%$ quickly as both $N$ and $T$ increase.
On the other hand, 2SIV appears to be severely biased. This is not surprising given that the pooled IV estimator is inconsistent under slope parameter heterogeneity. In particular, ARB increases with larger sample sizes and empirical size is heavily distorted throughout. These properties hold true regardless of the value of $\pi_{u}$.
The overidentifying restrictions test statistic (J test) associated with the optimal 2SIV estimator has satisfactory power to reject the null hypothesis of slope parameter homogeneity, especially in relatively larger samples. It is worth noting that power increases dramatically when more instruments with respect to $\mathbf{M}_{\widehat{\mathbf{F}}}\mathbf{X}_{i}$ are used. For example, when $\mathbf{M}_{\widehat{\mathbf{F}}} \mathbf{M}_{\widehat{\mathbf{F}}_{-3}}\mathbf{X}_{i,-3}$ is added to the existing set of instruments in Eq. (ref), then for $\pi_{u}=3/4$ and $ T=100\tau$, $N=25\tau$, $\tau=1,2,4$, power increases to $33.25\%$, $51.45\%$ and $94.7\%$, respectively, compared to $12.0\%$, $23.5\%$ and $75.1\%$ reported in Table (ref) under Panel A.\footnote{Detailed results for the case of this expanded set of instruments are available from the authors upon request.}
The results in Tables (ref)--(ref) are qualitatively similar to those in Table (ref), albeit the 2SIV pooled estimator exhibits smaller ARB for $\psi$ and $\beta_{2}$ compared to $\rho$. IVMG typically outperforms 2SIV in terms of ARB and has superior size properties.
In summary, MGIV performs well in all circumstances. In particular, the finite-sample bias of the estimator is negligible in almost all cases examined, and inferences are credible even with relatively small samples. Notably, the performance of the estimator appears to be robust to a wide range of values for $\pi_u$, which captures the proportion of variation in the total error is due to the idiosyncratic error. Thus, considering also the computational simplicity of our methodology, our IVMG estimator presents an attractive estimation approach in heterogeneous spatial dynamic panel data models with interactive effects when both $N$ and $T$ are large.
This section demonstrates our methodology by estimating a heterogeneous spatial growth model across EU regions over the period 2001-2018. The analysis focuses on uncovering the magnitude and influence of regional growth spillovers, a topic of extensive investigation in the literature.
Estimating the magnitude of growth spillovers has been a focal point in regional economic studies. ErturCoch2007 laid the groundwork with a spatially augmented neoclassical growth model, integrating productivity spillovers driven by capital investment. Their framework is supported by extensive empirical evidence on the role of knowledge and technological diffusion (e.g., AudretschFeldman2004, AutantBernardLeSage2011).
Further studies have refined methods to measure spillovers at the sub-national level, often focusing on urban regions. For example, BotazziPeri2003 and RodriguezRoseCrescenzi2008 analysed the spatial extent of innovation spillovers, while FunkeNiebuhr2005 assessed inter-regional economic interactions using spatial econometrics. Despite consensus on the existence of spillovers, their estimated magnitudes remain contentious. For instance, studies such as RamajoEtal2008, BenosEtal2015, OzyurtDees2018 and ElhorstEtal2024 report mixed findings, often influenced by methodological factors, including the choice of the weighting matrix and the temporal coverage of the dataset.\footnote{Differences in empirical methods also contribute to varying results; while growth regressions typically suggest convergence, tests of distribution dynamics point toward divergence (see e.g., DiVaioEtal2014 for a high-level discussion).}
Building on this literature, this section revisits the empirical specification proposed by ErturCoch2007 and ElhorstEtal2024, which include a spatial-time-lag, as well as spatial lags of the covariates. A key point of departure from these studies lies in allowing for a heterogeneous model that accommodates region-specific parameters. Accounting for slope parameter heterogeneity is crucial because European regions are inherently different. For example, industrially advanced regions such as Bavaria, may exhibit distinct economic dynamics compared to less industrialized or predominantly agricultural regions like Thessaly in Greece. Similarly, highly urbanized areas like Greater London may experience fundamentally different growth drivers and constraints compared to sparsely populated regions like Lapland in Finland.
Another significant departure from ErturCoch2007 and ElhorstEtal2024 is the inclusion of interactive fixed effects to capture unobserved common shocks and their heterogeneous impacts across regions. This is particularly relevant as the dataset spans the period 2001-2018, encompassing major economic events such as the 2008-2009 financial crisis, the Great Recession, and the European debt crisis (2010-2015).\footnote{The dataset is publicly available in Paul Elhorst's website: \href{https://spatial-panels.com/}{https://spatial-panels.com/}.} These events likely induced common shocks that affected multiple regions simultaneously, albeit to varying degrees. Purging the effects of these exogenous shocks is therefore essential for obtaining reliable estimates of growth parameters.
The empirical growth model is specified as:
where the dependent variable, $\Delta ln \left(y_{it}\right)$, represents the change in the natural logarithm of real GDP per capita (adjusted for purchasing power parity) in region $i$ at time $t$, with $i=1,\dots,266$ and $t=2001,\dots,2018$.\footnote{The analysis focuses on EU NUTS-2 regions, which provide a harmonized dataset on GDP and other macroeconomic indicators, ensuring cross-regional comparability and offering a higher resolution than national aggregates. While data at the NUTS-3 level exist, they are often incomplete or inconsistent, making NUTS-2 the most suitable choice.}
Thus, the growth rate of real GDP per capita in region $i$ at time $t$ is influenced by the initial level of GDP per capita in region $i$, the contemporaneous growth rate of GDP per capita in neighboring regions and its lag,\footnote{The lagged own growth rate was found to be statistically insignificant throughout and was therefore excluded from the model.} and a set of $k=4$ explanatory variables along with their spatial counterparts.
These covariates are defined as follows:
$x_{1it} \equiv ln\left(inv_{it}\right)$, where $inv$ denotes the investment rate (investment as a share of GDP). Investment drives capital accumulation and serves as a proxy for the resources allocated to future production capacity;
$x_{2it} \equiv ln\left(n_{it}\right)=ln\left(p_{it}+g+q\right)$, where $p_{it}$ represents the population growth rate for region $i$ at time $t$, while $g$ and $q$ denote the rates of technological progress and capital depreciation, respectively.\footnote{In line with the standard assumptions of the neoclassical growth framework, the rates of technological progress (g) and depreciation (q) are not indexed, as they are assumed to remain uniform across all regions and time periods. Following Islam1995, we set g + q = 0.05.} Population growth influences capital dilution and, indirectly, the region's convergence speed toward its steady state;
$x_{3it} \equiv ln\left(educ_{it}\right)$, where $educ$ denotes the share of the working-age population with tertiary education. This variable captures human capital, potentially a key driver of productivity and innovation in both neoclassical and endogenous growth frameworks; and
$x_{4it} \equiv ln\left(sci\&tech_{it}\right)$, where $sci\&tech$ represents the share of employment in science and technology. This variable serves as a proxy for regional engagement in innovation-driven activities, reflecting the intensity of knowledge-based economic development. The inclusion of $x_{3it}$ and $x_{4it}$ aligns with endogenous growth theory, which emphasizes the role of human capital and technological innovation as engines of long-term growth. In contrast, Solow's neoclassical framework considers these factors as exogenous.
The error term is composite and given by:
Conditional convergence occurs when $\theta_{i}<0$, implying that regions with lower (higher) initial GDP per capita tend to experience faster (slower) growth, holding all else constant. In contrast, conditional divergence arises when $\theta_{i}>0$, indicating that initial disparities in GDP per capita widen over time. A special case occurs when $\theta_{i}=0$, which YuEtAl2012 refer to as “spatial cointegration”. This scenario suggests that while GDP per capita growth rates may fluctuate across economies over the business cycle, regional economies ultimately remain on distinct growth trajectories throughout the sample period.
The spatial weights matrix is specified as:
where $d_{ij}$ represents the great-circle distance between regions $i$ and $j$.\footnote{This is measured in kilometers and constructed based on the latitude and longitude coordinates of the regions' centroids} Finally, $\zeta$ is the distance decay parameter. Following Ertur and Koch (2007), $\zeta$ is set to $0.02$, which prioritises interactions among proximate regions. Alternative specifications for $\mathbf{W}$, such as inverse squared distances, are explored in Section (ref).
The set of instruments is given by
where $\mathbf{X}_{i}=(\mathbf{x}_{i1},\cdots,\mathbf{x}_{iT})^{\prime }$ with $\mathbf{x}_{it}=\left(x_{1it},\dots,x_{4it}\right)^{\prime}$. Thus, a total of 20 instruments is utilized.
Table (ref) below provides results for the model specified in Eq. (ref), alongside several nested models that impose progressively stronger restrictions. Each column corresponds to a distinct specification, revealing the effects of these restrictions on the parameter estimates and their interpretations.
Column [6] represents the most general specification, allowing for spatial spillovers, interactive fixed effects, and slope-parameter heterogeneity. On the other hand, Column [1] represents the most restricted model. This specification imposes $u_{it}=\alpha_{i}+\tau_{t}+\varepsilon_{it}$, $\psi_{0,i}=\psi_{1,i}=\gamma_{\ell,i}=0$ for all $i$ and $\ell$, $\theta_{i}=\theta$ and $\beta_{\ell,i}=\beta_{\ell}$ for all $\ell$ in Eq. (ref).\footnote{The former restriction on the error term is equivalent to imposing one factor only in Eq. (ref) with the corresponding loadings being constant across $i$.} In essence, this formulation rules out spatial spillover effects, enforces slope-parameter homogeneity, and assumes that the error structure adheres to a two-way error components model. The columns in-between present intermediate cases. Column [2] extends Column [1] by incorporating latent common factors in the error term, relaxing the restriction of an additive two-way error components structure. Column [3] relaxes slope-parameter homogeneity while maintaining no spatial spillovers or interactive effects. Column [4] introduces spatial spillovers, while retaining the restrictions $u_{it}=\alpha_{i}+\tau_{t}+\varepsilon_{it}$, $\theta_{i}=\theta$, $\psi_{\tau,i}=\psi_{\tau}$ for all $\tau$ and $\beta_{\ell,i}=\beta_{\ell}$ for all $\ell$ in Eq. (ref). Finally, Column [5] builds on Column [4] by allowing for interactive fixed effects, although it still rules out slope-parameter heterogeneity.\footnote{The results in Columns [1]-[3] are obtained using the Stata command xtivdfreg, developed by KripfganzSarafidis2021, while those in Columns [4]-[6] are based on the Stata wrapper command spxtivdfreg, developed by KripfganzSarafidis2025.}
We start with the most restrictive specification first, Column [1]. The coefficient on $ln \left(y_{t-1}\right)$ is negative and highly significant, supporting the hypothesis of conditional convergence among regions. This finding aligns with the Solow growth model, which predicts that regions with higher initial GDP per capita grow more slowly as they approach their steady-state income levels. The coefficient on $ln\left(inv_{t}\right)$ is positive and significant, indicating that higher investment rates contribute positively to regional growth. Conversely, the coefficient on $ln\left(n_{t}\right)$ is negative and significant, consistent with neoclassical growth theory's prediction that higher population growth reduces per capita output by diluting capital per worker. The coefficients on educational attainment and share of employment in science and technology are both close to zero and statistically insignificant. The J test strongly rejects the null hypothesis of valid instruments, indicating potential model misspecification.
Allowing for interactive effects in the error process, as in Column [2], or slope-parameter heterogeneity, as in Column [3], leads to a more negative coefficient on $ln \left(y_{t-1}\right)$, which implies faster conditional convergence. This suggests that accounting for unobserved global shocks or shared regional trends and slope heterogeneity across regions is important.
Turning into the results incorporating spatial spillovers: The coefficients on $\mathbf{W} \Delta ln \left(y_{t}\right)$ and $\mathbf{W} \Delta ln \left(y_{t-1}\right)$ in Columns [4] and [5] imply the presence of a spatial unit root. Taken literally, this suggests that shocks to neighboring regions' growth rates persist indefinitely and propagate throughout the spatial network without decay, which is counterintuitive. At the same time, the coefficient of $ln \left(y_{t-1}\right)$ is close to zero, indicating the absence of conditional convergence under these specifications. This outcome, combined with the rejection of the null hypothesis of valid instruments using the J test, points to potential model mis-specification.
The most general specification presented in Column [6] yields a nuanced view of growth dynamics. Both $\mathbf{W} \Delta ln \left(y_{t}\right)$ and $\mathbf{W} \Delta ln \left(y_{t-1}\right)$ remain positive and highly significant, with their combined magnitude summing to $0.82$, which is strictly less than one. This result indicates that growth spillovers, while highly influential, dissipate over larger regional distances and do not propagate indefinitely through the spatial network. The coefficient on the investment rate of a region's neighbours is positive and statistically significant, highlighting substantial spillovers of regional investment activity. Notably, a similar pattern emerges for the coefficient on tertiary educational attainment, suggesting that higher educational levels in neighboring regions contribute positively to a region's economic growth through knowledge diffusion and human capital externalities.
Figure (ref) below depicts a map of European regions, shaded to reflect the magnitude and sign of $\theta_{i}$, the coefficient of the initial level of GDP per capita. Regions shaded in dark blue indicate estimates close to $-1$, while those in dark red represent estimates closer to $1$. Regions shaded in pale blue or pale red reflect estimates near $0$, which are not statistically significant. As shown, most wealthy regions in are dark blue, indicating that higher initial GDP per capita tends to result in slower growth, likely due to diminishing returns to capital. Conversely, dark blue regions in southern Europe, which are relatively poorer, suggest faster growth driven by low starting GDP levels and higher marginal returns to investment. In contrast, regions such as in central Spain and Portugal, in Wales and in northern Romania, exhibit $\theta_{i}$ values close to zero, indicating minimal evidence of conditional regional convergence. Note also that there is considerable heterogeneity in the magnitude of $\theta_{i}$ among blue-shaded regions, reflecting varying speeds of convergence.
Table (ref) below presents Mean Group direct, indirect (spillover), and total effects of $ln\left(inv_{t}\right)$, $ln\left(n_{t}\right)$, $ln\left(educ_{t}\right)$ and $ln\left(sci \& tech_{t}\right)$.
Specifically, by stacking the $N$ observations for each $t$ the model becomes:
where $\Delta \mathbf{y}_{(t)}=\left(\Delta y_{1t},\dots, \Delta y_{Nt}\right)^{\prime}$ is of dimension $N \times 1$, and similarly for the remaining variables, while $\boldsymbol{\Theta} \equiv diag \left(\theta_{1}, \dots, \theta_{N} \right)$, $\boldsymbol{\Psi}_{0} \equiv diag \left(\psi_{0,1}, \dots, \psi_{0,N} \right)$ etc.
Solving the model for $\Delta \mathbf{y}_{(t)}$ yields:
where $L$ denotes the lag operator The matrix of partial derivatives of the expected value of $\Delta \mathbf{y}_{(t)}$ with respect to the $\ell$th covariate is given by:
Building on the framework established by LesagePace2009 and DebarsyEtal2012, the Mean Group direct effect of a unit change in $\mathbf{x}_{\ell (t)}$ on $\Delta \mathbf{y}_{(t)}$ is calculated as the average of the diagonal elements of the matrix in Eq. (ref). The Mean Group indirect effect, on the other hand, is defined as the average of the off-diagonal column sums, capturing the spillover effects of changes in one region's covariate on other regions. The total effect is obtained as the sum of the direct and indirect effects. To analyze heterogeneity, the direct effects for a specific region $i$ correspond to the $(i,i)$th diagonal entry of the matrix in Eq. (ref), while the indirect effects for region $i$ are computed as the sum of the off-diagonal elements in the $i$th row of the same matrix.
The direct effect of $ln\left(inv_{t}\right)$ is small and statistically significant at the $10\%$ level when employing a one-tailed test, indicating that a region's own investment rate has a modest but discernible impact on its economic growth. On the other hand, the indirect effect, capturing the influence of neighboring regions' investment rates, is much larger and also highly significant. This finding is consistent with the conclusions of ElhorstEtal2024, who document that while the local impact of the investment rate tends to be small, its spillover effects are substantial. Overall, a $1\%$ increase in the investment rate is associated with a $0.265\%$ increase in the regional GDP per capita growth rate, with a remarkable $83.4\%$ of this total effect attributable to spillovers from neighboring regions. This substantial contribution of spillovers highlights the pivotal role of cross-regional investment linkages, such as shared infrastructure projects, inter-regional supply chains, and trade networks, in driving economic growth.
For $ln\left(n_{t}\right)$ the direct and spillover effects exhibit a striking contrast in both magnitude and direction. Specifically, a $1\%$ increase in a region's population growth rate corresponds to a $0.531\%$ decrease in its GDP per capita growth rate, reflecting the neoclassical growth theory's prediction of capital dilution effects. Conversely, the same increase in neighboring regions' population growth rates leads to a $0.696\%$ increase in the focal region's growth rate, likely due to enhanced labor market integration, demand spillovers, and knowledge diffusion across regional boundaries. The net result of these opposing forces is a near-cancellation of effects, which highlights an intricate balance of demographic dynamics, where the growth trajectory of a region is shaped not only by its internal factors but also by its interactions within the broader spatial network.
Turning to $ln\left(educ_{t}\right)$ and $ln\left(sci \& tech_{t}\right)$, the total effects for both variables are positive and statistically significant at the $10\%$ level under a one-tailed test.\footnote{Note here that the direct effect of $ln\left(sci \& tech_{t}\right)$ is (small and) positive even if the coefficient on $ln\left(sci \& tech_{t}\right)$ is (small and) negative. This phenomenon arises because a negative coefficient on a spatial lag of a covariate does not necessarily imply that the corresponding spillover effect is also negative. This issue is explained in LesagePace2009 (page 71).} Notably, these total effects are predominantly driven by their spillover components. For instance, the indirect effect of $ln\left(educ_{t}\right)$ is substantial, suggesting that a region benefits significantly from tertiary educational attainment of its neighbors. This outcome aligns with endogenous growth theories that emphasize the role of human capital spillovers in fostering innovation and productivity gains across regions. Similarly, the indirect effect of $ln\left(sci \& tech_{t}\right)$, though smaller in magnitude, indicates that regions appear to derive considerable benefits from innovation and knowledge diffusion originating in adjacent areas, reinforcing the argument for fostering cross-regional collaborations in research and development.
Figure (ref) below displays a map of all European regions, colored to represent the proportion of the total effect of the investment rate, $ln\left(inv_{t}\right)$, on growth attributable to the direct effect. Regions shaded in dark red indicate estimates near zero, suggesting that a substantial portion of the total effect of investment arises from spillover effects from neighboring regions. In contrast, regions shaded in dark blue correspond to estimates near unity, indicating that the majority of the total effect of investment on growth is driven by the region's own (direct) effect. As illustrated, many regions across Europe experience substantial spillover effects relative to the total. Notable exceptions include islands in the Mediterranean, such as Crete, Cyprus and Sicily, where the direct effect of a region's own investment is significantly more pronounced. Similar patterns are observed in northern Scotland, Northern Ireland, Thuringia in Germany (often referred to as “the green heart of Germany” due to its broad, dense forests), and a few capital regions, such as Lisbon in Portugal and Rome in Italy.
Taken together, these findings underscore the critical importance of incorporating both direct and spillover effects in growth models. Ignoring the spatial dimension of regional dynamics risks underestimating the true impact of key growth determinants and overlooking the complex interdependencies that shape economic growth. Moreover, the results highlight the necessity of accounting for slope parameter heterogeneity, as regions often exhibit distinct economic structures and resource endowments. Allowing for such heterogeneity ensures that growth models capture these diverse regional dynamics meaningfully.
From a policy perspective, the results emphasize the necessity of fostering inter-regional cooperation and investment. Policies that enhance connectivity --whether through infrastructure, trade facilitation, or collaborative innovation networks-- are likely to yield substantial aggregate benefits. Moreover, the prominence of spillover effects in education and technology underscores the value of region-wide initiatives to elevate human capital and technological capabilities.
Finally, incorporating interactive fixed effects into growth models offers a robust framework for addressing nonlinear unobserved heterogeneity, latent common shocks, and their region-specific impacts. This is particularly crucial in the context of global and economy-wide events, such as economic crises or technological revolutions, which may affect regions differently.
Table (ref) presents results from a series of robustness checks designed to test the stability of the estimated coefficients under varying assumptions and configurations.
Column [1] replicates the results from Column [6] of Table (ref), serving as the benchmark specification. Columns [2] and [3] retain the same model structure but modify the distance decay parameter $\zeta$ in Eq. (ref) to $0.015$ and $0.01$, respectively. Lowering $\zeta$ results in a slower exponential decay, assigning higher weights to more distant regions. Consequently, the influence of remote regions is amplified, leading to a weaker spatial localization effect. This adjustment reflects a scenario in which regional interactions are less constrained by geographical proximity and more broadly distributed across the spatial network.
In Column [4], the spatial weights matrix is specified as $w_{ij} = \frac{d_{ij}^{-2}}{\sum_{j=1}^{N} d_{ij}^{-2}}$, a specification also explored in Ertur and Koch (2007). This formulation imposes a quadratic decay in spatial interactions, whereby closer regions exert significantly larger weights relative to those farther away. Finally, Column [5] modifies the sample by excluding 10 island regions within the EU. These regions, due to their physical separation, may exhibit weaker integration into the mainland spatial network. This adjustment assesses whether the inclusion of such isolated regions disproportionately influences the estimated spatial dynamics.
Some key observations are as follows: Firstly, despite the variations in the weighting matrix and sample composition, the spatial spillover effects demonstrate remarkable consistency across specifications. In specific, the coefficients of $\mathbf{W} \Delta ln \left(y_{t}\right)$, $\mathbf{W} \Delta ln\left(y_{t-1}\right)$, $\mathbf{W} ln\left(inv_{t}\right)$ and $\mathbf{W} ln\left(educ_{t}\right)$ are similar in magnitude and significance across all configurations. This stability affirms the robustness of the estimates to different choices of the weighting matrix, and reinforces the conclusion that spillovers from neighboring regions play a substantial and persistent role in shaping a region's economic growth trajectory.
In addition, the estimated rate of conditional convergence, captured by the coefficient of $ln \left(y_{t-1}\right)$, also remains consistent across specifications. This invariance suggests that the process of regional income convergence is robust to changes in the spatial weighting structure and sample adjustments.
On the other hand, noticeable variations arise in the impact of a region's own investment rate and population growth: for $ln\left(inv_{t}\right)$, the direct effect appears diminished in some specifications, while its spatial counterpart increases in magnitude. This pattern highlights the dominant role of investment spillovers relative to localized investment impacts, particularly when greater weight is assigned to interactions with distant regions. The coefficient on $ln\left(n_{t}\right)$ also exhibits greater variability across specifications, reflecting the complex and context-dependent nature of demographic dynamics in influencing regional growth.
The robustness checks confirm the reliability of key findings, particularly the significance of spillover effects, while underscoring the need for carefully chosen spatial weighting matrices aligned with the model's theory and context. The consistent convergence rate further validates the framework's ability to capture regional growth dynamics.
This paper develops a Mean Group Instrumental Variables (MGIV) estimator tailored for spatial dynamic panel data models with interactive fixed effects. In contrast to existing methods that assume slope-parameter homogeneity, the MGIV estimator explicitly accommodates heterogeneity across cross-sectional units. Theoretical results establish its consistency and asymptotic normality under large $N$ and $T$ asymptotics. Furthermore, the estimator is asymptotically unbiased, enabling valid inferences without requiring bias correction. Its linear structure also renders it computationally efficient, making it suitable for large-scale applications.