Spatial function-on-function quantile regression

Technical Summary: Spatial Function-on-Function Quantile Regression

Problem Statement
The paper addresses the analysis of spatially indexed functional data, where both response and predictor variables are curves observed at multiple spatial locations. While existing literature offers functional linear models (FLM) and spatial autoregressive (SAR) models, a critical gap exists: current methods primarily focus on conditional means and often rely on Functional Principal Component Analysis (FPCA) for dimension reduction. These mean-based approaches fail to capture distributional dynamics such as tail behavior, heteroskedasticity, and state-dependent volatility, which are crucial in applications like environmental risk assessment (e.g., extreme pollution episodes). Furthermore, standard quantile regression methods for functional data typically assume independent and identically distributed (i.i.d.) curves, ignoring spatial dependence. The specific challenge addressed is modeling the conditional quantile of a functional response $Y_i(t)$ that depends on a functional predictor $X_i(s)$ and a spatially lagged functional response $\tilde{Y}_i(u)$, where the spatial lag introduces endogeneity.

Methodology
The authors propose a Penalized Spatial Function-on-Function Quantile Regression (SFoF-QR) framework. The methodology consists of three main components:

1. Model Specification:
   The model defines the conditional $\tau$-quantile of the response $Y_i(t)$ given the information set $F_i$ (containing the local predictor $X_i$ and the spatially lagged response $\tilde{Y}_i$) as:
   $$Q_\tau\{Y_i(t) | F_i\} = \int_{I_Y} \tilde{Y}_i(u)\rho_\tau(t, u)du + \int_{I_X} X_i(s)\beta_\tau(t, s)ds$$
   Here, $\rho_\tau(t, u)$ represents the spatial spillover surface, and $\beta_\tau(t, s)$ represents the predictor effect surface. Both surfaces are allowed to vary with the quantile level $\tau$. The spatial lag is defined as $\tilde{Y}_i(u) = \sum_{j=1}^n w_{ij}Y_j(u)$, where $W$ is a known spatial weight matrix.

2. Estimation Strategy (Two-Stage IV):
   Because the spatially lagged response $\tilde{Y}_i$ is endogenous (correlated with the structural error), a two-stage instrumental variable (IV) approach is employed:

   • Stage 1: The endogenous spatial lag $\tilde{Y}_i$ is regressed on exogenous instruments generated from the functional predictor $X_i$ and its spatial lags ($X_i, WX_i, W^2X_i, \dots$) using functional quantile regression. This yields instrumented (predicted) lag-score vectors.
   • Stage 2: The estimated coefficient surfaces $\rho_\tau$ and $\beta_\tau$ are obtained by minimizing a penalized quantile loss function using the instrumented lag from Stage 1.

3. Basis Expansion and Penalization:
   To handle the infinite-dimensional nature of the coefficient surfaces, the authors use tensor-product B-spline expansions rather than FPCA truncation. The surfaces are approximated as:
   $$\rho_\tau(t, u) \approx \sum_{\ell,m} \rho^{(\tau)}_{\ell m} \phi_\ell(t)\phi_m(u), \quad \beta_\tau(t, s) \approx \sum_{\ell,k} b^{(\tau)}_{\ell k} \phi_\ell(t)\psi_k(s)$$
   Roughness penalties (integrated squared second derivatives) are applied to the coefficient matrices to ensure smoothness and prevent overfitting. Smoothing parameters are selected via a Bayesian Information Criterion (BIC) grid search.

4. Asymptotic Theory:
   The paper establishes $\sqrt{n}$-asymptotic normality for the spline coefficient estimators under fixed spline dimensions. It further derives the induced finite-rank Gaussian-process limits for the reconstructed coefficient surfaces.

Key Contributions
The paper claims three primary contributions:

1. Novel Model: It introduces the first spatial function-on-function quantile regression model that jointly accounts for spatial correlation via a functional SAR structure while enabling inference on arbitrary conditional quantiles. This captures distributional dynamics beyond the conditional mean.
2. Estimation Framework: It develops a two-stage instrumental-variable estimation strategy that addresses endogeneity without relying on FPCA truncation. By using tensor-product B-splines with roughness penalties, the method avoids the information loss associated with selecting truncation levels and allows for the recovery of non-separable interaction surfaces.
3. Theoretical Foundation: It provides an asymptotic theory for the proposed IV-based estimators, including fixed-spline asymptotic normality and the convergence of the reconstructed surfaces to finite-rank Gaussian processes.

Results

• Monte Carlo Simulations: Extensive simulations demonstrate that the proposed SFoF-QR method (M6) significantly outperforms non-spatial competitors and mean-based spatial competitors (M1–M5) in recovering structural surfaces and predicting responses, particularly under conditions of signal-dependent heteroskedasticity and upper-tail contamination. The advantage of the spatial quantile approach increases with the strength of spatial dependence.
• Empirical Application: The method is applied to PM2.5 air quality data from 1,481 monitoring stations in Italy (2023–2024), using PM10 as the predictor.
  • The SFoF-QR model achieved the lowest in-sample Root Mean Squared Percentage Error (RMSE) of 6.316 and an $R^2$-type agreement of 0.970, outperforming non-spatial quantile models (RMSE 8.556) and matching the penalized spatial mean model.
  • Analysis of the estimated coefficient surfaces revealed pronounced state-dependent heteroscedasticity (stronger separation of quantile curves in winter vs. summer) and localized temporal effects that mean-based models could not capture.
  • The method successfully recovered smooth spatial spillover surfaces, whereas FPCA-based approaches showed less stability.

Significance and Claims
The authors position this work as a robust tool for environmental risk management and complex functional data analysis where tail behavior and spatial dependence are critical. They claim that their approach bridges the gap between spatial functional models and quantile regression, offering a flexible alternative to mean-based models that can handle asymmetric contamination and heteroskedasticity. The implementation is available in the SpatialFoFReg package. The paper modestly notes that while the method performs well in simulations and in-sample analysis, out-of-sample prediction in spatial autoregressive contexts requires specific strategies for handling spatial links between training and test sets, which were addressed via simulation but not fully explored in the real-data application due to data limitations. Future work is suggested regarding data-adaptive weight matrices, computational scalability, and joint estimation to prevent quantile crossing.