Factor-Augmented Spatial Autoregression

• FSAR is a hybrid modeling framework that combines spatial autoregression with low-dimensional factor structures to capture both local spillovers and global shocks.
• It encompasses multiple formulations—functional, panel, and high-dimensional—each using basis expansion, latent shocks, or diversified projections for factor estimation.
• FSAR improves model specification by decoupling spatial dependence via W from common variations, enhancing estimation accuracy and interpretability.

Searching arXiv for recent and foundational papers on FSAR and closely related spatial autoregressive models with factors. Factor-Augmented Spatial Autoregression (FSAR) denotes a class of spatial autoregressive models in which the SAR mechanism is combined with a lower-dimensional factor structure. In current arXiv usage, this combination appears in at least three technically distinct forms: a functional SAR whose basis coefficients induce a factor-like representation, a panel SAR with latent common factors in the disturbance, and a high-dimensional SAR in which latent factors are estimated from residual dependence and then fed back as regressors. Across these formulations, the common objective is to separate dependence transmitted through a spatial weights matrix $W$ from broader common variation captured by factors, loadings, or factor-like scores (Pineda-Ríos et al., 2016, Gefang et al., 25 Oct 2025, Shi et al., 31 Aug 2025).

1. Conceptual scope and model families

FSAR is not a single canonical specification. The literature represented here uses the label, or an equivalent interpretation, for models that retain a spatial autoregressive term while augmenting the model with a low-rank component. The spatial part governs local interdependence through $W$; the factor part captures dependence that is not adequately represented by local spillovers alone.

