Forecasting Lattice and Point Spatial Data: Comparison of Unilateral and Multilateral Sar Models

1. Introduction

Spatial data are nowaday present in many social and natural phenomena. They are mainly available in two formats: on regular lattices (as digital images and remotely-sensed data), and irregular polygons and sparse points (as in zonings and geological surveys). These data are representable on maps by geographic information systems (GIS), and statical models may be developed for practical purposes; such as filtering, prediction and decision. In digital images the main goals are denoising and sharpening (e.g. Grillenzoni, 2008) to support image segmentation and classification; in polygonal data the aims are structural analyzes, evaluation and out-of-sample prediction (see Barbiero and Grillenzoni, 2022).

One of the main features of spatial data is the auto-correlation (ACR), which arises from the interaction between the spatial units (pixels in images, and centroids in polygons). The general rule is that the minor is the distance of units, the greater is the size of ACR. This dependence must be properly represented in spatial regression models in order to satisfy the basic assumption of independent residuals. This, in turn, is necessary for having unbiased and efficient estimates of the parameters and their standard errors; that is, for performing statistical inference and prediction. It is also possible to show that correlating two independent, but strongly auto-correlated series, then spurious dependence between them may be found (Granger, 2012); this effect also occurs in spatial data.

On the other hand, ACR is useful for forecasting the character under study in areas where it is not measured, especially when exogenous covariates are not available. It follows that representing the ACR in spatial systems is one of the major concerns of regression modelling. As in time series analysis, it may be accomplished by introducing into the equations suitable lagged terms; i.e. the value of the dependent variable in nearby units. This leads to the spatial auto-regressive (SAR) models, which resemble the AR schemes of time series and dynamical systems. However, while in the time there is only a single direction (from past to present), in the plane there are almost infinite directions.

The specification of the models is twofold: unilateral or multilateral, depending on whether the linkage between the nearby units reflects or not a causal ordering. As an example, in lattice data one has the triangular and rook schemes in Figure 1a,b; in the first, the dynamic is unidirectional since starting from the upper-left corner it follows the sequence of writing a text (see Tjøstheim, 1983). Instead, in Figure 1b the dynamic is multi-directional as each cell involves parts of the text that are not yet written (Bustos et al., 2009). Figure 1c,b show the analogous situations in polygonal data, where the position of each unit is defined by its center (geometric or institutional). In Figure 1d the link is multi-directional with four nearest neighbors (NN), whereas in Figure 1c it is in the north direction only.

The statistical consequences of the two specifications for SAR models are important. Multidirection linkages violate the condition of independence between regressors and residuals, therefore make least squares (LS) estimates biased and inconsistent. To solve this problem, many alternative estimators have been proposed, such as maximum likelihood (ML, Anselin, 2003; LeSage and Pace, 2009), generalized method of moment (GM, Kelejian and Prucha, 1999), two-stage least squares (Lee, 2003) and recently indirect inference (Bao et al., 2020). These methods use iterative algorithms of optimization, which may not converge, and matrices of spatial contiguity which are hindered by the system/sample dimension. However, Lee (2002) showed that under suitable conditions of contiguity, the ordinary LS estimator is consistent even for multilateral models.

In this paper we focus on unilateral SAR for both lattice and polygonal data and we show their ability to preserve the optimal properties of LS estimates even with respect to the efficient ML and GM methods. This feature is important because the LS method is linear and can manage dataset of large dimensions, as it may avoid the direct use of contiguity matrices (i.e. the square arrays that connect all spatial units). Further, unilateral models applied to the data sorted in the same spatial direction as the models, have a recursive structure which enable effective forecasting functions. This feature confirms the close connection between parametric identifiability and model predictability.

Spatial prediction of missing data (inside the perimeter of the observed area) has been treated by LeSage and Pace (2004, 2008), Kelejian and Prucha (2007) and Goulard et al. (2017) for multilateral models, and Moijri et al. (2021) in the lattice case. They show that the naive predictor based on the reduced form of SAR models must be corrected with the best linear unbiased predictor (BLUP) of classical statistics. However, the recursive forecasts of unilateral models do not need the BLUP correction and outperform the others. We show evidence of these results with Monte Carlo experiments on synthetic data and out-of-sample forecasting on real data of geosciences, such as digital elevation models (e.g. USGS, 2017) and spatial diffusion of water isotopes (e.g. Moreiras Reynaga et al., 2021).

The paper is organized as follows: Section 2 deals with lattice models, it reviews the conditions of identification and consistency and evaluates their forecasts on random surfaces and digital images; Section 3 deals with SAR models for polygonal and sparse point data, it compares the forecasts of unilateral and multilateral models with various algorithms; two Appendices provide computational details.

2. Regular Lattice Data

Regular lattice data are mostly present in remote sensing and digital images. These data are in the form of rectangular arrays of the type $\mathit{Y}=\left\{y_{ij}\right\}$, with $i=1,2\cdots n$, $j=1,2\cdots m$ the indexes of position, which may be transformed into latitude and longitude. The values $y_{ij}$ are usually autocorrelated (e.g. clustered) and one of the main goals is to filter the array $\mathit{Y}$ with its own values, to obtain interpolates and forecasts: ${\widehat{y}}_{n+h,m+k}$. To accomplish this task, as in time series analysis (Fingleton, 2009), the SAR modeling puts each cell in relation to the contiguous ones, such as $y_{ij}=g\left(y_{i\pm h,j\pm k}\right)+e_{ij}$, where $h,k=1,2\cdots p$ are spatial lags, p is the order of dependence and $\left\{e_{ij}\right\}$ is an unpredictable sequence.

In time series the unidirectional (past-present) ordering of data is also called causal and allows to define the causality between stochastic processes (e.g. Herrera et al., 2012). In lattice data, the causal ordering can be established by following the lexicographic way of reading/writing a text; i.e. processing the cells of $\mathit{Y}$ starting from the upper-left corner. Unilateral dynamics which satisfy this feature are the row-wise$y_{i,j-k}$, the half-plane$y_{i-h,j\pm k},h>0$ and the one-quadrant (or triangular) with $y_{i-h,j-k},(h+k)>0$ in Fig. 1a (see Tjøstheim, 1983, Bustos et al., 2009). Although causality is naturally related to the dynamic of events and aspects of human learning, digital filters of denosing and sharpening may follow non-causal models, such as the rook scheme in Fig. 1b (e.g. Grillenzoni, 2008). However, the unidirectional approach remains the favorite solution as allows properties of sequentiality and prediction to the models.

Under linear and unilateral constraints, the triangular SAR(p) model becomes

$$y_{ij}=\sum _{h=0}^p\sum _{\begin{array}{c}k=0\\ (h+k)>0\end{array}}^p{\varphi }_{hk}{y}_{i-h,j-k}+e_{ij},e_{ij}\sim \mathrm{IN}(0,{\sigma }_e^2),$$

(1)

with independent normal (IN) residuals. In statistics, the model (1) is usually estimated with the maximum likelihood (ML); this requires the vectorization of the data matrix $\mathit{y}={\mathrm{vec}}_j\left(\mathit{Y}\right)$, of length $N=n\times m$ and the inversion of auto-covariance matrices ${\Gamma }_{\mathit{y}}$ of size $N\times N$, see Basu and Reinsel (1993, Sect. 4) and Awang and Shitan (2006). This solution is statistically efficient but computationally expensive and can be implemented only for moderate dimensions of the lattice field.

Instead, the LS method can be applied even for large values of $n,m$; rewriting the model (1) in regression form as

$$\begin{array}{cc}\hfill y_{ij}& ={\mathit{\varphi }}^{\prime }{\mathit{x}}_{ij}+e_{ij},\hfill \\ \hfill {\mathit{x}}_{ij}& ={\left[y_{i,j-1}\cdots y_{i-h,j-k}\cdots y_{i-p,j-p}\right]}^{\prime },\hfill \\ \hfill \mathit{\varphi }& ={\left[{\varphi }_{01}\cdots {\varphi }_{hk}\cdots {\varphi }_{pp}\right]}^{\prime },\hfill \end{array}$$

the LS estimator becomes

$$\begin{array}{cc}\hfill {\widehat{\mathit{\varphi }}}_N& ={\left(\sum _{i=p+1}^n\sum _{j=p+1}^m{\mathit{x}}_{ij}{\mathit{x}}_{ij}^{\prime }\right)}^{-1}\sum _{i=p+1}^n\sum _{j=p+1}^m{\mathit{x}}_{ij}{y}_{ij}\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill & =\mathit{\varphi }+{\left(\sum _{i=p+1}^n\sum _{j=p+1}^m{\mathit{x}}_{ij}{\mathit{x}}_{ij}^{\prime }\right)}^{-1}\sum _{i=p+1}^n\sum _{j=p+1}^m{\mathit{x}}_{ij}{e}_{ij}\hfill \end{array}$$

