A Simple Test of Spatial Autocorrelation for Centered Variables

We present a simple test of spatial autocorrelation based on the skedastic structure of the spatial series. Its distribution function is known for all sample sizes. Moreover, it is very simple to obtain, specially in a case of small samples where the new GQsp test has great power, higher than other alternatives existing in the literature. Article History: Received: May 13 2020 / Revised: May 25 2020 / Accepted: August 13 2020


Introduction
One of the purposes of the analysis of spatial series is the detection of cross-sectional relations among the observation of the series, as a previous step to build, for example, an econometric model. This is what we call spatial autocorrelation. It is also well-known that spatial observations lack of a natural order, which difficult the analysis. The solution in the field of spatial econometrics consists in the specification of the so-called weighting or contiguity matrix, W. The purpose of this matrix is to reflect the network of cross-sectional dependencies, which is, per se, a non-observable phenomenon. There is great flexibility in the definition of the weights (see, for example, Griffith, 1996;Anselin, 2002) which can be symmetric or not, binary or continuous, etc; the only restrictions are that they must be non-negative and the main diagonal of this matrix should be made of zeros (a region does not interact with itself). Usually, these weights are specified using the criterion of geographical proximity in the sense that, for example, a pair of regions physically contiguous, or separated by a certain distance, are supposed to be neighbours, able to interact. We can use 1 for these pairs of regions and 0 otherwise. The result will be a binary symmetric matrix that may also be row-standardized so that each row sums to one; and that is non symmetric in this case (see Harris et al., 2011, for a more extended discussion on W).
The mechanisms of spatial dependence can be of different types (Cliff and Ord, 1981) but the most popular are the autoregressive, SAR, and moving average, SMA, equations: where B sar = (I − δW) and B sma = (I + δW), y is the (n × 1) vector of observations of the series, u is a (n × 1) vector of error terms (that for simplicity we assume u ∼ N (0, σ 2 I)), and δ is a parameter of spatial autocorrelation.
The assumption of independence implies, in both cases, that δ = 0, which is the null hypothesis in most of the spatial autocorrelation tests that exist in the literature (Anselin and Florax, 1995;Kelejian and Piras, 2017). Especially outstanding among this is the Moran's I test (Moran, 1950), whose characteristics (simplicity, reliability) gave him a prominent role among practitioners. However, it is very difficult to obtain its distribution function in finite samples (Sen, 1990;Tiefelsdorf and Boots, 1995;Tiefelsdorf, 2000), which forces to the researchers to use the (asymptotic) normal approximation (Cliff and Ord, 1981;Kelejian and Prucha, 2001).
In this paper we present a new test of spatial autocorrelation, called GQ sp in reference to the traditional Golfeld-Quandt test of heteroskedasticity (Goldfeld and Quandt, 1965), whose distribution function is known for all sample sizes and which appear to have better properties than the Moran's I.

