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.

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 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 the LS estimator becomes

At computational level, the double summation in Eq. (2) can be implemeted in recursive form as follows 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 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 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 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.

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 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 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

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 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.

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 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$.

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.