(2b)

At computational level, the double summation in Eq. (2) can be implemeted in recursive form as follows

$$\begin{array}{cc}\hfill {\mathit{R}}_{i,m}& ={\mathit{R}}_{i-1,m}+{\mathit{P}}_{im},{\mathit{P}}_{ij}={\mathit{P}}_{i,j-1}+{\mathit{x}}_{ij}{\mathit{x}}_{ij}^{\prime },{\mathit{P}}_0={\mathit{O}}_{(p^2-1)},\hfill \\ \hfill {\mathit{s}}_{i,m}& ={\mathit{s}}_{i-1,m}+{\mathit{q}}_{im},{\mathit{q}}_{ij}={\mathit{q}}_{i,j-1}+{\mathit{x}}_{ij}{y}_{ij},{\mathit{q}}_0={\mathit{0}}_{(p^2-1)},\hfill \\ \hfill {\widehat{\mathit{\varphi }}}_N& ={\mathit{R}}_{nm}^{-1}{\mathit{s}}_{nm},{\widehat{\Sigma }}_N={\mathit{R}}_{nm}^{-1}{\widehat{\sigma }}_e^2,\hfill \end{array}$$

with N=$(n-p)(m-p)$. Unlike the ML approach, this algorithm only involves the inversion of the matrix $\mathit{R}$ of size ($p^2$–1), for any values of $n,m$. As regards the statistical properties, the LS estimator (2) is unbiased and consistent under the assumption of unilateral dynamic of the model (1) because E$\left({\mathit{x}}_{ij}{e}_{ij}\right)=\mathit{0}$.

SAR(1). As an example, let p=1, then the model (1) is $y_{ij}=({\varphi }_1{y}_{i,j-1}+{\varphi }_2{y}_{i-1,j}+{\varphi }_3{y}_{i-1,j-1})+e_{ij}$. Applying the formula of Basu and Reinsel (1993, p.634) to the lagged term $y_{i,j-1}$, one obtains the moving average (MA) representation

$$y_{i,j-1}=e_{i,j-1}+\sum _{\begin{array}{c}h=1\\ (h+l)<i\end{array}}^i\sum _{\begin{array}{c}k=1\\ (k+l)<j\end{array}}^{j-1}\sum _{l=0}^{max(i,j-1)}{\psi }_{hkl}{e}_{i-h-l,j-1-k-l},$$

(3)

from which E$\left(y_{i,j-1}{e}_{ij}\right)$=0. Similarly, in Eq. (2b) it follows that E$\left({\mathit{x}}_{ij}{e}_{ij}\right)$=$\mathit{0}$, because all entries of ${\mathit{x}}_{ij}={[y_{i,j-1},y_{i-1,j},y_{i-1,j-1}]}^{\prime }$ do not depend on $e_{ij}$, and the LS estimator (2a) is unbiased. This result does not apply to multilateral SAR models, because the decomposition (3) would involve $e_{i+h,j+k}$; for the rook-type scheme in Fig. 1b the ML estimator is necessary. Finally, as in time series, the consistency of LS does not need the stability of the model (1); rather, it improves with signal-to-noise ratio ${\sigma }_y^2/{\sigma }_e^2$ (e.g. Bhattacharyya et al., 1997).

SAR(p,q,r). Empirical models often have a subset structure, i.e. have missing lagged regressors or they aggregate $y_{i-k,j-h}$ under a common parameter. The proper definition of such models is SAR($p,q,r$), where p=maximum lag of regressors, q=number of spatial units and r=number of parameters. Two parsimonious models that will be extensively applied in the paper are the triangular SAR(1,3,1) and the rook SAR(1,4,1) with drift $\alpha$, defined as

$$\begin{array}{cccc}\hfill y_{ij}& =\alpha +\varphi \left(y_{i-1,j}+y_{i,j-1}+y_{i-1,j-1}\right)/3+e_{ij},\hfill & \hfill y_{i0}& =y_{0j}=0,\hfill \end{array}$$

(4a)

$$\begin{array}{cccc}\hfill y_{ij}& =\alpha +\varphi \left(y_{i,j-1}+y_{i,j+1}+y_{i-1,j}+y_{i+1,j}\right)/4+e_{ij},\hfill & \hfill y_{i,m+1}& =y_{n+1,j}=0,\hfill \end{array}$$

(4b)

see Fig. 1a,b. These models can also be written as $y_{ij}=\alpha +\varphi {\overline{y}}_{ij}+e_{ij}$, where $\overline{y}$ is an average of lagged local terms.

The models (4) are first-order both in the spatial lag and the number of coefficients; their condition of stability is given by $\left|\varphi \right|<1$, as a Corollary of the Proposition 1 of Basu and Reinsel (1993). The rook model (4b) has a multilateral structure, with problems of LS estimation and forecasting; however, it may represent circular dynamics and is a natural competitor of the triangular (4a). This admits an MA representation similar to Eq. (3) with weights

$${\psi }_{hkl}={\displaystyle \frac{(h+k+l)!}{h!k!l!}}{\left({\displaystyle \frac{\varphi }{3}}\right)}^{(h+k+l)},$$

see Bustos et al. (2009). If $e_{ij}\sim \mathrm{IN}(0,{\sigma }_e^2)$, it follows that E$\left(y_{11}\right)=\alpha$; the definition of the border values is crucial in simulating the rook model, as it involves missing values. For this reason, we assume the general starting condition ${\mathit{Y}}_0=\alpha$ which is the simplest out-of-sample forecast.

2.1. The Vector Form

Previous representations are termed raster in GIS softwares; when the lattice dimension is moderate (e.g. $n\times m<1000$), it may be useful to pass at the vector form. Basu and Reinsel (1993), Awang and Shitan (2006) and Mojiri et al. (2021) consider complex forms of vectorization for ML estimation. We follow the spatial econometric approach which vectorizes the data-matrix by columns and build the contiguity matrix of all cells, see LeSage and Pace (2009). Thus let $\mathit{y}={\mathrm{vec}}_j\left(\mathit{Y}\right)$ the $N=n\times m$ data vector; the first-order vector model with drift is

$$\mathit{y}=\alpha {\mathit{1}}_N+\varphi \mathit{W}\mathit{y}+\mathit{e},\mathit{e}\sim \mathrm{N}(\mathit{0},{\sigma }_e^2{\mathit{I}}_N),$$

(5)

where ${\mathit{1}}_N$ is a unit vector of length N and $\mathit{W}$ is an $N\times N$ contiguity matrix. Its structure depends on the spatial dynamics and the marginal values: the models (4a,b) have matrices $\mathit{W}$ with sub-diagonals <1 and sparse 1 on the main diagonal, in correspondence of the border values $y_{i1},y_{1j},y_{im},y_{nj}$. For n=4, m=7 one has the arrays in Figure 2 (see Appendix A); the distinctive feature of the unilateral model (4a) is that $\mathit{W}$ is lower triangular.

The arrays in Figure 2 are obtained by solving the static system $\mathit{y}=\mathit{W}\mathit{y}$, where m determines the number of blocks of size $n\times n$, with matrices ${\mathit{I}}_n$ at the corners (see Appendix A). However, when $\mathit{W}$ is inserted in the dynamic model (5), its diagonal elements $w_{ll}$=1 must be replaced by 0, obtaining ${\mathit{W}}_0$. Each sub-diagonal corresponds to a specific regressor $y_{i-h,j-k}$ of models (4); thus, the array in Fig. 2a can be decomposed as ${\mathit{W}}_0=(W_1+W_2+W_3)$, providing a SAR(3) model. The data-generation of $\left\{y_l\right\}$ proceeds by rows, starting from an initial $\mathit{y}={\mathit{y}}_0$, as

$$y_l=\alpha +\varphi {\mathit{w}}_l^{\prime }\mathit{y}+e_l,l=1,2\cdots N,$$

(6)

where ${\mathit{w}}_l^{\prime }$ is the l-th row of ${\mathit{W}}_0$. Unilateral processes are insensitive to ${\mathit{y}}_0$, since ${\mathit{W}}_0$ triangular fills $\mathit{y}$ recursively; instead, the model (4b), having the forward terms $y_{i,j+1},y_{i+1,j}$, needs a non-null initial vector, e.g. ${\mathit{y}}_0=\alpha {\mathit{1}}_N$.

Using matrix algebra, one may obtain the reduced (MA) form of the system (5)

$$\begin{array}{cc}\hfill \left({\mathit{I}}_N-\varphi \mathit{W}\right)\mathit{y}& =\alpha {\mathit{1}}_N+\mathit{e},\hfill \\ \hfill \mathit{y}& ={\left({\mathit{I}}_N-\varphi \mathit{W}\right)}^{-1}\left(\alpha {\mathit{1}}_N+\mathit{e}\right),\hfill \end{array}$$