Formulation	Core representation	Factor role
Functional SAR	$$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$	Basis scores from functional covariates act as factors
Panel SAR with latent common factors	$(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$	Latent common shocks enter the error term
High-dimensional FSAR	$$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$	Estimated latent factors are added as regressors

A useful unifying interpretation is that FSAR augments the SAR mean or disturbance structure with a low-dimensional block. In the functional construction, the block is generated by basis scores of curves. In the panel construction, it is a conventional factor model $A f_t$ in the errors. In the high-dimensional construction, factors are extracted from the residual covariance of multivariate spatial responses and then reintroduced as exogenous covariates. This suggests that FSAR is best understood as a modeling principle rather than a single estimator class (Pineda-Ríos et al., 2016, Gefang et al., 25 Oct 2025, Shi et al., 31 Aug 2025).

2. Structural formulations

The functional SAR model starts from the conventional SAR regression

$$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$

and extends the regressor to a curve $X_i(t)$. Its functional specification is

$$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$

with $\varepsilon \sim N(0,\sigma^2 I)$, $W$0 a symmetric spatial proximity matrix, and $W$1. Solving yields

$W$2

After basis expansion and truncation, the model becomes

$W$3

which is structurally of SAR form (Pineda-Ríos et al., 2016).

The panel formulation with latent common factors is

$W$4

so that

$W$5

At the unit level,

$W$6

This is explicitly described as a factor-augmented spatial autoregression in panel form, with large $W$7, small $W$8, and a known number $W$9 of factors (Gefang et al., 25 Oct 2025).

The high-dimensional formulation begins from componentwise SAR equations

$$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$0

under diverging dimensions $$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$1 and $$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$2. A latent factor structure is imposed on the $$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$3-dimensional error vector: $$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$4 Substitution gives the FSAR system

$$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$5

The defining move is that the fixed-dimensional latent factors are estimated consistently by diversified projections and then fed back into the SAR model as exogenous covariates (Shi et al., 31 Aug 2025).

3. Functional covariates as factor-generating mechanisms

In the functional construction, each spatial unit $$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$6 carries a square-integrable stochastic process $$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$7. With an orthonormal basis $$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$8,

$$Y = \rho W Y + \mathbf{A}\mathbf{b}_k + \varepsilon$$9

and therefore

$(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$0

Truncating at $(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$1 yields

$(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$2

with $(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$3 the matrix of scores and $(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$4 the coefficient vector.

This representation admits an explicit FSAR interpretation. The paper states that the matrix $(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$5 is analogous to factor scores in a factor model or functional principal component scores, while $(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$6 acts as loadings mapping the functional scores to a scalar predictor. Under this reading, one can set $(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$7 and $(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$8, so that the truncated model becomes

$(I_N-W)y_t = X_t\theta + A f_t + \varepsilon_t$9

The paper does not label the model as factor-augmented, but it states that technically the representation is an FSAR where the factor space is generated by the functional data via a basis (Pineda-Ríos et al., 2016).

The factor interpretation is especially transparent when the basis is the Karhunen–Loève basis. Then $$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$0 with uncorrelated scores $$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$1 of variance $$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$2, and the expansion becomes a direct factor decomposition of the functional data. This clarifies a common misconception: factor augmentation in spatial models need not mean latent macro shocks only. It may also arise from deterministic dimension reduction of observed functional covariates. A plausible implication is that FSAR can unify functional regression and spatial econometrics whenever the basis scores are treated as low-rank regressors (Pineda-Ríos et al., 2016).

4. Latent common factors and unrestricted spatial interdependence in panels

The panel FSAR with latent common factors places the factor structure in the disturbance rather than directly in the mean. The additive block $$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$3 represents common shocks that are uncorrelated with $$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$4, while the SAR term $$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$5 captures contemporaneous spatial spillovers. This separation is methodologically important because it distinguishes local network dependence from broad cross-sectional comovement (Gefang et al., 25 Oct 2025).

A distinctive feature of this formulation is the treatment of $$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$6. “Unrestricted spatial interdependence” means that $$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$7 is not fixed ex ante from geographic distance or contiguity rules; off-diagonal entries can be positive or negative; and no parametric structure ties all spatial parameters together. The imposed restrictions are zero diagonal, invertibility of $$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$8, and shrinkage priors encouraging many small entries. The paper requires $$Y_j = \rho_j W Y_j + X_{(j)}\beta_{(j)} + Z b_j + \Omega_j$$9 so that $A f_t$0 is invertible and the Neumann-series expansion is valid. In FSAR terms, this moves the model away from predetermined geography toward data-driven network estimation (Gefang et al., 25 Oct 2025).

The empirical application to 202 NUTS2 regions in the European Union illustrates how the factor and spatial components are interpreted jointly. The model includes a $A f_t$1 unrestricted spatial weights matrix and four known factors: an EU factor, a North factor, a South factor, and an East factor. Estimated coefficients indicate negative initial GVA and lagged GVA growth, consistent with convergence, and positive scientist, capital, employment, and population effects. The estimated $A f_t$2 shows strong country-level clusters, relatively weak German connections to other countries, strong interlinkages among Austria, Belgium, Spain, France, Italy, and Portugal, and a small number of negative spatial relationships suggestive of crowding out. Yet the factor variance decomposition shows that latent factors together explain 22.18% of the residual variance, corresponding to 7.32% of the total variance in GVA growth; most variation is instead explained by observable regressors and spatial spillovers. This directly counters the view that factor augmentation necessarily dominates the dependence structure in FSAR applications (Gefang et al., 25 Oct 2025).

5. High-dimensional FSAR, diversified projections, and penalized estimation

The high-dimensional FSAR model addresses settings in which both the response dimension $A f_t$3 and the covariate dimension $A f_t$4 diverge with $A f_t$5. The paper describes FSAR as a special case of MSAR with a novel factor structure imposed on the high-dimensional random error vector. Because the latent factor dimension $A f_t$6 is fixed, the factors can be estimated consistently by diversified projections as long as the dimension of the multivariate response is diverging (Shi et al., 31 Aug 2025).

The error factor model is

$A f_t$7

with $A f_t$8, $A f_t$9, and idiosyncratic noise $$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$0. The diversified projections method chooses a projection matrix $$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$1 and estimates factors through projected residuals: $$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$2 Under Conditions C1–C9, the factor estimator satisfies

$$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$3

The estimated factors are then added back to each componentwise SAR equation, producing the factor-augmented quasi-log-likelihood

$$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$4

with $$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$5 (Shi et al., 31 Aug 2025).

Variable selection is handled by a SCAD-penalized quasi-likelihood. The penalty derivative is

$$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$6

Theorem 4 states uniform selection consistency when

$$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$7

namely

$$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$8

For tuning, the paper proposes

$$Z(\mathbf{s}) = X(\mathbf{s})\beta + e(\mathbf{s}), \qquad e(\mathbf{s}) = \rho W e(\mathbf{s}) + v,$$9

and Theorem 5 establishes asymptotically correct support recovery for the BIC-selected model (Shi et al., 31 Aug 2025).

The computational motivation is explicit. Conventional MSAR models can require on the order of $X_i(t)$0 parameters, whereas FSAR reduces this to approximately $X_i(t)$1, linear in $X_i(t)$2. Estimation is componentwise and parallelizable. In simulations over DIM, SBM, and LSM networks with $X_i(t)$3, $X_i(t)$4, $X_i(t)$5, and $X_i(t)$6, the factor-augmented maximum likelihood estimator has substantially smaller mean absolute error than CMLE, with relative improvement margins typically 30–50%, and 95% confidence intervals show coverage close to the nominal level, about 94–96% (Shi et al., 31 Aug 2025).

6. Estimation strategies, asymptotics, and interpretation

Despite their different constructions, the FSAR variants share a common estimation logic: spatial dependence and low-rank structure are estimated jointly or sequentially, with the low-rank block either fixed by basis expansion, recovered from residuals, or estimated as a latent error component.

In the functional SAR model, Gaussian likelihood follows from

$X_i(t)$7

Using the truncated approximation $X_i(t)$8, the approximate log-likelihood is

$X_i(t)$9

Because analytical solutions for $$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$0 and $$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$1 are not available, the paper recommends numerical maximization and describes an alternating procedure: initialize $$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$2, estimate $$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$3, update $$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$4 numerically, re-estimate $$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$5, and iterate until convergence. The reconstructed coefficient function is

$$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$6

The paper claims convergence in probability and almost sure convergence of the MLE of $$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$7 and $$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$8 under the stated assumptions (Pineda-Ríos et al., 2016).

In the panel FSAR with unrestricted $$\begin{cases} Y = \displaystyle\int_{\mathcal{T}} X(t)\beta(t)\,dt + \nu,\ \nu = \rho W \nu + \varepsilon, \end{cases}$$9, estimation is two-phase and Bayesian. Phase 1 estimates $\varepsilon \sim N(0,\sigma^2 I)$0 and $\varepsilon \sim N(0,\sigma^2 I)$1 equation by equation by instrumental-variable regression with Dirichlet–Laplace global-local shrinkage priors and Variational Bayes. Phase 2 estimates factors and loadings from residuals

$\varepsilon \sim N(0,\sigma^2 I)$2

using an ordering-invariant prior based on the singular-value decomposition of $\varepsilon \sim N(0,\sigma^2 I)$3,

$\varepsilon \sim N(0,\sigma^2 I)$4

The conditional posteriors of $\varepsilon \sim N(0,\sigma^2 I)$5 and $\varepsilon \sim N(0,\sigma^2 I)$6 are Gaussian, which is central to computational feasibility when $\varepsilon \sim N(0,\sigma^2 I)$7. Reported computation times are approximately 2 seconds per stage for $\varepsilon \sim N(0,\sigma^2 I)$8, about 2 minutes per stage for $\varepsilon \sim N(0,\sigma^2 I)$9, and about 10 minutes per stage for $W$00. In Monte Carlo experiments, similarity between estimated and true $W$01 is summarized by correlation and SSIM; without factors, $W$02 at $W$03, $W$04 at $W$05, and $W$06 at $W$07; with factors, $W$08, $W$09, and $W$10, respectively. Estimated factors correlate with true factors above 0.955 and often around 0.99 (Gefang et al., 25 Oct 2025).

Several interpretation issues follow directly from these results. First, factor augmentation does not have a unique location in the model: it may enter as a functional mean component, an error factor, or an estimated regressor block. Second, the spatial weights matrix need not be exogenous or purely geographic; one recent FSAR construction estimates $W$11 entirely from the data. Third, factor augmentation improves specification even when its variance contribution is modest, because it separates unobserved common shocks from spatial spillovers. A plausible implication is that debates about whether cross-sectional dependence is “spatial” or “factor-driven” are often specification questions about how local and global dependence are partitioned rather than mutually exclusive modeling choices (Pineda-Ríos et al., 2016, Gefang et al., 25 Oct 2025, Shi et al., 31 Aug 2025).

7. Assumptions, limitations, and research directions

The assumptions underlying FSAR differ by formulation but have a common structure. Functional FSAR requires square-integrable covariates $W$12, finite second moments $W$13, and invertibility of $W$14. The panel model with unrestricted $W$15 assumes zero diagonal, $W$16, factors uncorrelated with $W$17, and a known number of factors. The high-dimensional model assumes fixed factor dimension $W$18, bounded loadings, sub-Weibull tails for factors and idiosyncratic errors, sparsity of each $W$19, and growth-rate conditions linking $W$20, $W$21, and $W$22 (Pineda-Ríos et al., 2016, Gefang et al., 25 Oct 2025, Shi et al., 31 Aug 2025).

The limitations are similarly model-specific. The panel Bayesian approach assumes a known number of latent common factors and a static $W$23; the paper notes possible extensions to time-varying $W$24, non-Gaussian errors, endogenous regressors, and nonlinear spatial effects. The high-dimensional framework is restricted to continuous responses, static cross-sectional spatial data, and a fixed number of factors, though the concluding discussion suggests relaxing the continuity assumption, allowing a diverging number of latent factors, modeling cross-response spillover explicitly, extending to dynamic spatial panel data or spatio-temporal factor-augmented models, and exploring alternative penalties such as adaptive LASSO or MCP. In the functional setting, the paper’s asymptotic argument depends on truncation $W$25 and control of the remainder

$W$26

which makes approximation quality sensitive to basis choice and truncation level (Pineda-Ríos et al., 2016, Gefang et al., 25 Oct 2025, Shi et al., 31 Aug 2025).

Taken together, these models establish FSAR as a family of spatial autoregressive specifications enriched by low-dimensional structure. In one branch, factors are functional basis scores derived from observed curves. In another, they are latent common shocks used to isolate unobserved cross-sectional dependence from unrestricted spatial spillovers. In a third, they are consistently estimated from high-dimensional residual dependence and then used to regularize multivariate SAR estimation. The shared message is that spatial autoregression and factor modeling are complementary rather than competing devices for representing dependence, and FSAR provides the formal mechanism for combining them (Pineda-Ríos et al., 2016, Gefang et al., 25 Oct 2025, Shi et al., 31 Aug 2025).