Heteroskedasticity and Spatial Autocorrelation
A peculiarity of the spatial framework to which this paper addresses is that series with a spatial structure will not only show features of cross-sectional dependencies, but also of heteroskedasticity which is evident in the covariance matrix of the series in (1) and (2). However, until now, the literature has mostly focused on the first aspect. In fact, Moran's I tests if the covariance between the series and its spatial lag is not statistically different from zero. The information about the skedastic variances is not used at all.
Matrix W is known because it has been supplied by the user. Let us assume that W is a binary symmetric matrix based, for example, in first order contiguity. This type of matrices can be decomposed in the usual way using the matrices of eigenvectors, Q, and eigenvalues, Λ which is diagonal, so that W = QΛQ ′ ; both matrices are (n × n). In a way similar to Griffith (1996Griffith ( , 2000, we may use the Q matrix to filter so thatỹ = Q ′ y. The result of this simple transformation is that the spatial cross-correlation will be totally removed but the skedastic nature will be appreciated more clearly. For example, in the case of the SAR series of (1) we obtain: beingũ ∼ N (0, σ 2 I) because of the property of orthogonality of the eigenvectors, Q ′ Q = I. In sum,ỹ ∼ N (0, σ 2 ∆ −2 ). It is clear that the filtered series is spatially independent. Moreover, the skedastic function that rules its variances is V (y r ) = σ 2 (1−δηr) 2 , being η r the r-th eigenvalue of W. A similar result is obtained for the SMA case of (2): whereỹ ∼ N (0, σ 2 ∆ 2 ). The difference is that, now, the skedastic function is proportional to the eigenvalues, V (y r ) = (1 − δη r ) 2 . What we propose is to use this information to develop a test of spatial autocorrelation capable of exploiting the skedastic structure of the spatially filtered series.
To progress in this direction we can use the well-known Golfeld-Quandt statistic (Goldfeld and Quandt, 1965) with good properties as a heteroskedasticity test but, in this case, employed as a test of spatial autocorrelation with the following null and alternative hypotheses: For this it is only necessary to: 1. Order the values of the filtered series,ỹ r , according to the values of the associated eigenvalues, η r , with r = 1, 2, . . . , n. The ordering may be ascending or descending, it does not matter.
2. Remove the m central observations to increase the power of the test, so that we have two sub-samples of equal size with n−m 2 observations in each one. The applied literature suggest removing one third of the sample: m ≈ n 3 .

Solve the LS estimation of the equation:
in each sub-sample.τ r is an artificial variable obtained asτ = Q ′ τ , being τ a (n × 1) vector of ones, whose purpose is to account for the existence of a constant term in the right-hand side of equations (1)-(2). Let us call the respective sum of squares of each regression as SS j , with j = 1, 2. If we are sure that the original series is centered around zero, we can avoid the LS estimation of (6) and obtain the sum of squares directly fromỹ.
4. Obtain the Golfeld-Quandt as usual: where k 1 and k 2 are the respective degrees of freedom, under the null hypothesis of (5), of both sum of squares. If the original variable was centered, it is immediate to obtain: Given that, under the null hypothesis, there is no spatial correlation in the sample and the two sub-samples formed in the second step are independent, the distribution function of the GQ sp statistic of (7) is F ( n−m 2 ; n−m 2 ) . If the original series was not centered, the sum of squares would have been obtained using the LS residual from the equation (6) and we would have to discount one degree of freedom from the F distribution which is F ( n−m 2 −1; n−m 2 −1) . Then, the rule of decision is: being F the abscissa of the F distribution, with the corresponding degrees of freedom, and with a probability mass of to its right.
Appendix A contains a more detailed discussion in relation to the sources of power for the GQ sp test. Next section continues with a Monte Carlo experiment.

A Monte Carlo Study
The GQ sp statistic has been introduced in the previous section as a test of spatial autocorrelation for a given spatial series, but using the filtered series. In fact, the matrix of eigenvectors applied to the original series removes the spatial autocorrelation structure strengthening the skedastic feature in the filtered series. That is the reason to adapt the traditional Golfeld-Quandt statistic to the problem of detecting spatial autocorrelation. We expect a good behaviour of the GQ sp statistic in cases of small sample sizes. To assess the properties of this new test, we compare its results against the well-known Moran's I, one of the most popular test in this field.
The main characteristics of the experiment are the following: a. Spatial autocorrelation processes: Two types of equations, SMA and SAR, have been used in the study: where τ is a (n×1) vector of ones and µ is a parameter that indicates if the series is centered on zero or not. e. Error term: In a spatial setting, it is interesting to check for the robustness, in our case, of these two test to the assumption of normality so we have simulated two cases of Normality, f. Variance of the error term: Three possible cases of interest have been simulated, according to: • Spatially structured heteroskedasticity: V [u s ] = σ 2 s ; σ s = ∑ s≠r b sr φ r ; φ s ∼ U (0, 1), being b sr a sequence of non-negative weights with b ss = 1, s = 1, 2, . . . n, and b rs = λω rs , i.e., that recreate a SMA process. g. Number of draws: The experiment consists of 1000 draws for each case.
In sum, the Monte Carlo implies two different processes (SMA and SAR equations), 2 different scale factors, 5 sample sizes, 2 distribution functions, 3 types of variances and 51 different values of the spatial autocorrelation coefficient; i.e., 6120 different configurations in total. Moreover, the W matrix corresponds to a first order contiguity matrix (consequently, symmetric) on a hexagonal lattice.

Results of the MC
To facilitate the discussion, we group the results in four cases: Case I which is the standard case of a SMA or SAR process with normal and homoskedastic disturbances; Case II corresponds to Log-normal errors whereas; Cases III and IV have heteroskedastic disturbances, with random and spatially structured variance. Moreover, we omit the results corresponding to large sample sizes (i.e., n = 400, 900) because the power is almost uniformly 100% without size problems. Same applies for the scale factor where we did not notice any remarkable difference depending on the value of µ. Thus, we concentrate on small/medium sample sizes (i.e., n = 16, 25, 100), with a unitary factor of scale, µ = 1. Table 1 shows the estimated sizes for the 4 cases and 2 processes. Overall, these results confirm the tendency of the Moran's I to underestimate the size of the test for samples of small size. Case I corresponds to the standard case and the estimated size; for sample sizes of n = 16 and n = 25 is well below the statistical limits for a significance level of 5%. Only for a sample of medium size, n = 100, this estimate produces acceptable values. Log-normality has a strong impact on these results, reducing significantly these estimates. The situation repeats for the case of heteroskedastic disturbances, although a random structure in the variance seem to have a stronger impact than a spatially structured variance. The type of spatial process, whether SMA or SAR, does not make any remarkable difference. In sum, only in 5 cases, up to a total of 24, the estimated size pertains to the 5% significance interval for a theoretical value of 5%, which amounts to a poor 21%. The average estimated size for the 24 cases is just 2.1%.
The percentage of correct estimates rises to 67% in the case of the GQ sp test whereas the average estimated size is 4.3%, inside the significance interval. All the estimates are correct for Cases I and IV. Once again, the log-normal errors is the worst situation for this test, with a strong tendency to underestimate the theoretical size. Moreover, there are hardly differences depending on the spatial process, and a random variance has worst consequences than a spatially structured variance. In sum, it is clear that GQ sp has less problems to correctly estimate the size, even in cases of very small sample sizes (n = 16).      1. Most of the power functions are strongly asymmetrical, specially for small sample sizes (n = 16) and SMA processes. The GQ sp test has a good behaviour for the range of negative values in the spatial autocorrelation coefficient, sometimes better than in the range of positive values. This is not the case with the Moran's I, whose performance for the case of negative spatial autocorrelation is pretty disappointing.
2. Sample size has a beneficial impact in the behaviour of the two tests, so that both become highly credible for medium sample sizes, n = 100, no matter the type of spatial process. Results are also interesting for a sample size of n = 25, but only for SAR processes.
3. The estimated power function of the GQ sp test is usually above that of the Moran's I. The major differences, greater than 50 points, occur for very small sample sizes, SMA processes and negative spatial autocorrelation. Those differences reduce as the sample size increases and also for SAR processes and positive spatial autocorrelation. The two power functions become practically indistinguishable for a medium sample size, n = 100, no matter the type of spatial process or the sign of the spatial autocorrelation.
4. The introduction of anomalies in the data generation process has clear impact in the power functions of the two tests. These consequences are very strong for the case of small sample sizes, n = 16, but it also seriously reduces the credibility of the two tests for medium-sized samples where the power functions are more open.
5. The case of log-normality shown in Figure 2 is, probably, the worst. The power of both tests is negligible for SMA processes and small sample sizes, n = 16, although its impact is smaller for SAR processes. Sample size helps to improve these deficiencies. Moreover, we should remind the acute size problems caused by log-normality, which makes the behavior of tests even more irregular.
6. Something similar occurs due to heteroskedasticity, with strong anomalies at the extremes of the stability interval especially for the case of small sample cases, n = 16 and SMA processes. The negative impact is mitigated with the sample size, although the power functions of both tests are more open, which denotes a worsening of the estimated power. Note also that, in the two cases of random and spatially structured variance, the GQ sp had less problems with the estimated size than Moran's I, which improves its reliability.
7. In sum, the GQ sp is more robust to departures in the data generation process, either in terms of heteroskedaticity or log-normality.

Conclusions and Main Findings
This paper presents a new test of spatial autocorrelation, especially adequate for samples of small size (by samples of small size we mean 25 observations or even less). We call it GQ sp test, because it is a reformulation of the classical Golfeld-Quandt test of heteroskedasticity (Goldfeld and Quandt, 1965).
The key point to build the test is that a series with spatial autocorrelation, according to SMA or SAR processes, comes also with a very well defined heteroskedastic structure in the variance. The process is not immediate because, first, the series must be filtered by using the matrix of eigenvectors associated to the W matrix; then, the filtered series is free from spatial autocorrelation but maintains the skedastic nature depending on the roots of W. In these conditions, is relatively simple to adapt the Golfeld-Quandt test.
The GQ sp test needs to obtain previously the eigenvector matrix associated to W. However, the cost of this additional calculus alleviates because of the framework for which we propose the new test: the case of small sample sizes. The Monte Carlo solved in the paper confirms the adequacy of the Golfeld-Quandt to the new situation where it is clearly superior to the popular Moran's I, both in what respects to size and power.
Appendix A. Sources of Power for the GQ sp Test The filtered series used in the GQ sp test, assuming normality and a zero constant in the right-hand side of the equation can be written as: where σ 2 s = σ 2 (1 + δη r ) 2 in the SMA case and σ 2 s = σ 2 (1 − δη r ) −2 for SAR processes. The filtered series,ỹ has been ordered according to the eigenvalues (increasing or decreasing; it does not matter) of W. We have excluded the m central observations, obtaining two sub-samples, each one with K = n−m 2 observations. Moreover, we have obtained the sum of squares in each sub-sample, SR j = ∑ s∈jỹ 2 s ; j = 1, 2, and finally build the GQ sp statistic, GQ sp = SR 1 SR 2 , whose distribution under the null hypothesis that H 0 ∶ δ = 0 is a central F (K,K) .
Under the alternative hypothesis H A ∶ δ ≠ 0, the statistic GQ sp continues to be a centered F distribution but the distribution function are not (K, K), which confers power to the statistic. Let us assume a SMA process like that in (1), then is immediate to obtain: The negative eigenvalues will always appear in the first sub-sample, contributing to SR 1 , whereas the positive η s will appear in the second sub-sample (recall that the sum of the eigenvalues of W is zero, so there will be both positive and negative). For the case of δ < 0, the terms (1 + δη s ) in the first sub-sample will be predominantly greater than one, (1 + δη s ) > 1, whereas the opposite will occur in the second sub-sample where this terms will be, predominantly, smaller than one, 0 < (1 + δη s ) < 1. Note also that this sequence of terms are non-negative, (1 + δη s ) ≥ 0, and that they sum to n, ∑ s=n s=1 (1 + δη s ) = n. So, for δ < 0 it will happen that: In this case, we expect quotients for the GQ sp statistic substantially greater than one, far away from the expected value of the statistic under the null hypothesis of no correlation which is one. Table A.1 presents the four cases depending on the type of process and the sign of the spatial autocorrelation coefficient.
Given that E[ỹ s ] = 0, the F distribution for the GQ sp statistic will have a non-centrality parameter of zero. If, however E[ỹ s ] ≠ 0, the distribution of the statistic of (7) will present a non-centrality parameter which can be avoided just by centering the variable,ỹ s −ỹ, in which case the degrees of freedom of the F distribution are (K − 1; K − 1) or by using the LS residuals of a previous regression on k exogenous variables, in which case the degrees of freedom of the F distribution are (K − k; K − k) as it is traditional in the Golfeld-Quandt literature.
In any case, the sequence of observations included in each sum of squares are independent both under the null and under the alternative hypotheses, which allows us to obtain the power function of the GQ sp statistic by solving the probability: being q α the abscissa of the F distribution, with the corresponding degrees of freedom, with a probability mass of α to its right. In the case of a SMA series, the second probability can be written as: where k * i = σ 2 ∑ s∈j (1 + δη s ) 2 ; i = 1, 2 are the degrees of freedom of the respective sum of squares under the alternative hypothesis and q * = k * 2 k * 1 q . Similarly, the first probability of (A.4) leads us to: P r [q 1− ≥ GQ sp ] = P r q * 1− ≥ F (k * 1 ,k * 2 ) δ ≠ 0 . (A.6) All the elements in (A.5) and (A.6) are known by the user (they depend on the weighting matrix and the F distribution) so the probabilities can be easily computed. Note that the behaviour of the GQ sp depends crucially on the weighting matrix W selected by the user, through its eigenvalues, η s ; s = 1, 2, . . . n.
Finally, our experience with the GQ sp test tends to confirm the suggestion of Harvey and Phillips (1974) of removing one third of the central observations in order to guarantee an adequate equilibrium between power and reliability for the Golfeld-Quandt test.