(7)

this provides an automatic way to generate SAR data, independent of ${\mathit{y}}_0$. Differences with the AR method (6) arise from the choice of the border values and the contiguity matrices; however, in models with $\mathit{W}$ triangular and diag($\mathit{W}$)=0 the algorithms (6) and (7) are equivalent. In this case, the impulse-response matrix $\left({\mathit{I}}_N-\varphi \mathit{W}\right)$ is always invertible and triangular, and leads to the MA decomposition

$$\begin{array}{cc}\hfill \mathit{M}& ={\left({\mathit{I}}_N-\varphi \mathit{W}\right)}^{-1}={\mathit{I}}_N+\sum _{k=1}^{\infty }{\varphi }^k{\mathit{W}}^k,\hfill \\ \hfill \mathit{y}& ={\displaystyle \frac{\alpha }{(1-\varphi )}}{\mathit{1}}_N+\sum _{k=1}^{\infty }{\varphi }^k{\mathit{W}}^k\mathit{e}+\mathit{e},\hfill \end{array}$$

(8)

which converges when $\left|\varphi \right|<1$. The expression of the constant term follows from the row-normalization, where $\mathit{W}{\mathit{1}}_N={\mathit{1}}_N$ and also ${\mathit{W}}^k{\mathit{1}}_N={\mathit{1}}_N$. Further, each element of $\mathit{y}$ becomes a linear combination of $\mathit{e}$, and in unidirectional models, they are function of one-quadrant errors only, providing the MA decomposition

$$y_{ij}={\displaystyle \frac{\alpha }{(1-\varphi )}}+\sum _{h=1}^{i-1}\sum _{k=1}^{j-1}{\psi }_{hk}{e}_{i-h,j-k}+e_{ij},$$

where the weights ${\psi }_{hk}\to 0$ are complex functions of ${\varphi }^k{\mathit{W}}^k$.

Finally, the representation (5) enables to define the LS estimator of parameters $\mathit{\delta }={[\alpha ,\varphi ]}^{\prime }$ of the models (4) using the entire sample $\mathit{Z}=\left[{\mathit{1}}_N,\mathit{W}\mathit{y}\right]$ as

$$\begin{array}{cc}\hfill {\widehat{\mathit{\delta }}}_N& ={\left(\mathit{Z}{\mathit{Z}}^{\prime }\right)}^{-1}\mathit{Z}\mathit{y},\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill {\widehat{\Sigma }}_N& ={\left(\mathit{Z}{\mathit{Z}}^{\prime }\right)}^{-1}{\widehat{\sigma }}_e^2,\hfill \end{array}$$

(9b)

this improves the estimator (2) based on $N-p$ observations. The matrix ${\widehat{\Sigma }}_N$ is the dispersion and if $\mathit{W}$ is triangular, as in Fig. 2a, then LS (9) is consistent; its analysis will be developed in Sect. 3.

Spatial Forecasting. Prediction is one of the central goals of SAR modeling; typical examples are restoring parts of remote sensing data which are hidden by clouds or extrapolating the image out of its perimeter. However, existing literature on lattice data has mostly concerned with in-sample filtering and interpolation, focusing on techniques of image sharpening (Grillenzoni, 2008), robust denosing (Bustos et al., 2009) and trend estimation (Mojiri et al., 2021), also in conjunction with non-parametric smoothing.

This section aims to forecast data segments which are external of the estimation perimeter $n\times m$; e.g. on the right hand side as $\left\{y_{i,m+k}\right\}$ or, in vector form ${\mathit{y}}_{N+l}$ on units placed beyond N. By, defining $h,k,l=1,2\cdots H,K,L$ the prediction indexes, the forecast function depends on the SAR representation used, see Kelejian and Prucha (2007). Here, we have the lattice form (4), the vector AR (6) and the reduced MA (7); for the triangular model (4a), the 3 predictors are

$$\begin{array}{cc}\hfill {\widehat{y}}_{n+h,m+k}& =\mathrm{E}\left(y_{n+h,m+k}|y_{n-i+1,m-j+1}\right),h,k=1,2\cdots H,K,\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\widehat{y}}_{n+h,m+k}& =\alpha +\varphi \left({\widehat{y}}_{n+h,m+k-1}+{\widehat{y}}_{n+h-1,m+k-1}+{\widehat{y}}_{n+h-1,m+k}\right)/3,\hfill \end{array}$$

(10a)

$$\begin{array}{cc}\hfill {\widehat{y}}_{N+l}& =\alpha +\varphi {\mathit{w}}_{N+l}^{\prime }\widehat{\mathit{y}},{\widehat{\mathit{y}}}^{\prime }=\left[y_{11},y_{21}\cdots y_{nm},{\widehat{y}}_{N+1}\cdots {\widehat{y}}_{N+l-1}\cdots 0\right],\hfill \end{array}$$

(10b)

$$\begin{array}{cc}\hfill \widehat{\mathit{y}}& ={\left({\mathit{I}}_{N+L}-\varphi {\mathit{W}}_{N+L}\right)}^{-1}\alpha {\mathit{1}}_{N+L},l=1,2\cdots L,\hfill \end{array}$$

(10c)

where ${\mathit{w}}_{N+l}^{\prime }$ is the $(N+l)$-th row of the augmented matrix ${\mathit{W}}_{N+L}$ and the vector $\widehat{\mathit{y}}$ in (10b) is updated with the forecasts. The predictor (10c) is nearly automatic, but in absence of exogenous variables provides constant forecasts; on the other hand, the functions (10a,b) need border values to start. These can be set as ${\widehat{y}}_{1,m+k}^{\left(0\right)}=\widehat{\alpha }$ for all k, and may be improved by r-iterating the forecasts with the nearby values as ${\widehat{y}}_{1,m+k}^{(r+1)}={\widehat{y}}_{2,m+k-1}^{\left(r\right)}$, etc.. This approach can be extensively applied to the rook model (4b), to provide initial values also of its forward elements $y_{ij+1}$.

The variance of predictors is useful for building confidence intervals of the forecasts (10) and testing hypotheses; as in time series models, the process variance ${\sigma }_y^2=\mathrm{var}\left(y_{ij}\right)$ represent an upper bound. Applying the formula of Mojiri et al. (2018, p.187) to the triangular model (4a), one obtains

$${\sigma }_y^2={\sigma }_e^2/\left\{1-2{\displaystyle \frac{{\varphi }^2}{9}}\left[{\displaystyle \frac{3}{2}}+\rho (1,1)+\rho (0,1)+\rho (1,0)\right]\right\}\le {\sigma }_e^2/\left(1-{\varphi }^2\right),$$

where $\left|\rho \right(.,.\left)\right|<1$ are autocorrelation coefficients. The last formula is the variance of an AR(1) process; hence, for the model (4a), the approximate variance of the predictor (10a) may be ${\sigma }_e^2\left[1-{\varphi }^{2(h+k)}\right]/\left(1-{\varphi }^2\right)$. A more general result can be obtained for the representations (5) and (10c), namely

$$\mathrm{V}\left(\widehat{\mathit{y}}\right)={\left({\mathit{I}}_{N+L}-\varphi {\mathit{W}}_{N+L}\right)}^{-1}{\left({\mathit{I}}_{N+L}-\varphi {\mathit{W}}_{N+L}^{\prime }\right)}^{-1}{\sigma }_e^2,$$

whose diagonal elements increase with the index l.

2.2. Simulations and Applications

We develop simulation experiments to test the performance of LS (9) and ML (19) estimators (which are based on the same matrices $\mathit{W}$), applied to unilateral and multilateral SAR models. We consider the models (4a,b) with parameters ${\alpha }_0=1$, ${\varphi }_0=0.5,0.75,1.0,1.025$ and border values $(\alpha +e_{ij})$, i.e. $y_{00}=0$ etc.. Owing to the forward terms $y_{i+1,j},y_{i,j+1}$, the generation of the rook process (4b) is not fully possible with the AR Equation (6) and needs the MA reduced form (7), with the matrix in Fig. 2b, without the diagonal entries. Instead, for the unilateral model (4a), the generating methods (6) and (7) are equivalent.

200 replications on a $32\times 32$ lattice (N=1024 cells) are generated with Normal disturbances and ${\sigma }_e^2=1$. Performance statistics are the relative bias $|\widehat{\varphi }-{\varphi }_0|/{\varphi }_0$, the relative root mean squared errors (the square root of MSE($\widehat{\varphi }$)/${\varphi }_0^2$) and the p-values of the Jarque and Bera (1987) Normality test. Since the relative statistics of $\widehat{\alpha },\widehat{\varphi }$ do not depend on the parameter size, the overall performance of LS and ML methods can be compared by summing the index values.

