Spherical Spatial Autoregressive Model

• Spherical spatial autoregressive models are advanced frameworks that capture spatial dependencies on spherical domains using spectral harmonics and optimal transport techniques.
• They employ isotropic kernels and LASSO regularization to estimate harmonic coefficients, achieving dimensionality reduction and improved predictive accuracy.
• These models provide rigorous inference, uncertainty quantification, and practical applications in fields such as geochemistry and demography.

Spherical spatial autoregressive models generalize classical spatial autoregression to data residing on spherical domains, such as $\mathbb{S}^2$ or more generally a unit sphere in a separable Hilbert space. These models accommodate rotational invariance and the non-Euclidean geometry intrinsic to spherical data, providing a framework for rigorous spatial modeling, inference, and prediction in geochemistry, demography, and other disciplines. Three principal formulations have emerged: spectral-domain operator-based SPHAR models (Caponera et al., 2019), parametric functional SAR models (Caponera et al., 2021), and optimal transport geometric SSAR models (Xu et al., 23 Jan 2026).

1. Model Specification and Structural Foundations

The canonical order-1 spherical SAR model posits a mean-zero, stationary-in-time Gaussian random field on the sphere, isotropic in space:

$$X_t(\cdot) = \Psi[X_{t-1}(\cdot)] + \epsilon_t(\cdot)$$

where $\Psi:L^2(\mathbb{S}^2)\to L^2(\mathbb{S}^2)$ is a bounded (Hilbert–Schmidt) linear operator or autoregressive kernel, and $\{\epsilon_t\}$ is isotropic spherical white noise with strictly positive angular power spectrum (Caponera et al., 2019, Caponera et al., 2021). $\Psi$ is typically represented as convolution with an isotropic kernel $K(\langle s,y\rangle)$, or equivalently, via expansion in spherical harmonics or Legendre polynomials.

In the optimal-transport based SSAR framework, observations $y_1,\ldots,y_n\in S$ are points on the unit sphere $S=\{\nu\in H:\|\nu\|_H=1\}$ in a Hilbert space $H$, with model residuals defined by the logarithmic Fréchet difference $q_i=y_i\ominus \mu_c$ (where $\mu_c$ is the population Fréchet mean). The spatial dependence is enforced by:

$$q_i - \bar q = \rho\sum_{j=1}^n w_{ij}(q_j - \bar q) + \epsilon_i$$

where $W_n=(w_{ij})$ is a symmetric, row-standardized spatial-weights matrix constructed from spherical neighborhood structure (Xu et al., 23 Jan 2026).

2. Spectral Diagonalization and Parametric Dynamics

Isotropy of the field implies $\Psi$ commutes with rotations, and is thus diagonalizable by spherical harmonics $Y_{\ell m}(s)$. For any $f\in L^2(\mathbb{S}^2)$,

$$\Psi(f)(s) = \sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell \phi_{\ell,m}\langle f, Y_{\ell,m}\rangle_{L^2} Y_{\ell,m}(s)$$

with eigenvalues $\phi_{\ell,m}$ controlling amplification/attenuation of spatial frequencies. For strictly isotropic cases $\phi_{\ell,m}\equiv\phi_\ell$ (Caponera et al., 2019, Caponera et al., 2021). Stationarity and existence of the solution require $\sup_{\ell,m}|\phi_{\ell,m}|<1$.

Parametric extensions enforce regular variation on the spectral parameters, e.g.

$$\phi_\ell=G\ell^{-\alpha},\quad G\in(-1,1)\setminus\{0\},\ \alpha>1,$$

with noise spectrum $C_{\ell;Z}=H\ell^{-\gamma},\ H>0,\ \gamma>2$ (Caponera et al., 2021). The harmonic-domain dynamics decouple into

$$a_{\ell m}(t)=\phi_\ell a_{\ell m}(t-1) + a_{\ell m;Z}(t)$$

and estimation targets $(G,\alpha)$ underlying power-law spectral decay.

3. Estimation Procedures and Oracle Properties

Three estimation regimes prevail. In the spectral SPHAR model, LASSO regularization is imposed on harmonics:

$$\widehat\phi=\arg\min_\phi \left\{ \frac{1}{T}\sum_{t=1}^T \|X_t-\Psi_\phi X_{t-1}\|^2_{L^2} + \lambda \sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell}|\phi_{\ell,m}| \right\}$$

with $\lambda\asymp\sqrt{\log N/T}$, $N=\sum_{\ell=0}^{L-1}(2\ell+1)$. Coordinate descent (soft-thresholding) solves the decoupled scalar AR(1) regressions in the harmonic domain (Caponera et al., 2019).

Parametric models employ nonlinear least-squares, formulating a reduced objective in $\alpha$ by profile minimization in $G$, with closed-form $G^*(\alpha)$ and