The results are displayed in Table 1 and Figure 3; the main findings are as follows: LS (9) is uniformly better than ML (19) when the matrix $\mathit{W}$ is triangular. Its efficiency improves as ${\varphi }_0\to 1$, as a consequence of the super-consistency property of LS estimates in AR models (e.g. Grillenzoni, 2004). Its performance slightly worsens when using LS (2), due to the smaller sample size, which excludes border values. The performance of LS estimates in multilateral SAR models is biased for $\varphi <1$, but still benefits of the super-consistency. It is interesting to note that, unlike time series, the Normal distribution of LS estimates of the AR parameter ${\widehat{\varphi }}_N$ holds even in the presence of the unit-root${\varphi }_0=1$ (see Figure 3b).

As regard the ML method (19), we use the implementation of LeSage and Pace (2009), which is computationally demanding and time consuming. However, it significantly improves the estimates of multilateral models, providing good levels of unbiasedness and efficiency when ${\varphi }_0\ll 1$. On the other hand, it requires conditions of stationarity and its performance in triangular models is negative, meaning that ML is mainly designed for sparse weight matrices. We will resume this topic in Sect. 3 for non-lattice data; as a conclusion, we can state that LS is suitable for unilateral models, whereas ML must be used in multilateral SAR.

We also carry out two empirical applications to test the performance of various SAR models in out-of-sample forecasting of spatial data. In comparing predictions, one cannot proceed by simulations as in the experiments of Tab. 1, because simulated data depend on the ground models, which will influence the final results. Hence, we fit SAR models to data generated from external sources.

Random Surface. The first application considers rough random surfaces (Sheffield, 2022); these techniques are used in physics to study electromagnetic, fluid and plasma phenomena (see Sjodin, 2017), and are applied in engineering (e.g. 3D printers). We follow an approach based on trigonometric series, where a matrix of random numbers is convolved with a Gaussian filter to achieve spatial autocorrelation. Thus, let the fast Fourier transform (FFT) of the bivariate random sequence $u_{ij}$ on a $m\times n$ regular lattice be

$$F\left(u_{ij}\right)=\sum _{h=0}^{n-1}\sum _{k=0}^{m-1}{u}_{hk}{e}^{-ı2\pi hi/n}{e}^{-ı2\pi kj/m},u_{ij}\sim \mathrm{IN}(0,{\sigma }_u^2),$$

where ı is the imaginary unit. Let $R_{ij}\left(\delta \right)=e^{-0.5(i^2+j^2)/{\delta }^2}$ an exponential ACR function with $\delta >0$; then the model surface generator becomes

$$y_{ij}=\alpha +\Re \left\{F^{-1}\left[R_{ij}\left(\delta \right)F\left(u_{ij}\right)\right]\right\}+e_{ij},e_{ij}\sim \mathrm{IN}(0,{\sigma }_e^2),$$

(11)

where $\Re (·)$ is the real part, $F^{-1}(·)$ is the inverse FFT and $e_{ij}$ is a white noise.

Despite the randomness of $u_{hk}$, the rough surface $g_{ij}=\Re (·)$ in Eq. (11) is still smooth, hence we add an external noise $e_{ij}$ to increase its realism; the resulting series $y_{ij}$ resemble real terrain models of geostatistics (see Figs. 4a,b). In the experiment we set n=m=32 (which yields N=1024), $e,u$ independent Gaussian with ${\sigma }_u$=1, ${\sigma }_e$=0.067, and $\alpha$=1, $\delta$=2.5 provide mild spatial ACR. The goal is to forecast the last K=7 columns (about 22%) of the surface $\mathit{Y}$ with the algorithms (10). Table 2 reports the in-sample parameter estimates (with LS and ML methods) of the models (4a,b) and the mean absolute percentage forecast errors statistics

$$\mathrm{MAPE}={\displaystyle \frac{1}{nK}}\sum _{i=1}^n\sum _{j=m-K+1}^m\left|{\displaystyle \frac{{\widehat{y}}_{ij}-y_{ij}}{y_{ij}}}\right|.$$

(12)

A triangular SAR(1,3,3) model (with 3 distinct AR parameters) is also fitted to the data, but it does not improve the forecasts of the constrained (4a). The rook (4b) performs poorly with the (biased) LS estimates, whereas it significantly improves with the ML ones; anyway, the path of the predicted surface in Figure 4d remains disappointing. The performance of the vector predictor (10b) is similar to that of the lattice method (10a) and is reported in Appendix B, together with that of the function (10c).

Remote Sensing. The second application regards a real case study; we consider high-resolution elevation data obtained with aerial laser scanners. The sample array $\mathit{Y}$ comes from USGS (2017) and has 340×455 pixels (see Fig. 5); this would yield a vector SAR model with 154,700 rows. Further, the image extrapolation may require the inclusion of spatial trend components (e.g. Grillenzoni, 2004), such as bivariate polynomials $g_d(i,j)$ of degree $d\ge 1$. Given the implied computational issues, we only perform analysis with the lattice implementation (2); in the case of a triangular model with quadratic trend, the forecasting function becomes

$$\begin{array}{cc}\hfill {\widehat{y}}_{i+h,j+k}& =\alpha +g_2(i+h,j+k)+\varphi \left({\widehat{y}}_{i+h-1,j+k-1}+{\widehat{y}}_{i+h-1,j+k}+{\widehat{y}}_{i+h,j+k-1}\right)/3,\hfill \\ \hfill g_2(i+h,j+k)& ={\beta }_1(i+h)+{\beta }_2(j+k)+{\beta }_3{(i+h)}^2+{\beta }_4{(j+k)}^2+{\beta }_5(i+h)(j+k),\hfill \end{array}$$

with $h,k=1,2\cdots H,K$, which can be sequentially managed by cells. The AR part of the rook model involves forward values in the lower-right side; these may be provided by the polynomial itself: ${\tilde{y}}_{i+h+1,j+k+1}={\widehat{g}}_2(i+h+1,j+k+1)$.

In the application to the USGS image, we left the right portion of 55 columns (about 12 % of the array) as out-of-sample data to forecast. LS estimates of the parameters are reported in Table 3 and graphical results are displayed in Figure 5. Despite the better in-sample fitting (with $R^2$=0.98), the rook model has a disappointing forecasting performance both in terms of MAPE statistic and visual appearance (Figure 5b). The reason is partly due to the (biased) estimates of its parameters, with $\widehat{\alpha }<0,\widehat{\varphi }>1$; however, adjusting their values reduces the MAPE statistics but do not change the (uniform) path of the predicted surface.

In real life applications, where small and inner portions of images must be restored, the spatial polynomials $g(i,j)$ may be fitted locally (around missing pixels) or may be replaced by nonparametric smoothers. In any event, the role of the AR components remains fundamental, as the coefficient $\widehat{\varphi }$ is the most significant in both models of Table , with statistic $T_{\widehat{\varphi }}>100$.

3. Point and Polygonal Data

When the data are irregularly distributed in space (either as point processes or polygonal areas), notations and representations have to change. The spatial process is now defined as $\left\{y_s\right\},s=1,2,\cdots N$; where the index s also applies to the planar coordinates ${\mathit{x}}_s^{\prime }=[i_s,j_s]$ (latitude and longitude) of the polygon centers (geometric or political), see Figure 1c,d. These variables are irregularly distributed in the plane but have a fixed (non-stochastic) nature; moreover, they may influence the level of the process $y_s$, according to spatial trends. Thus, the SAR(1) model becomes

$$\begin{array}{cc}\hfill y_s& =\alpha +\varphi y_{s-1}+{\mathit{\beta }}^{\prime }{\mathit{x}}_s+e_s,s=1,2\cdots N,\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill \mathit{y}& =\alpha {\mathit{1}}_N+\varphi \mathit{W}\mathit{y}+\mathit{X}\mathit{\beta }+\mathit{e},\mathit{e}\sim \mathrm{N}\left(\mathit{0},{\sigma }_e^2{\mathit{I}}_N\right),\hfill \end{array}$$

(13b)

$$\begin{array}{cc}\hfill & ={\left({\mathit{I}}_N-\varphi \mathit{W}\right)}^{-1}\left(\alpha {\mathit{1}}_N+\mathit{X}\mathit{\beta }+\mathit{e}\right),\hfill \end{array}$$

(13c)

where $s-1$ refers to the unit which is closest to the s-th term, and the vector of regressors ${\mathit{x}}_s^{\prime }=[i_s,j_s,\cdots ]$ may include exogenous covariates.