$$\widehat{\alpha}_N = \arg\max_{\alpha\in[a_1,a_2]} \left\{ 2\log\widehat{U}_N(\alpha) - \log\widehat{D}_N(\alpha) \right\}$$

yielding consistent and asymptotically normal estimators (Caponera et al., 2021).

The SSAR model adopts a GMM strategy, leveraging moment conditions based on empirical covariances of Hilbert–Schmidt residuals and instruments $P_n$ (typically $P_n=W_n$). The estimate

$$\widehat\rho = \arg\min_\rho \left[m_n(\rho)\right]^2$$

is computed via robust 1D optimization, with asymptotic normality established under broad regularity (Xu et al., 23 Jan 2026).

4. Asymptotic Theory and Statistical Testing

SPHAR(1) LASSO estimators obey nonasymptotic oracle bounds of order $s\sqrt{\log N/T}$ for $s$-sparse coefficient vectors, and analogous prediction error rates. The parametric nonlinear least-square estimator is weakly consistent and asymptotically normal:

$$\sqrt{N}(\widehat{\alpha}_N-\alpha_0) \xrightarrow{d} \mathcal{N}(0,\sigma^2(\theta_0))$$

where $\sigma^2(\theta_0)$ is an explicit function of fourth-order moments of the harmonic coefficients and model parameters (Caponera et al., 2021).

In the SSAR context, Wald-type test statistics

$$T_w = \frac{[n/\eta_n]\widehat\rho^2}{\widehat\sigma_\rho^2}$$

are asymptotically $\chi^2_1$ under the null hypothesis $H_0:\rho=0$, with bootstrap refinement enhancing finite-sample calibration. The GMM estimator is unbiased with RMSE decaying at $n^{-1/2}$ for fixed neighbor number $k$. Laws of large numbers for spatially dependent arrays and delta-method arguments underpin the convergence proofs (Xu et al., 23 Jan 2026).

5. Uncertainty Quantification and Prediction

SSAR models support non-Euclidean, distribution-free, conformal prediction for uncertainty quantification. The split-conformal method partitions data into training/calibration, fits the SSAR, computes conformity scores (residual norms in $C(H)$), and defines prediction sets on $S$ via empirical quantiles. Theoretical guarantees ensure $(1-\alpha)$ coverage up to $o_p(1)$ (Xu et al., 23 Jan 2026). Jackknife+ variants achieve improved, nearly exact finite-sample coverage. The prediction back to the sphere uses the Rodrigues/exponential map to lift operator residuals to points on $S$, maintaining interpretability in terms of geodesic radii.

In the spectral SPHAR context, LASSO-induced sparsity results in prediction using only a subset $\hat S$ of dominating harmonics:

$$\widehat{X}_t(s) = \sum_{(\ell,m)\in\hat S} \widehat\phi_{\ell,m} \alpha_{\ell,m}(t-1) Y_{\ell,m}(s)$$

yielding computational efficiency and improved accuracy when underlying spatial dynamics are low-dimensional (Caponera et al., 2019).

6. Practical Implementation and Applications

Construction of the spatial-weights matrix $W_n$ in SSAR requires explicit consideration of spherical geometry. Neighbors are chosen in the physical coordinate domain; kernels such as uniform or Gaussian, with bandwidth $h$, define weights and are typically row-standardized. In high-dimensional regimes, principal component analysis is recommended for error covariance estimation, though bootstrap is preferable in very large $C(H)$ spaces (Xu et al., 23 Jan 2026).

Real-data applications—Spanish GEMAS geochemical composition and Japanese age-at-death distributions—demonstrate pronounced spatial autoregressive effects ($\rho\approx 0.68$ and $0.41$ respectively). SSAR and SRMSAR yield lower geodesic and compositional errors compared to Euclideanized alternatives, with significantly reduced computation times. In some scenarios, inclusion of exogenous covariates via SRMSAR improves out-of-sample performance; in others, purely spatial models are sufficient (Xu et al., 23 Jan 2026).

7. Interpretative Implications and Model Comparison

A central interpretative feature of spectral SPHAR is harmonic-domain sparsity: only a minority of spherical harmonics contribute to space-time dependence (multipolar selection). This regularization guards against overfitting high-frequency noise and achieves dimensionality reduction ($s\ll N$). The SSAR class leverages optimal transport geometry, extending autoregressive concepts to intrinsic spherical domains, and provides a unified theory for estimation and inference in both finite- and infinite-dimensional settings.

Ignoring spherical geometry in model specification—i.e., by “Euclideanizing” data—induces substantial prediction errors and inefficient inference. Spherical SAR frameworks adapt both kernel construction and spatial regression principles to the manifold context, combining computational feasibility with comprehensive theoretical guarantees (Caponera et al., 2019, Caponera et al., 2021, Xu et al., 23 Jan 2026).