Unlike the previous section, the adjacency matrix $\mathit{W}$ has an irregular structure, which depends on the rule of contiguity. The most common rule is to put each observation $y_s$ in relation to its nearest neighbor (NN) term $y_r$, according to the Euclidean distance: ${min}_r\{d_{sr}={[{(i_s-i_r)}^2+{(j_s-j_r)}^2]}^{1/2}\}$. Furthermore, under the unilateral constraint, the north-west (NW) direction may be followed as in the lattice case. However, polygonal data have not a lexicographic order; hence, the unilateral constraint may simply be defined along the north direction. In this setting, the observations $\{y_s,{\mathit{x}}_s\}$ are ordered according to the shortest distance from the northern border, i.e. according to the inverse of the latitude $1/i_s$. Such ordering may be denoted by $\left(i_s\right)$ or simply $\left(s\right)$, and each equation of the model (13a) can be written in sequential form, as a time series model

$$y_{\left(s\right)}=\alpha +\varphi y_{(s-1)}+{\mathit{\beta }}^{\prime }{\mathit{x}}_{\left(s\right)}+e_{\left(s\right)},y_{\left(0\right)}=0,$$

(14)

where $y_{(s-1)}$ is the north NN of $y_{\left(s\right)}$. The property E$\left[y_{(s-1)}{e}_{\left(s\right)}\right]=0$ as in time series is the basis for unbiased LS estimation and forecasting. However, it does not hold in the simple NN linkage, even when data are ordered along the latitude.

As an illustration of NN matrices, we simulate N=30 random points with $[i_s,j_s]\sim \mathrm{IU}{(0,1)}^2$ and consider q=3 contiguous terms, see Figure f6. Under the north ordering of data, the entries of matrices are concentrated around the null main diagonal diag($\mathit{W}$)=$\mathit{0}$, and with the unilateral constraint the array is lower triangular. In order to obtain non-sparse matrices, LeSage and Pace (2004) ordered data according to the sum $(i_s+j_s)$, but just used the simple NN rule. The matrices ${\mathit{W}}_k$ in Figure f6 yield SAR(3) models, but if they are combined as $\overline{\mathit{W}}=3^{-1}{\sum }_{k=1}^3{\mathit{W}}_k$ they provide constrained SAR(3,1). The analogous of the lattice model (4a) then becomes

$$y_{\left(s\right)}=\alpha +\varphi \left(q^{-1}\sum _{k=1}^q{y}_{(s-k)}\right)+{\mathit{\beta }}^{\prime }{\mathit{x}}_{\left(s\right)}+e_{\left(s\right)}.$$

(15)

The econometric literature often discusses models with autocorrelated errors, such as $\mathit{u}={\left({\mathit{I}}_N-\theta \mathit{W}\right)}^{-1}\mathit{e}$, e.g. Kelejian and Prucha (2007). Apart from estimation difficulties, it should be noted that they are encompassed by unconstrained SAR(2); indeed, the combined system $\left({\mathit{I}}_N-\theta \mathit{W}\right)\left({\mathit{I}}_N-\varphi \mathit{W}\right)\mathit{y}=\left({\mathit{I}}_N-\theta \mathit{W}\right)\left(\alpha +\mathit{X}\mathit{\beta }\right)+\mathit{e}$ involves just 2 spatial lags. In place of $\mathit{u}$, it is better to insert into the model (13b) a lagged term on the exogenous part $+\mathit{W}\mathit{X}\mathit{\gamma }$, that corresponds to $\mathit{\gamma }=\theta \mathit{\beta }$ in the model with $\mathit{u}$ errors. This is the Durbin model (LeSage and Pace, 2009), which is useful when ${\mathit{x}}_{\left(s\right)}$ are autocorrelated, i.e. $(i_s,j_s)$ have a non-random pattern.

When the matrix $\mathit{W}$ is lower triangular, the transfer function $({\mathit{I}}_N-\varphi \mathit{W})$ is always nonsingular and its inverse admits a convergent power expansion

$$\mathit{y}={\displaystyle \frac{\alpha }{(1-\varphi )}}{\mathit{1}}_N+\mathit{X}\mathit{\beta }+\sum _{k=1}^{\infty }{\varphi }^k{\mathit{W}}^k\mathit{X}\mathit{\beta }+\sum _{k=1}^{\infty }{\varphi }^k{\mathit{W}}^k\mathit{e}+\mathit{e}.$$

This expression shows that the process $y_{\left(s\right)}$, in the location s-th, depends on the values of the inputs ${\mathit{x}}_{\left(r\right)}$ and shocks $e_{\left(r\right)}$ in the north locations ($r<s$), and with weights ${\varphi }^k{\mathit{W}}^k$ which decline with the distance. To appreciate this diffusion effect, note that when $\mathit{W}$ is strictly lower triangular, then also ${\mathit{W}}^2$ is lower triangular but with nonzero entries shifted away from the main diagonal. Hence, ${\mathit{W}}^2$ reflects the second-order contiguity (the neighbors of the neighbors), and so forth for ${\mathit{W}}^k\to \mathit{O}$ the null matrix. The same feature, however, does not hold in non-triangular matrices, such as those based on the simple NN.

Figure 6. Contiguity matrices and planar graphs of N=30 random centroids on the unit square [0,1]. Panels: (a,c) North nearest neighbors (NN); (b,d) Multidirection NN. Matrices: ${\mathit{W}}_1$ first NN (blue); ${\mathit{W}}_2$ second NN (red); ${\mathit{W}}_3$ third NN (green).

Figure 6. Contiguity matrices and planar graphs of N=30 random centroids on the unit square [0,1]. Panels: (a,c) North nearest neighbors (NN); (b,d) Multidirection NN. Matrices: ${\mathit{W}}_1$ first NN (blue); ${\mathit{W}}_2$ second NN (red); ${\mathit{W}}_3$ third NN (green).

The unilateral model (14) corresponds to a time-series AR with unequally spaced observations. This feature is more clear in the presence of multiple lagged terms $y_{(s-k)}$; the suitable treatment is to weight the observations by the inverse Euclidean distance $1/d_{(s,s-k)}$. The analog of the constrained model (15), with north NN linkage and with row normalization, is then given by

$$\begin{array}{cc}\hfill y_{\left(s\right)}& =\alpha +\varphi \left(\sum _{k=1}^q{w}_{\left(sk\right)}{y}_{(s-k)}\right)+{\mathit{\beta }}^{\prime }{\mathit{x}}_{\left(s\right)}+e_{\left(s\right)},\hfill \\ \hfill w_{\left(sk\right)}& =d_{(s,s-k)}^{-1}/\sum _{k=1}^q{d}_{(s,s-k)}^{-1},d_{(s,s-k)}={\left[{\left(i_{\left(s\right)}-i_{(s-k)}\right)}^2+{\left(j_{\left(s\right)}-j_{(s-k)}\right)}^2\right]}^{1/2},\hfill \end{array}$$

(16)

this also admits a vector representation as (13b) with triangular matrix $\mathit{W}$, having q decreasing coefficients per row.

Alternatively, since the inverse distances decay too fast, one can use the exponential weights $w_{\left(sk\right)}\propto {\lambda }^{d_{(s,s-k)}}$, with 0$<\lambda <$1. In vector form, one has to define the diagonal matrices ${\Lambda }_k={\lambda }^{{\mathit{D}}_k}$ where ${\mathit{D}}_k=\mathrm{diag}[0\cdots 0,d_{(k,1)},$$d_{(k+1,2)}\cdots d_{(N,N-k)}]$ contain the k-th north NN distances of the centroids. These arrays must be multiplied by the k-lag contiguity matrices ${\mathit{W}}_k$, then summed as $\mathit{W}={\sum }_k{\Lambda }_k{\mathit{W}}_k$ and finally normalized by rows. The resulting model is nonlinear in the parameters and requires iterative algorithms; however, the value of $\lambda$ may be selected a-priori in the range [.9, 1), then the estimators compute the suitable value of $\varphi$.

3.1. Estimation Methods

Writing the models (14)-(16) in vector form, as $y_{\left(s\right)}={\mathit{\delta }}^{\prime }{\mathit{z}}_{\left(s\right)}+e_{\left(s\right)}$, with $\mathit{\delta }={\left[\alpha ,\varphi ,\mathit{\beta }\right]}^{\prime }$ and ${\mathit{z}}_{\left(s\right)}={\left[1,{\overline{y}}_{(s-1)},{\mathit{x}}_{\left(s\right)}\right]}^{\prime }$, where ${\overline{y}}_{(s-1)}$ is the average of q-north NN terms, then the LS estimator of the parameters $\mathit{\delta }$ is given by

$$\begin{array}{cc}\hfill {\widehat{\mathit{\delta }}}_N& ={\left(\sum _{s=q+1}^N{\mathit{z}}_{\left(s\right)}{\mathit{z}}_{\left(s\right)}^{\prime }\right)}^{-1}\sum _{s=q+1}^N{\mathit{z}}_{\left(s\right)}{y}_{\left(s\right)},\hfill \\ \hfill & =\mathit{\delta }+{\left(\sum _{s=q+1}^N{\mathit{z}}_{\left(s\right)}{\mathit{z}}_{\left(s\right)}^{\prime }\right)}^{-1}\sum _{s=q+1}^N{\mathit{z}}_{\left(s\right)}{e}_{\left(s\right)},\hfill \end{array}$$

(17)

its unbiasedness follows from E$\left[{\mathit{z}}_{\left(s\right)}{e}_{\left(s\right)}\right]=\mathit{0}$ by the sequential structure of the models, in particular E$\left[e_{\left(s\right)}{\overline{y}}_{(s-1)}\right]=0$ in the unilateral model.

To see what happens in detail, consider the LS estimator in matrix form

$${\widehat{\mathit{\delta }}}_N={\left({\mathit{Z}}^{\prime }\mathit{Z}\right)}^{-1}{\mathit{Z}}^{\prime }\mathit{y}=\mathit{\delta }+{\left({\mathit{Z}}^{\prime }\mathit{Z}\right)}^{-1}{\mathit{Z}}^{\prime }\mathit{e}$$

(18)

with regressors $\mathit{Z}=\left[{\mathit{1}}_N,\mathit{W}\mathit{y},\mathit{X}\right]$, the crucial term is $\mathit{W}\mathit{y}$. From Eq. (13c)

$$\begin{array}{cc}\hfill \mathrm{E}\left({\mathit{Z}}^{\prime }\mathit{e}\right)& =\mathrm{E}\left[{\left(\mathit{W}\mathit{y}\right)}^{\prime }\mathit{e}\right]=\mathrm{E}\left\{{\left[{({\mathit{I}}_N-\varphi \mathit{W})}^{-1}(\alpha {\mathit{1}}_N+\mathit{X}\mathit{\beta }+\mathit{e})\right]}^{\prime }{\mathit{W}}^{\prime }\mathit{e}\right\}\hfill \\ \hfill & =\mathrm{E}\left[{\mathit{e}}^{\prime }{({\mathit{I}}_N-\varphi {\mathit{W}}^{\prime })}^{-1}{\mathit{W}}^{\prime }\mathit{e}\right]=\mathrm{E}\left({\mathit{e}}^{\prime }{\mathit{H}}^{\prime }\mathit{e}\right),\left(\mathrm{say}\right)\hfill \end{array}$$

this term vanishes only if the trace (tr) of ${\mathit{H}}^{\prime }={\left({\mathit{I}}_N-\varphi {\mathit{W}}^{\prime }\right)}^{-1}{\mathit{W}}^{\prime }$ is 0, because

$$\mathrm{E}\left({\mathit{y}}^{\prime }{\mathit{W}}^{\prime }\mathit{e}\right)=\mathrm{E}\left[\mathrm{tr}\left({\mathit{e}}^{\prime }{\mathit{H}}^{\prime }\mathit{e}\right)\right]=\mathrm{tr}\left[\mathit{H}\mathrm{E}\left(\mathit{e}{\mathit{e}}^{\prime }\right)\right]=\mathrm{tr}\left(\mathit{H}\right){\sigma }_e^2{\mathit{I}}_N.$$

Now, $\mathrm{tr}\left(\mathit{H}\right)=0$ in general holds if $\mathit{W}$ is strictly lower triangular, because $\mathit{H}$ is itself triangular with null diagonal.

As regards the convergence of LS estimator, multiplying Eq. (17) by $1/N$, under regularity conditions (stability and existence of second-order moments) one has that ${\mathrm{Plim}}_{N\to \infty }\left({\mathit{Z}}^{\prime }\mathit{Z}/N\right)={\Sigma }_z$, the covariance matrix of ${\mathit{z}}_{\left(s\right)}$. Whereas

$${\mathrm{Plim}}_{N\to \infty }\left({\mathit{Z}}^{\prime }\mathit{e}/N\right)=\mathrm{E}\left[{\mathit{z}}_{\left(s\right)}{e}_{\left(s\right)}\right]=\mathrm{tr}\left(\mathit{H}\right){\sigma }_e^2,$$

which is 0 if $\mathit{W}$ is triangular. Finally, under regularity conditions (${\Sigma }_{\mathit{z}}$ positive definite, $\left|\varphi \right|<1$ and ${\mathit{x}}_s$ stationary) one has

$$\sqrt{N}\left({\widehat{\mathit{\delta }}}_N-\mathit{\delta }\right)\underset{N\to \infty }{\overset{\mathrm{D}}{\to }}\mathrm{N}\left(\mathit{0},{\Sigma }_{\mathit{z}}^{-1}{\sigma }_e^2\right),$$

where the convergence is in distribution (D), see Yao and Brockwell (2006).

For multilateral SAR models ($\mathit{W}$ non-triangular), the LS is generally biased and inconsistent, and alternative estimators must be sought. The ML method may be suitable if the sample size N is moderate and the observations have Gaussian distribution. This estimator maximizes the log-likelihood function

$$\begin{array}{cc}\hfill \ell (\alpha ,\varphi ,\mathit{\beta },\sigma )& \propto -Nlog\left(\sigma \right)+log\left[det\left({\mathit{I}}_N-\varphi \mathit{W}\right)\right]-\left({\mathit{e}}^{\prime }\mathit{e}\right)/2{\sigma }_e^2,\hfill \\ \hfill \mathit{e}(\alpha ,\varphi ,\mathit{\beta })& =\left({\mathit{I}}_N-\varphi \mathit{W}\right)\mathit{y}-\left(\alpha {\mathit{1}}_N+\mathit{X}\mathit{\beta }\right),\hfill \end{array}$$

(19)

which requires $det(·)>0$, i.e. the stability condition $\varphi \in [0,1)$ (LeSage and Pace, 2009, p.63). This requirement is not necessary in the LS method; moreover, numerical issues may arise in the optimization of Eq. (19), when computing $logdet(·)={\sum }_{i=1}^{N}log\left(1-\varphi {\lambda }_i\right)$, with the eigenvalues ${\lambda }_i$ of $\mathit{W}$. This techniques cannot be applied to triangular matrices, because ${\lambda }_i$=0 for all i.

Another method which is suitable for parameter estimates in multilateral $\mathit{W}$ is that of generalized moments (GM, Kelejian and Prucha, 1999); it arises by applying the instrumental variables (IV) method to the LS estimator (17). The first step is the definition of valid IVs: on the basis of the predictor decomposition $\widehat{\mathit{y}}={\sum }_k{\varphi }^k{\mathit{W}}^k\left(\mathit{X}\mathit{\beta }+\alpha \right)$, they may be given by $\mathit{M}=[\mathit{X},\mathit{W}\mathit{X},{\mathit{W}}^2\mathit{X}\cdots ]$. This provides the projection matrix $\mathit{P}=\mathit{M}{\left({\mathit{M}}^{\prime }\mathit{M}\right)}^{-1}{\mathit{M}}^{\prime }$ and the regressors $\mathit{V}=\mathit{P}\mathit{Z}$, which lead to the IV estimator ${\tilde{\mathit{\delta }}}_N={\left({\mathit{V}}^{\prime }\mathit{Z}\right)}^{-1}{\mathit{V}}^{\prime }\mathit{y}$. This solution is also of GM-type, as minimizes the functional $Q_N\left(\mathit{\delta }\right)={\mathit{m}}^{\prime }{\mathit{S}}^{-1}\mathit{m}$, with the empirical moments $\mathit{m}={\mathit{M}}^{\prime }\mathit{e}/N$ and $\mathit{S}={\mathit{M}}^{\prime }\mathit{e}{\mathit{e}}^{\prime }\mathit{M}/N$. In the following, we check the effectiveness of ML and GM estimators with simulation experiments.

3.2. Forecasting Functions

In forecasting, we have L out-of-sample units with value ${\mathit{y}}_L=\mathit{0}$, whose coordinates and regressors ${\mathit{x}}_l^{\prime }=\left[1,i_l,j_l,x_{1l}...x_{ml}\right]$ are known for all $l=1,2\cdots L$. Their locations $(i_l,j_l)$ may be inside or outside the observed region; in both cases they are placed at the end of the data matrix ${\mathit{X}}_{N+L}$. If the data are ordered in a certain direction (e.g. north-south with $i_{\left(s\right)}\ge i_{(s+1)}$), and all L units are outside the observed region, and in the same direction, then the weight matrix is nearly block-diagonal: ${\mathit{W}}_{N+L}=\mathrm{diag}\left[{\mathit{W}}_N,{\mathit{W}}_L\right]$, where ${\mathit{W}}_L$ is overlapped to ${\mathit{W}}_N$ for the contiguity of L units with the in-sample ones. In the other cases, it has a more complex structure, with triangular sub-matrices under the unilateral constraint

$${\mathit{W}}_{N+L}=\left[\begin{array}{c}{\mathit{W}}_{NN}{\mathit{O}}_{NL}\\ {\mathit{W}}_{LN}{\mathit{W}}_{LL}\end{array}\right]$$

In any case, the estimation of parameters $\mathit{\delta }$ is performed with in-sample data ${\mathit{y}}_N$ and the matrix ${\mathit{W}}_{NN}$; the other blocks will be used in forecasting.

The forecasting function of ${\mathit{y}}_L$ depends on the SAR representation used for $\mathit{y}$; Kelejian and Prucha (2007) also considered models with autocorrelated errors $\mathit{u}$. In the reduced (MA) form (13c), the in-sample (fitted values) and out-of-sample forecasts are jointly computed as

$$\begin{array}{cc}\hfill {\tilde{\mathit{y}}}_L& =\mathrm{E}\left({\mathit{y}}_L|{\mathit{X}}_{N+L},{\mathit{W}}_{N+L}\right),\hfill \\ \hfill {\tilde{\mathit{y}}}_{N+L}& ={\left({\mathit{I}}_{N+L}-\varphi {\mathit{W}}_{N+L}\right)}^{-1}\left(\alpha {\mathit{1}}_{N+L}+{\mathit{X}}_{N+L}\mathit{\beta }\right),\hfill \\ \hfill & ={\left[{\tilde{\mathit{y}}}_N^{\prime },{\tilde{\mathit{y}}}_L^{\prime }\right]}^{\prime },\hfill \end{array}$$

(20)

where the contiguity matrix may have q-entries per row (the spatial lags) and may use non-uniform (IDW) weights. The matrix solution (20) is nearly automatic, but in absence of exogenous variables it provides constant forecasts.

The second predictor comes from the structural (AR) representation (13b); it is not automatic and must be managed sequentially by rows as Eq. (13a)

$$\begin{array}{cc}\hfill {\widehat{\mathit{y}}}_L& =\mathrm{E}\left({\mathit{y}}_L|{\mathit{y}}_N,{\mathit{X}}_{N+L},{\mathit{W}}_{N+L}\right),\hfill \\ \hfill {\widehat{y}}_{N+l}& =\alpha +\varphi {\mathit{w}}_{N+l}^{\prime }{\widehat{\mathit{y}}}_{N+l-1}+{\mathit{x}}_{N+l}^{\prime }\mathit{\beta },l=1,2\cdots L,\hfill \\ \hfill {\widehat{\mathit{y}}}_{N+l-1}& ={\left[y_1,y_2,\cdots ,y_N,{\widehat{y}}_{N+1}\cdots {\widehat{y}}_{N+l-1},0,\cdots ,0\right]}^{\prime },\hfill \end{array}$$

(21)

where ${\mathit{w}}_l^{\prime },{\mathit{x}}_l^{\prime }$ are the l-th rows of the matrices ${\mathit{W}}_{N+L},{\mathit{X}}_{N+L}$. Note that for non-triangular $\mathit{W}$, Eq. (21) involves missing values in the running vector ${\widehat{\mathit{y}}}_{N+l-1}$. If these values are provided by Eq. (20), then the vector of observations and forecasts at l-th step becomes ${\widehat{\mathit{y}}}_{N+l-1}={\left[y_1\cdots y_N,{\widehat{y}}_{N+1}\cdots {\widehat{y}}_{N+l-1},{\tilde{y}}_{N+l}\cdots {\tilde{y}}_{N+L}\right]}^{\prime }$ and, at the end $l=L$, it will only contain the improved forecasts (21).

LeSage and Pace (2004) and Goulard et al. (2017) also discussed solutions based on the best linear unbiased predictor (BLUP) of Goldberger. This approach arises from the conditional mean and variance of Gaussian random vectors:

$$\begin{array}{cc}\hfill \mathrm{E}\left({\mathit{y}}_L|{\mathit{y}}_N\right)& ={\mathit{\mu }}_L+{\Sigma }_{LN}{\Sigma }_{NN}^{-1}\left({\mathit{y}}_N-{\mathit{\mu }}_N\right),\hfill \\ \hfill {\widehat{\tilde{\mathit{y}}}}_L& ={\tilde{\mathit{y}}}_L+{\Sigma }_{LN}{\Sigma }_{NN}^{-1}\left({\mathit{y}}_N-{\tilde{\mathit{y}}}_N\right),\hfill \end{array}$$

(22)

where ${\tilde{\mathit{y}}}_L$ is a subvector the predictor in Eq. (20) and ${\Sigma }_{LN},{\Sigma }_{NN}$ come from the partition of the joint covariance matrix of ${\mathit{y}}_{N+L}$, given by

$${\Sigma }_{N+L}={\left({\mathit{I}}_{N+L}-\varphi {\mathit{W}}_{N+L}\right)}^{-1}{\left({\mathit{I}}_{N+L}-\varphi {\mathit{W}}_{N+L}^{\prime }\right)}^{-1}{\sigma }_e^2=\left[\begin{array}{cc}{\Sigma }_{NN}& {\Sigma }_{NL}\\ {\Sigma }_{LN}& {\Sigma }_{LL}\end{array}\right].$$

(23)

The solution (22)-(23) requires a double matrix inversion: first ${\left({\mathit{I}}_{N+L}-\varphi {\mathit{W}}_{N+L}\right)}^{-1}$ and next ${\Sigma }_{NN}^{-1}$; in the presence of large N, this may introduce approximation errors. To reduce numerical instability, one may partition the simpler array ${\Omega }_{N+L}={\Sigma }_{N+L}^{-1}$ and exploit the algebric relationship ${\Sigma }_{LN}{\Sigma }_{NN}^{-1}=-{\Omega }_{LL}^{-1}{\Omega }_{LN}$ (see LeSage and Pace, 2008). The latter involves a single matrix inversion ${\Omega }_{LL}^{-1}$, of a submatrix of lower size $L\ll N$. Finally, the conditional variance $\mathrm{V}\left({\mathit{y}}_L|{\mathit{y}}_N\right)$ of the BLUP estimator may provide the dispersion of the forecasts (22)

$$\mathrm{V}\left({\widehat{\tilde{\mathit{y}}}}_L\right)={\Sigma }_{LL}-{\Sigma }_{LN}{\Sigma }_{NN}^{-1}{\Sigma }_{NL},$$

this matrix requires the estimation of the variance ${\sigma }_e^2$, which can be carried out with the in-sample residuals ${\widehat{\mathit{e}}}_N=\left({\mathit{y}}_N-{\widehat{\mathit{y}}}_N\right)$ of Eqs. (20).

As a final remark, notice that all formulae displayed until now can be extended to unconstrained SAR(q) models, by replacing the IR matrix with one of order q-th. For the predictors (20),(21) this yields the following dispersion

$$\mathrm{V}\left(\tilde{\mathit{y}}\right)={\left({\mathit{I}}_{N+L}-\sum _{k=1}^q{\varphi }_k{\mathit{W}}_k\right)}^{-1}{\left[{\left({\mathit{I}}_{N+L}-\sum _{k=1}^q{\varphi }_k{\mathit{W}}_k\right)}^{-1}\right]}^{\prime }{\sigma }_e^2,$$

where the k-th lag contiguity matrices ${\mathit{W}}_k$ have size $N+L$.

3.3. Simulations and Applications

To compare the various estimators, we perform simulation experiments of SAR models (13)-(16) on a random grid of N=150 points, unformly distributed in the unit square: $[i_s,j_s]\sim \mathrm{IU}{(0,1)}^2$, as in Figure f6c. We generate realizations of the process $\left\{y_s\right\}$ with independent inputs $x_s\sim \mathrm{IU}(0,2)$ Uniform and $e_s\sim \mathrm{IN}(0,1)$ Normal; parameters $\alpha =1,\varphi =0.65,0.95,\beta =-1$; contiguity matrices $\mathit{W}$ with q=1,3 lags, with north NN and simple NN links. Subsequently, M=500 replications are fitted with LS, ML, GM estimators and relative biases, relative RMSE and p-values of Normality test were computed for the coefficents $\widehat{\alpha },\widehat{\varphi },\widehat{\beta }$ and then averaged. Note that relativization (e.g. bias$\left({\widehat{\varphi }}_N\right)/\left|{\varphi }_0\right|$) allows to sum the individual statistics so as to directly compare the various methods.

The results are reported in Table Section 2.2, where "variable grid" means that at each replication the centroids change and "variable weights" indicates that IDW are used in the matrices $\mathit{W}$. The Durbin model includes a lagged term on the exogenous part, as ${\beta }_1\mathit{x}+{\beta }_2\mathit{W}\mathit{x}$, with coefficients $[{\beta }_1,{\beta }_2]=[1,-.65]$. The software used for ML and GM estimatess is the Matlab package of LeSage and Pace (2009); it fits the Durbin model only with the ML method. Main conclusions from Table Section 2.2 are that ML and GM methods do not improve the LS estimates in the case of unidirectional (triangular) $\mathit{W}$ matrices. However, they significantly outperform LS in the multidirectional (e.g. simple NN) case; in particular, the ML method is the best. However, these results are not homogeneous as regards unbiasedness and efficiency of estimates, where the GM estimator may locally have smaller bias.

Forecasting. Unlike estimation, prediction evaluation of SAR methods cannot be based on pre-specified models and contiguity matrices, but it must refer to data generated from autonomous sources. As in the lattice case, we considered the random surfaces (11), sampled at random points $i_s\sim \mathrm{IU}[1,n],j_s\sim \mathrm{IU}[1,m]$, with $s=1,2\cdots N$; however, the results are too favourable to unilateral SAR models. We then turn to the real USGS (2007) data in Figure 4, with n=340, m=455; we sample N=150 cells and withhold for forecasting L=30 outside, see Figure 7a. The fitted model is SAR (15) with q=3 and ${\mathit{x}}_s={[i_s,j_s]}^{\prime }$, and forecast functions (20)-(22); the unilateral model is estimated with LS and the multilateral one with ML.

As a final application, we consider original point data concerning the measurement of stable isotopes of oxygen (${\delta }^{18}$O) and hydrogen (${\delta }^2$H) in the ground water. Mapping isoscapes is useful for physical monitoring of the hydrological cycle and for archeological and forensic investigations, regarding the path of people movement. We consider the datasets of Gautam et al. (2020), recorded in South-Korea in 2010, with N=130; and Moreiras Reynaga et al. (2021), recorded in Mexico in 2007, with N=234. We estimate SAR model (15) with q=3 for $y_s={\delta }^{18}\mathrm{O},{\delta }^2\mathrm{H}$ and ${\mathit{x}}_s$ = latitude, longitude, to forecast about 20% of observations; the results are in Tab. 6 and Fig. 8. They confirm the better performance of unilateral SAR models with the AR algorithm (21), especially when the data have significant autocorrelation (i.e. large $\left|\varphi \right|$) and marked spatial trends (significant ${\beta }_{1,2}$). The reduction of MAPE staristics in Tab. 6 ranges from -18% in Korea to -38% in Mexico.

Table 6. Parameter estimates (and T-statistics) and MAPE forecast statistics of the SAR model (15) with q=3 applied to the water isotope data of South Korea ${\delta }^2$H (first block) and Mexico ${\delta }^{18}$O (second block).

Model	Method	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\varphi }}$	${\widehat{\mathit{\beta }}}_1$	${\widehat{\mathit{\beta }}}_2$	${\widehat{\mathit{\sigma }}}_{\mathit{e}}$	$\begin{array}{c}\mathit{MAPE}\left(20\right)\end{array}$	$\begin{array}{c}\mathit{MAPE}\left(21\right)\end{array}$	$\begin{array}{c}\mathit{MAPE}\left(22\right)\end{array}$
Unilateral	LS (17)	123.7	0.337	-0.557	-2.51	4.28	0.069	0.065	0.065
T-stat., $R^2$	"	(1.70)	(4.72)	(-1.04)	(-3.37)	0.478	.	.	.
Multilat.	ML (19)	97.9	0.622	-0.376	-1.98	3.53	0.080	0.079	0.079
T-stat., $R^2$	"	(1.69)	(9.46)	(-0.86)	(-3.41)	0.645	.	.	.
Unilateral	LS (17)	4.30	0.636	0.089	0.093	1.28	0.377	0.369	0.370
T-stat., $R^2$	"	(1.69)	(8.17)	(3.19)	(2.72)	0.336	.	.	.
Multilat.	ML (19)	-0.85	0.479	0.053	0.101	1.26	0.600	0.597	0.588
T-stat., $R^2$	"	(-0.36)	(7.41)	(1.88)	(2.89)	0.351	.	.	.

Table 6. Parameter estimates (and T-statistics) and MAPE forecast statistics of the SAR model (15) with q=3 applied to the water isotope data of South Korea ${\delta }^2$H (first block) and Mexico ${\delta }^{18}$O (second block).

Model	Method	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\varphi }}$	${\widehat{\mathit{\beta }}}_1$	${\widehat{\mathit{\beta }}}_2$	${\widehat{\mathit{\sigma }}}_{\mathit{e}}$	$\begin{array}{c}\mathit{MAPE}\left(20\right)\end{array}$	$\begin{array}{c}\mathit{MAPE}\left(21\right)\end{array}$	$\begin{array}{c}\mathit{MAPE}\left(22\right)\end{array}$
Unilateral	LS (17)	123.7	0.337	-0.557	-2.51	4.28	0.069	0.065	0.065
T-stat., $R^2$	"	(1.70)	(4.72)	(-1.04)	(-3.37)	0.478	.	.	.
Multilat.	ML (19)	97.9	0.622	-0.376	-1.98	3.53	0.080	0.079	0.079
T-stat., $R^2$	"	(1.69)	(9.46)	(-0.86)	(-3.41)	0.645	.	.	.
Unilateral	LS (17)	4.30	0.636	0.089	0.093	1.28	0.377	0.369	0.370
T-stat., $R^2$	"	(1.69)	(8.17)	(3.19)	(2.72)	0.336	.	.	.
Multilat.	ML (19)	-0.85	0.479	0.053	0.101	1.26	0.600	0.597	0.588
T-stat., $R^2$	"	(-0.36)	(7.41)	(1.88)	(2.89)	0.351	.	.	.

Figure 8. Water isotopes data and SAR forecasts (22): Red=out-of-sample data, Green= forecasts with $\mathit{W}$ triangular; Black=forecasts with $\mathit{W}$ multilateral. Panels: a) South-Korea sample locations; b) Latitudinal view of Hydrogen isotope; c) Mexico sample locations; d) Longitudinal view of Oxigen isotope.

Figure 8. Water isotopes data and SAR forecasts (22): Red=out-of-sample data, Green= forecasts with $\mathit{W}$ triangular; Black=forecasts with $\mathit{W}$ multilateral. Panels: a) South-Korea sample locations; b) Latitudinal view of Hydrogen isotope; c) Mexico sample locations; d) Longitudinal view of Oxigen isotope.

4. Conclusions

In this paper we have compared unilateral and multilateral SAR models in forecasting spatial data of both lattice and polygonal type. SAR systems are natural extensions of classical autoregressions, their difference is in the treatment of the spatial dependence: while unilateral models choose a single direction in space, as in time series, multilateral models consider multiple directions. These approaches lead to different contiguity matrices: triangular and sparse, usually computed with the nearest-neighbor approach. While triangular SAR models can be consistently estimated with least squares, multilateral SAR require maximum likelihood or method of moments; in the latter cases, numerical complexity increase with the dimension of the contiguity matrices. Instead, LS is not sensitive to the size of the lattice or the number of spatial units involved and can manage even large scale systems, see Eqs. (9) and (17). Simulation experiments on small-medium scale systems have shown that ML and GM are suitable for multilateral models, but LS is preferable for the unilateral ones, especially in presence of unstable systems ($\left|\varphi \right|\ge 1$).

As in time series, unilateral models with data ordered in a certain direction, yield recursive (sequential) calculations, which allow simple and unconditional forecasting functions. Instead, multilateral models involve forward values which may be pre-estimated only with the reduced form (20). In the applications we have seen that out-of-sample forecasting of multilateral models benefits from the BLUP improvement (22), but the performance of the unilateral models with (21) remains better, see Tabs. 5 and 6. As regards the practical usage of SAR point forecasts, we mention the possibility to use them as alternatives to non-parametric smoothers, such as Kriging and Kernels, which are particularly sensitive to the border effects. Further, in order to overcome the limits of unidirectionality, one can estimate unilateral models in various directions (e.g. NS, SN, WE, EW), and then combine their forecasts. This solution is suitable for points inside the perimeter of the investigated area, but it can also be applied to external points by using diagonal paths (e.g. NE).

These are open topics for further research; other important areas of application of triangular contiguity matrices are to space-time (ST) data; in particular, to STAR and dynamic panel-data models (see Grillenzoni, 2004 and Copiello and Grillenzoni, 2021). In the latter reference it is shown how the unidirectional specification improves the performance even of robust (M-type) estimators, designed for data with outliers. Finally, the unilateral specification is also suitable for testing the Granger causality between spatial processes (e.g. Herrera et al., 2012) because it is related to diffusion effects which follow specific directions in the space.