Next Article in Journal
A Bibliometric Review of the Ordered Weighted Averaging Operator
Next Article in Special Issue
Advancing Green TFP Calculation: A Novel Spatiotemporal Econometric Solow Residual Method and Its Application to China’s Urban Industrial Sectors
Previous Article in Journal
State of Health Prediction of Electric Vehicles’ Retired Batteries Based on First-Life Historical Degradation Data Using Predictive Time-Series Algorithms
Previous Article in Special Issue
Effect of Collaborative Innovation on High-Quality Economic Development in Beijing–Tianjin–Hebei Urban Agglomeration—An Empirical Analysis Based on the Spatial Durbin Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

A Survey of Spatial Unit Roots

1
Center for Policy Research and Department of Economics, Syracuse University, 426 Eggers Hall, Syracuse, NY 13244-1020, USA
2
Department of Economics, Syracuse University, Syracuse, NY 13244-1020, USA
*
Author to whom correspondence should be addressed.
Mathematics 2024, 12(7), 1052; https://doi.org/10.3390/math12071052
Submission received: 26 January 2024 / Revised: 26 March 2024 / Accepted: 29 March 2024 / Published: 31 March 2024
(This article belongs to the Special Issue Mathematical Economics and Spatial Econometrics)

Abstract

:
This paper conducts a brief survey of spatial unit roots within the context of spatial econometrics. We summarize important concepts and assumptions in this area and study the parameter space of the spatial autoregressive coefficient, which leads to the idea of spatial unit roots. Like the case in time series, the spatial unit roots lead to spurious regression because the system cannot achieve equilibrium. This phenomenon undermines the power of the usual Ordinary Least Squares (OLS) method, so various estimation methods such as Quasi-maximum Likelihood Estimate (QMLE), Two Stage Least Squares (2SLS), and Generalized Spatial Two Stage Least Squares (GS2SLS) are explored. This paper considers the assumptions needed to guarantee the identification and asymptotic properties of these methods. Because of the potential damage of spatial unit roots, we study some test procedures to detect them. Lastly, we offer insights into how to relax the compactness assumption to avoid spatial unit roots, as well as the relationship between spatial unit roots and other models, such as the Spatial Dynamic Panel Data (SDPD) model and Lévy–Brownian motion.

1. Introduction

There is an extensive literature using spatial statistics that deals with cross-sectional correlation, and is popular in regional science, urban economics and geography, to mention a few; see Anselin [1] for a nice introduction to this literature. Unlike time-series, there is typically no unique natural ordering for cross-sectional data. Spatial dependence models may use a metric of economic distance that provides cross-sectional data with a structure similar to that provided by the time index in time series. Examples in economics usually involve spillover effects or externalities due to geographical proximity. For example, the productivity of public capital, like roads and highways, on the output of neighboring states. Also, the pricing of welfare in one state that pushes recipients to other states. In a linear regression model, this spatial correlation may be in the disturbances and is called the spatial error model (SEM), or modeled on the dependent variable itself and named the spatial autoregression (SAR) model, or the spatial lag model. Unlike the autoregressive model or lagged model in time series where there is a natural ordering across time and lagged values are well defined, in cross-sections, this is dealt with using neighbors whose shocks or disturbances are affected by their neighbor’s shocks or disturbances in the SEM. For the SAR, the dependent variable is affected by neighbors, like house price being affected by neighboring house prices.This leads to the construction of a weight matrix that defines one’s neighbors by distance or contagion; see Anselin [1]. So, rather than lagged house prices as in an autoregressive time series model, own house price is related to a weighted average of neighboring house prices. Spatial dependence has been extended from cross-section to panel data; see Chapter 13 of Baltagi [2] or Elhorst [3] for a textbook treatment of the subject.
OLS yields inconsistent estimates in the SAR model due to the endogeneity of the spatial lagged dependent variable and the disturbances. For the SEM model, OLS yields an unbiased but inefficient estimator. Because of the limitation of OLS, (Q)MLE is often used to estimate the spatial models as in Ord [4], Anselin [1] and Lee [5]. However, MLE is sometimes computationally intensive especially for large sample sizes because of the requirement to compute the Jacobian term in the likelihood function. Ord [4] proposes a simplified computational procedure only requiring the eigenvalues of the spatial weight matrix, but computing accurate eigenvalues is increasingly difficult for large n. Kelejian and Prucha [6] suggest a Generalized Method of Moments (GMM) estimator for the SEM with SAR structure. Alternatively, constructing instrumental variables (IVs) from the exogenous variables, Kelejian and Prucha [7] propose the GS2SLS estimator and Lee [8] discuss the best GS2SLS estimator by using more efficient IVs. Lee [9] considers a GMM estimator for the SAR models with exogenoues variables and showed it is more efficient than the 2SLS estimator and is as efficient as the ML estimator asymptotically. Spatial panel data models (with dynamic terms) are also estimated using MLE, 2SLS or GMM as in Yu et al. [10], Baltagi and Liu [11], Kapoor et al. [12].
However, for all of these estimation methods considered above, they constrain the parameter space of the spatial coefficient to limit the degree of correlation between units. This is because when such spatial correlation is too strong, the spatial echoes passing through each unit do not die out and so the system cannot achieve an equilibrium. When the spatial weight matrix is row-normalized by convention, such a constraint requires the spatial coefficient to be smaller than 1 in absolute value. In this survey, we summarize the developments that relax this constraint and allow the spatial coefficient to equal or sufficiently approach unity, which is known as the spatial unit root. This is of practical relevance because there are many cases where the spatial coefficients are close to 1. For example, Keller and Shiue [13] detect the inter-regional trade of Chinese rice and find that rice prices for different provinces are highly related with spatial coefficients lying between 0.9 to 0.95.
When (near) spatial unit roots exist, the standard estimation procedures are not necessarily reliable and statistical inference is invalid. To remedy such a case, Fingleton [14] avoids the circularity of the spatial weight matrix and conducts a Monte Carlo simulation to explore the performance of OLS estimation. Alternatively, Lee and Yu [15] artificially let the spatial autoregression coefficient sufficiently approach the unit roots and derive the asymptotic behavior of QMLE and 2SLS estimators. The spatial unit roots have also been generalized to Spatial Dynamic Panel Data (SDPD) model by Yu and Lee [16].
In order to investigate the possible spatial unit roots, several test procedures have been proposed. Fingleton [14] suggests a “very high” value of the Moran’s I statistic could be useful for testing spatial unit roots. Lauridsen and Kosfeld [17,18] propose a two-stage Lagrange Multiplier (LM) tests that distinguish the spatial unit roots from the stationary positive spatial correlation. By the fact that spatial impulses do not die out under spatial unit roots and lead to explosive variance, Beenstock et al. [19] numerically calculate the critical value even under an irregular spatial weight matrix. These tests have been extensively used in the literature; see Yesilyurt and Elhorst [20], Olejnik [21], Machado et al. [22], Beenstock and Felsenstein [23].
This paper is organized as follows. Section 2 introduces some basic concepts in spatial econometrics. Next, the parameter space of the spatial autoregressive coefficient and the corresponding singular points are considered. Stationarity and spatial cointegration concepts derived from spatial unit roots are introduced. General assumptions and the corresponding implications in spatial econometrics are discussed. Section 3 introduces the potential problems with the existence of spatial unit roots: spurious and nonsense regression. Section 4 investigates estimation methods and inference under spatial unit roots. Section 5 discusses how to test for spatial unit roots. Section 6 discusses spatial unit roots in the SAR model while Section 7 concludes.

2. Basic Concepts in Spatial Econometrics

Spatial models study the spatial dependence between units. In practice, the spatial weight matrix is used to describe such dependence. Let W n * be the n × n contiguity-based spatial weight matrix, i.e., w i j * = 1 if unit i and j are contiguous and 0 otherwise. Also, the diagonal elements are set to 0 by convention. In practice, the spatial weight matrix is generally row-normalized by w n , i j = w n , i j * j w n , i j * , so the row sum of W n will be one.
Different spatial model specifications have different implications. The SEM with spatial autoregressive (SAR) structure on the n × 1 error vector u n can be expressed as u n = λ 1 W n u n + ϵ n = I n λ 1 W n 1 ϵ n , where λ 1 is known as the spatial coefficient and satisfies some assumptions that will be introduced later, and ϵ n is the n × 1 independent and identically distributed (i.i.d.) innovations with variance σ ϵ 2 . The error covariance matrix for the u n with SAR structure is
Ω S A R = E u n u n = σ ϵ 2 I n λ 1 W n 1 I n λ 1 W n 1 = σ ϵ 2 B n B n 1 ,
where u n = ϵ n ( I n λ 1 W n ) 1 and B n = I n λ 1 W n . Though W n may be sparse, B n B n 1 is not necessarily so, thus the spatial covariance structure induced by such SEM model is classified as global. Conversely, a spatial moving average (SMA) specification for the error vector u n can be expressed as u n = λ 2 W n ϵ n + ϵ n = I n + λ 2 W n ϵ n , and the corresponding covariance matrix will be
Ω S M A = E u n u n = σ ϵ 2 I n + λ 2 W n + W n + λ 2 2 W n W n ,
including only W n and W n W n which are first and second order neighbors if W n is defined as first-order contiguity. Hence, such a model is generally classified as local. See Baltagi et al. [24].
Kelejian and Prucha [7] consider a “cross-sectional (first-order) autoregressive spatial model with (first-order) autoregressive disturbances” (SARAR) and is labeled as spatial autoregressive combined (SAC) model. If the right-hand side includes both the independent variable and the spatially lagged dependent variable then it is termed the mixed regressive, spatial autoregression (MRSAR) or mixed SAR model. The spatial Durbin model (SDM) includes both the spatial lagged dependent variable and independent variables. A full model labeled as general nesting spatial (GNS) model given in Elhorst [3] is
Y n = λ 1 W n Y n + α ι n + X n β + W n X n θ + u n , u n = λ 2 W n u n + ϵ n ,
where λ 1 is referred as spatial autoregressive coefficient (SAC) and λ 2 is called the spatial autocorrelation coefficient.

2.1. Parameter Space of the Spatial Autoregressive Coefficient

Consider the pure SAR model with data generating process (DGP):
Y n = λ W n Y n + ϵ n ,
where Y n is an n × 1 vector of observations on the dependent variable, W n is an n × n spatial weight matrix and ϵ n is an n × 1 vector of disturbances which are assumed to be i.i.d. ( 0 , σ ϵ 2 ) . The reduced form equation of Y n can be written as
Y n = S n 1 ϵ n ,
where S n = I n λ W n . So to guarantee the system achieves equilibrium, a crucial assumption is that the absolute value of the spatial coefficient λ is strictly less than 1 (see Kelejian and Prucha [6] (Assumption 2), Kelejian and Prucha [7] (Assumption 2)) to ensure the nonsingularity of S n . This assumption follows from a sufficient condition for the invertibility matrix in Horn and Johnson [25] (Corollary 5.6.16, p. 351):
Theorem 1. 
An n × n matrix A n is nonsingular if there exists a matrix norm · such that I n A n < 1 . If this condition is satisified, A n 1 = k = 0 ( I n A n ) k .
Thus, S n is invertible if there exists a matrix norm such that I n ( I n λ W n ) = λ W n < 1 . It is also well known that any norm of a matrix is larger than all of its eigenvalues. Let X n be the eigenvector matrix and ρ i , i = 1 , , n , be the egienvalues of W n , then
| ρ i | X n = ρ i X n = W n X n W n X n .
So it is easy to see that | λ | < 1 W n and therefore
1 ρ m i n < λ < 1 ρ m a x ,
because i n ρ i = tr ( W n ) = 0 so that ρ m i n < 0 < ρ m a x . A useful result is given in Ord [4]:
Theorem 2. 
If W n has eigenvalues ρ 1 , , ρ n , | ρ I n W n | = i = 1 n ρ ρ i . Then for S n = I n λ W n , det ( S n ) = | I n λ W n | = i = 1 n 1 λ ρ i .
Moreover, the log-likelihood function for Y n , given Y n = y in (4) is
( λ , σ 2 ) = n 2 ln 2 π σ 2 1 2 σ 2 y S n S n y + ln | S n | ,
and | S n | = i = 1 n ( 1 λ ρ i ) > 0 , of which a sufficient condition is λ ρ i < 1 , for all i. Again, since ρ m i n < 0 < ρ m a x , we obtain the range of λ that is given in (7).
However, either by Theorem 1 or 2, (7) is a sufficient condition for the invertibility of S n . The singular points of S n are 1 ρ 1 , , 1 ρ n by Theorem 2, and the number of these singular points are at most countably many as n . This raises a problem as stated in Kelejian and Robinson [26]. These singular points can be determined generally by the nth polynomial numerically, and to avoid inconsistency, they should be removed from the possible values of λ . Griffith [27] (p. 19) states that this condition also ensures stationarity, but Kelejian and Robinson [26] give a counter example showing that it does not when W n is a row-normalized double queen weight matrix.
When the matrix is row-normalized, Kelejian and Prucha [6] (Note 8, p. 120) show that ρ m a x = 1 by Geršgorin’s theorem, and typically, | ρ m i n | < 1 [26]. We assume | λ | < 1 , and this is why λ is generally interpreted as the spatial autocorrelation coefficient similar to its counterpart in time series.

2.2. Stationarity, Order of Integration and Cointegration

Stationarity is a key assumption in time series. Similarly, in spatial econometrics, when the stationary assumption does not hold, spurious or nonsense regression appears, as will be shown in the next section. Stationarity is tightly connected with S n . The formal definition of stationarity is given in Anselin [1] (p. 42):
Definition 1. 
A process is strictly stationary if any finite subset x i , x j , , x n from the stochastic process x i , i I has the same joint distribution as the subset x i + s , x j + s , , x n + s for any s, where s represents an uniform shift in time, space or time–space.
But we generally consider a weaker version, covariance stationarity. For the intuition of stationarity and the connection with the inverse of S n , see Beenstock et al. [19]. Consider the pure SAR model in (4) and (5), the covariance matrix for Y n is
Var ( Y n ) = E Y n Y n = σ ϵ 2 S n 1 S n 1 σ ϵ 2 B n ,
where B n = S n 1 S n 1 and S n 1 is defined in (5). Letting b k j be the ( k , j ) element of the matrix B n , by normalizing σ ϵ 2 = 1 , Y k has variance b k k and covariance b k j with Y j . By Definition 1, stationarity requires that b k k and b k j remain unchanged asymptotically (this implicitly assumes that both location j and k are far away from the edge). Note that by Theorem 1, we have
S n 1 = k = 0 ( I n S n ) k = k = 0 ( λ W n ) k = I n + λ W n + λ 2 W n 2 + .
So the stationarity assumption is equivalent to λ m k W n m k 0 , for all k as m . Since m < n and m , m represents the “remote” area, and W n m k is the m k step neighborhood of unit k. Thus, intuitively, stationarity requires that the shocks from far away locations will asymptotically not affect the epicenter area.
Another two concepts tightly connected with unit roots and stationarity are the order of integration and cointegration. The order of integration is originally a concept in time series that describes the minimum number of differences that a non-stationary process needs to be (covariance) stationary. Cointegration, on the other hand, describes the minimum order of integration of a combination of two or more series with the same order of integration. The formal definitions of the order of integration and cointegration are given in Hamilton [28]:
Definition 2. 
A time series is integrated of order d, denoted I ( d ) , if ( 1 L ) d X t is a stationary process, where L is the lag operator and 1 L is the first difference.
Definition 3. 
Time series X and Y are cointegrated of order I ( d , b ) , if both of them are I ( d ) , and there exists a cointegrated vector ( a , b ) such that a X + b Y I ( d b ) .
In time series, the lag operator L is defined by L X t = X t 1 because of the natural order of temporal data. In spatial econometrics, we regard W n as the spatial lag operator and I n W n as the first order spatial difference, see Anselin [1] (pp. 22–26). We also use I ( d ) and I ( a , b ) to refer to spatial integration and spatial cointegration respectively.
More specifically, for the pure SAR model in (4) with a row-normalized weight matrix, if λ = 1 , Y n I ( 1 ) since ϵ n is stationary. Also, suppose both Y n and X n are I ( 1 ) , but they have a long-term equilibrium relationship Y n = X n β + ϵ n , then obviously, ( X n , Y n ) are I ( 1 , 1 ) with cointegrated vector ( β , 1 ) .

2.3. Some Fundamental Assumptions

Different assumptions are made in spatial econometrics for different estimation methods. The most common ones are listed here, and the implications are explained.
Assumption 1. 
The disturbances ϵ 1 , , ϵ n are i.i.d. for all n (so uniformly) with zero mean and finite variance σ 2 . Additionally, fourth moments exist.
Assumption 2. 
The elements of the exogenous variables X n are uniformly bounded for all n. The lim n X n X n n exists and is nonsingular.
Assumption 3. 
The matrix S n is nonsingular.
The existence of up to the fourth moment of disturbances is needed to apply the central limit theorem for (a system) of the linear-quadratic form (see Kelejian and Prucha [29] (p. 226) and Kelejian and Prucha [30] (p. 63)). The nonsigularity of S n makes sure the system achieves an equilibrium as well as ensuring that the mean and variance of Y n exist.
Assumption 4. 
The matrices W n and S n 1 = ( I n λ W n ) 1 are uniformly bounded (UB) in both row and column sums for all n. (We say a matrix is UB in row (column) sums if its maximum row (column) sum is finite. This property preserves under finite matrix multiplication.)
The UB condition for W n implicitly assumes a limited number of neighbors for all units even as n , so the weight matrix W n is sparse for large n. This assumption is relaxed in Lee et al. [31] by introducing dominant (popular) units. In practice, the spatial units have a limited number of neighbors. Though sometimes w n , i j may be defined as the inverse of the distance between i, j physically or economically, w n , i j tends to be 0 between far away units as n increases. So in general this assumption is satisfied. The UB of S n 1 is to ensure the covariance matrix Var ( Y n ) in (9) is still UB, which limits the correlation between two different units since the UB property preserves under matrix multiplication.
Other assumptions to ensure identification conditions or the derivation of asymptotic distributions of estimators will be mentioned when needed.

3. Spurious Regression When (Near) Unit Roots Exist

The variance of Y n explodes when unit roots exist, and OLS estimation may perform unsatisfactorily: the estimators are inconsistent, the test statistics do not have familiar distributions, and may even converge to a constant. These phenomena have been studied extensively in time series, and similar symptoms occur in spatial econometrics.

3.1. Spurious Regression of Driftless Series and Spatial Integration

Fingleton [14] studies unit roots and spatial cointegration in spatial econometrics. Using Monte Carlo simulations, he finds that spatial unit roots will lead to a spurious regression and proves that when two vectors are spatially cointegrated, even running a regression on the error-correction model yields inconsistent estimates. Beenstock et al. [19] distinguish between the terms spurious regression and nonsense regression and argue that, Fingleton [14] refers to nonsense regression instead of spurious regression. When Y n and X n are driftless random walks, the nonsense regression occurs because of the increased variances of Y n and X n over time. On the other hand, the spurious regression occurs when Y n and X n are independent random walks with drift, which causes their means to increase over time. See also Mur and Trívez [32].
To run the simulation, two independent pure SAR processes Y n and X n containing spatial unit roots are generated separately as in (5). But as discussed in Section 2.1, S n 1 does not exist under a row-normalized weight matrix when λ = 1 . To avoid this, Fingleton introduces the “unconnected central cell”, which manually sets one row of the spatial weight matrix equal to 0 to avoid circularity. This is a time-series analogy because there is always a starting point in temporal data ( t = 1 ). By doing so, the singular point is slightly larger than 1, and the existence of S n 1 is ensured [32]. Regressing X n on Y n , the t-statistic and coefficient of determination R 2 show the significance of the parameter between two unrelated variables when spatial unit roots exist. Letting e be the OLS residuals, Moran’s I, defined as
I M o r a n = n i = 1 n j = 1 n w n , i j e W n e e e = e W n e e e ( when W n is row - normalized ) ,
is the spatial version of the Durbin–Watson statistic, and thus is a measure of spatial autocorrelation. The simulation results of Moran’s I show a high level of positive spatial autocorrelation in the residuals and evidence for the presence of a spurious regression.
To remedy this situation, spatial differencing is introduced to the SAR process with unit roots:
Δ Y n = γ Δ X n + ϵ n ,
where Δ Y n = Y n W n Y n and Δ X n = X n W n X n . When both Y n and X n are spatial I ( 1 ) processes, we have Δ Y n = ϵ Y and Δ X n = ϵ X . The regression of the first-order spatial difference variable is equivalent to regressing two independent I ( 0 ) processes, which should theoretically yield γ ^ = 0 .
Next, spatially cointegrated series are considered. To generate ( X n , Y n ) I ( 1 , 1 ) , the “error-correction representation” is used. The idea is adopted from Robert [33] that the existence of error-correction representation is a necessary and sufficient condition for a cointegrated time series. The spatial analogy is
Y n = W n Y n + c ( W n X n W n Y n ) + e 1 n , X n = W n X n + d ( W n X n W n Y n ) + e 2 n , e 1 n N 0 , σ 1 2 I n , e 2 n N 0 , σ 2 2 I n ,
where W n X n W n Y n is the equilibrium error assumed (for simplicity) stationary and hence the name “error-correction”. The spatial unit root series X n and Y n have a long-term equilibrium X n = Y n . Note that (13) has two equations and two unknowns and W n is a noncircular matrix.
Moran’s I statistic may act as a useful indicator for cointegration because the cointegration regression (regress Y n on X n ) involves endogenous variables. Also, the first-order regression is inappropriate because of omitted variable bias concerning the equilibrium error W n X n W n Y n . Rearranging (13) yields the appropriate specifications:
Y n W n Y n = c ( W n X n W n Y n ) + e 1 n , X n W n X n = d ( W n X n W n Y n ) + e 2 n .
But OLS estimation for either c or d is inconsistent because of the presence of a spatially lagged dependent variable, which is different from the traditional time series counterpart.

3.2. Spurious Regression with Deterministic Trends

Fingleton [14] studies the effect of spatial unit roots by simulation while Mur and Trívez [32] show that the variance of the spatial unit roots series explodes. For the DGP in (4)
Y n = λ W n Y n + ϵ n Y n = ( I n λ W n ) 1 ϵ n = S n 1 ϵ n ,
since the contiguity-based spatial weight matrix is symmetric, W n can be decomposed as W n = Q n Λ n Q n 1 no matter whether it is row-normalized or not. Thus S n = I n λ W n = Q n ( I n λ Λ n ) Q n 1 = Q n Δ n Q n 1 where Λ n is the eigenvalue matrix of W n and thus Δ n is also diagonal. So, Mur and Trívez [32] derive the variance of Y n as
Var ( Y n ) = σ 2 B n B n 1 = σ 2 Q n Δ n 1 Q n 1 Q n 1 Δ n 1 Q n = σ 2 Q n Δ n 1 Q n Q n 1 Δ n 1 Q n ,
with
Q n Δ n 1 = q 11 q 1 n q n 1 q n n diag 1 1 λ ρ 1 , , 1 1 λ ρ n = q 11 1 λ ρ 1 q 1 n 1 λ ρ n q n 1 1 λ ρ 1 q n n 1 λ ρ n .
Let m i j be the element of row i and column j of the matrix Q n Q n 1 , and the variance of the r-th element of Y n is then
Var Y n , r = σ 2 i n j n q i r q j r 1 λ ρ i 1 λ ρ j m i j .
If W n is row-normalized, at least one ρ i = 1 , which means that if the spatial unit roots exist and λ = 1 , Var ( Y n ) explodes. If W n is not row-normalized so that Q n is symmetric and orthogonal with Q n Q n = I n , Var ( Y n ) = σ 2 Q n Δ n 2 Q n , then Mur and Trívez [32] show the variance of the observation at r reduces to
Var Y n , r = σ 2 i = 1 n q i r 2 1 λ ρ i 2 = , when λ = 1 ρ i .
But when λ is not one of the singular points of S n , Var Y n is not necessarily a function of n, i.e., the variances of Y n do not increase as the sample size grows [32]. This is in line with the discussion of stationarity in Section 2.2 and reveals the possible source for nonsense regression concerned with the spatial unit root SAR series [34] (p. 303).
Mur and Trívez [32] focus on the spurious regression when a spatial deterministic trend exists and show that under such circumstances similar symptoms related to unit roots occur. Consider the DGP
Y n = δ ι n + λ W n Y n + ϵ n .
Y n = ( I n λ W n ) 1 ( δ ι n + ϵ n )
= δ ι n + S n 1 ϵ n ,
where ι n is an n × 1 unit vector and ι n = S n 1 ι n . Comparing the term δ ι n in (22) with the time trends in a typical time series model, “ ι n ” is similar to the time trend “t”: t is different in terms of the relative position in time and the element of ι n is different in terms of its relative position in space. Also, the presence of such a trend term in the SAR process leads to spurious regression. Consider the simple regression
Y 1 n = α + Y 2 n β + μ n ,
where Y 1 n and Y 2 n are unrelated SAR processes generated by (20), respectively
Y 1 n = δ 1 ι n + λ 1 W n 1 Y 1 n + ϵ 1 n Y 1 n = δ 1 S 1 n 1 ι n + S 1 n 1 ϵ 1 n , Y 2 n = δ 2 ι n + λ 2 W n 2 Y 2 n + ϵ 2 n Y 2 n = δ 2 S 2 n 1 ι n + S 2 n 1 ϵ 2 n .
Assuming ϵ 1 n and ϵ 2 n are independent white noises, we expect the estimate of β in (23) to be 0. However, this is generally not the case which means spurious regression occurs. This can be seen from the fact that the correlation coefficient between Y 1 n and Y 2 n given in Mur and Trívez [32]
r Y 1 n , Y 2 n = r Y 1 n , r Y ¯ 1 n Y 2 n , r Y ¯ 2 n r Y 1 n , r Y ¯ 1 n 2 r Y 2 n , r Y ¯ 2 n 2 1 , n .
Though Fingleton argues that a high value of Moran’s I statistics may be a good indicator for the existence of spatial unit roots and spatial cointegration, he cannot distinguish between them, or even from the (genuine) positive spatial autocorrelation case. Some testing methods are developed and summarized in Section 5. The trend SAR series proposed by Mur and Trívez [32] seems to receive less attention, which may be due to the fact that when the mixed SAR process contains only a constant exogenous variable and W n is row-normalized, multicollinearity occurs; see Kelejian and Prucha [7] (p. 105), Lee [35] (p. 258), Lee [5] (p. 1907).

3.3. Spurious Regression under the near Unit Roots with a Row-Normalized, Circular Weight Matrix

The spurious regression considered in the previous two sections is under a row-normalized, noncircular weight matrix, which implicitly assumes an unconnected central unit in the spatial system, and is too restrictive to be used in empirical applications. Thus, Lee and Yu [34] study the spurious regression under a circular, row-normalized spatial weight matrix. The DGP process for the (mixed) SAR series is
Y j n = λ j n W j n Y j n + Z j n γ j + ϵ j n , j = 1 , , m .
In this case, λ j n 1 , since the unit roots are singular points of W j n . They study the consequence when λ n approaches 1, namely
λ n 0 = 1 1 ψ n ,
where ψ n as n .

3.3.1. Decomposition of Y n

Though the variance of Y n explodes as λ n 1 , Y n can be decomposed into a stable part and an unstable part by the decomposition of the weight matrix W n . This decomposition given in Lee and Yu [15] is used in Yu et al. [10,36], Yu and Lee [16] to study unit roots in a spatial dynamic panel data (SDPD) model. Because of its importance, this procedure is summarized here. See Baltagi et al. [37] and Lee and Yu [15,34] for more information.
Theorem 3. 
Suppose that W n is a row-normalized weight matrix from a symmetric matrix C n , i.e., W n = Λ n 1 C n , where Λ n is a diagonal matrix with its diagonal elements formed by the row sums of C n . Then (i) the eigenvalues of W n are all real; and (ii) W n is diagonalizable.
(i) can be easily seen from the fact that all symmetric matrices have real eigenvalues. For (ii)
W n = Λ n 1 C n Λ n 1 2 W n Λ n 1 2 = Λ n 1 2 C n Λ n 1 2 .
Let D n * be the eigenvalue matrix of Λ n 1 2 C n Λ n 1 2 , and R n * be the corresponding orthogonal eigenvector matrix, i.e., Λ n 1 2 C n Λ n 1 2 = R n * D n * R n * . Lee and Yu [15] show W n can be expressed as:
W n = Λ n 1 2 R n * D n * R n * Λ n 1 2 = Λ n 1 2 R n * D n * R n * 1 Λ n 1 2 = Λ n 1 2 R n * D n * Λ n 1 2 R n * 1 .
Let R n = Λ n 1 2 R n * , D n = D n * . By definition, R n is the eigenvector of W n and D n = D n * is the corresponding eigenvalue, so the eigendecomposition of W n is
W n = R n D n R n 1 .
Moreover, the largest eigenvalues of a row-normalized matrix are 1 in absolute value. Without loss of generality, Lee and Yu [15] assume there are m n eigenvalues equal to 1 and let
D n = diag 1 m n , d n , m n + 1 , , d n n ,
where 1 m n is 1 × m n vector of ones and d n i 1 for all i. So the eigenvalue matrix D n can be decomposed into two parts:
D n = J n + D ˜ n ,
where J n = diag 1 m n , 0 , , 0 and D ˜ n = diag 0 , , 0 , d n , m n + 1 , , d n n . Accordingly,
W n = R n D n R n 1 = R n J n + D ˜ n R n 1 = R n J n R n 1 + R n D ˜ n R n 1 = W n u + W ˜ n ,
where W n u = R n J n R n 1 and W ˜ n = R n D ˜ n R n 1 .
Lee and Yu [15] note that J n J n = J n and J n D ˜ n = 0 , so
W n u W n u = W n u , W n u W ˜ n = 0 , W n W n u = W n u .
Denoting S n ( λ ) = I n λ W n = R n I n λ D n R n 1 , and S n = S n ( λ n 0 ) where λ n 0 is the true value of λ , they obtain S n 1 ( λ ) = R n I n λ D n 1 R n 1 and thus
I n λ n 0 D n = diag 1 λ n 0 , , 1 λ n 0 , 1 λ n 0 d n , m n + 1 , , 1 λ n 0 d n n ,
I n λ n 0 D n 1 = diag ψ n , , ψ n , 1 1 λ n 0 d n , m n + 1 , , 1 1 λ n 0 d n n ,
since λ n 0 = 1 1 ψ n . Similarly,
I n λ n 0 D ˜ n 1 = diag 1 m n , 1 1 λ n 0 d n , m n + 1 , , 1 1 λ n 0 d n n .
Comparing the first m n diagonals of I n λ n 0 D n 1 and I n λ n 0 D ˜ n 1 , it is easy to obtain I n λ n 0 D n 1 = ψ n λ n 0 J n + I n λ n 0 D ˜ n 1 , thus
S n 1 λ n 0 = ψ n λ n 0 W n u + I n λ n 0 W ˜ n 1 .
Denote G n = W n S n 1 , by (38) and (34):
G n = ψ n λ n 0 W n u + W n I n λ n 0 W ˜ n 1 .
Thus, Lee and Yu [15] decompose Y j n = S j n 1 ( λ j n ) ( Z j n γ j + ϵ j n ) into two parts by (38):
Y j n = ψ j n Y j n u + Y ˜ j n ,
where Y j n u = λ j n W j n u Z j n γ j + ϵ j n and Y ˜ j n = I n λ j n W ˜ j n 1 Z j n γ j + ϵ j n .
It can be seen from (38) that when λ n 0 1 , S n 1 λ n 0 is ill-conditioned because ψ n and hence the variance of Y n explodes. This is caused by the first unstable term in (40) (the second term is stable). Thus, Y j n is of order ψ j n , which may grow too fast to yield useful asymptotic analysis. Thus, a rate-adjusted factor 1 ψ j n is needed to maintain a controllable rate. A similar idea applies to QMLE and 2SLS methods for estimation as shown later.

3.3.2. Spurious Regression of OLS under Near Unit Root

To study spurious regression, following Fingleton [14], Lee and Yu [34] consider the DGP that is similar to (26) but without exogenous variables:
Y j n = λ j n W j n Y j n + ϵ j n , j = 1 , , m .
Denote Y 1 , n = Y 2 n , , Y m n , X n = ι n , Y 1 , n where ι n is the n × 1 vector of ones. Let λ be a scalar, β = β 2 , , β m be an ( m 1 ) × 1 vector and δ λ , β . OLS, which may yield spurious regression, is then
Y 1 n = α ι n + Y 1 , n β + V n = X n δ + V n ,
where V n is an n × 1 vector of disturbances. To make sure the variable of interest is under controllable order, scale Y j n and S j n 1 S j n 1 as
Y j n * = 1 ψ j n Y j n , S j n = 1 ψ j n 2 S j n 1 S j n 1 , j = 1 , , m .
Lee and Yu [34] introduce a sufficient condition for Assumption 4 that ensures the UB of S n 1 , by (38):
Assumption 5. 
I n λ j n W ˜ j n 1 and W j n u are U B .
Under Assumption 5, S j n is UB. To study the properties of OLS estimation, it is sufficient to show the asymptotic behaviors of 1 n X n * Y 1 n * and 1 n Y i n * Y j n * , where X n * = ι n , Y 1 , n * . Proofs of these properties are about orders of matrices and random vectors as well as first and second moments of quadratic forms (some useful lemmas can be found at https://www.asc.ohio-state.edu/lee.1777/wp/sar-qml-r-appen-04feb.pdf, accessed on 28 March 2024). The OLS estimates for α and β : α ^ n , β ^ n = X n X n 1 X n Y 1 n can be expressed in terms of 1 n X n * Y 1 n * and 1 n X n * X n *
n α ^ n β ^ n = ψ 1 n Y m 1 1 n X n * X n * 1 1 n X n * Y 1 n * ,
where
Y m = 1 0 0 Y 2 m and Y 2 m = ψ 2 n 0 0 0 ψ 3 n 0 0 0 ψ m n .
Lee and Yu [34] notice that the scaling factor Y m is needed because the columns of X n have different orders ( 1 , ψ 2 n , , ψ m n ) and Y 1 n is of order ψ 1 n . Based on this fact, they give the following result:
1 n X n * X n * = D n , x x * + O p 1 n ,
where D n , x x * diag 1 , σ 2 2 1 n tr S 2 n , , σ m 2 1 n tr S m n , O p ( 1 n ) means the remaing terms of 1 n X n * X n * are at most of order 1 n and S j n symbols are defined in (43), for all j = 2 , , m .
Moreover, 1 n X n * Y 1 n * has limiting variance matrix:
Σ m * = σ 1 2 diag lim n 1 n l n S 1 n l n , σ 2 2 lim n 1 n tr S 2 n S 1 n , , σ m 2 lim n 1 n tr ( S m n S 1 n ) .
Lee and Yu [34] also adjust Assumption 2 to ensure 1 n X n * X n * is of full rank when n and obtain the asymptotic distributions of the OLS estimators:
Assumption 6. 
lim n 1 n tr S j n 0 for j = 1 , 2 , , m .
Denoting
D x x * lim n D n , x x * = diag 1 , σ 2 2 lim n 1 n tr S 2 n , , σ m 2 lim n 1 n t r S m n ,
then
n ψ 1 n Y m λ ^ n , β ^ n d N 0 , D x x * 1 Σ m * D x x * 1 .
Since Y 1 , n is independent of Y 1 n , one may expect insignificant β ^ in (42). However, whether β ^ j converges to 0 in probability or not depends on the factor n ψ 1 n ψ j n for 2 j m (note that Y m in (49) is diagonal): if n ψ 1 n ψ j n , β ^ j n p 0 is n ψ 1 n ψ j n -consistent; if n ψ 1 n ψ j n c < , β ^ j n is asymptotically normal because its limiting variance does not converge to 0; if n ψ 1 n ψ j n , β ^ j n is not bounded in probability and diverges. Intuitively, spurious regression will not occur if ψ 1 n approaches faster than ψ j n , or equivalently, λ 1 n approaches 1 more quickly than λ j n [34].

3.3.3. Other Test Statistics

It is also important to discuss the statistical properties of other test statistics based on (42), which could be potentially useful for distinguishing spurious regression. Lee and Yu [34] give the following theorem as a prerequisite:
Theorem 4. 
Under Assumptions 1 and 4, for any nonstochastic UB square matrix B n ,
1 n Y 1 n * B n P n Y 1 n * = O p 1 n .
where P n = X n X n X n 1 X n is the projection matrix of X n
Based on this theorem, Lee and Yu [34] show the order of the estimated variance of the disturbances σ ^ n 2 is
1 ψ 1 n 2 σ ^ n 2 = 1 ψ 1 n 2 e n e n n m = σ 1 2 tr S 1 n + O p 1 n ,
and
t β j d N 0 , lim n n · tr S j n S 1 n tr S j n · tr S 1 n , for every 2 j m ,
F = d 1 ( m 1 ) U m U m ,
where U m is an m 1 random vector with its elements u j m N 0 , lim n n tr S j n S 1 n tr S j n tr S 1 n . Thus t β j does not have a familiar asymptotic standard normal distribution, and the F-statistic has no familiar χ 2 ( m 1 ) distribution either.
Even though the t- and F-statistic are not reliable, Lee and Yu [34] suggest the combination of R 2 and Moran’s I could be a good indicator for spurious regression under near unit roots. Let M n 0 = I n 1 n l n l n , where l n is an n × 1 vector of ones, the coefficient of determination R 2 is
R 2 : = 1 e n e n Y 1 n M n 0 Y 1 n = 1 e n e n / n ψ 1 n 2 Y 1 n M n 0 Y 1 n / n ψ 1 n 2 = O p 1 n p 0 .
And the Moran’s I statistic is
I Moran : = 1 n e n W n e n 1 n e n e n = tr S 1 n * 1 W n S 1 n * 1 tr S 1 n * 1 S 1 n * 1 + O p 1 n p 1 , when W n = W 1 n .

3.3.4. Constant Terms in the DGP of Y j n ’s

Lee and Yu [34] also study the constant term and unit roots at the same time. It could be shown that the estimation of β is the same as in (49) after reparameterization. Consider the DGP of series Y j n c with near unit roots and a constant term c j n as
Y j n c = λ j n W j n Y j n c + c j n ι n + ϵ j n , for every j = 1 , , m .
Regress Y 1 n c on ι n and Y 1 , n c :
Y 1 n c = α c ι n + Y 1 , n c β c + V n c .
where α c is a constant and Y 1 , n c is similarly defined as Y 1 above. Since W j n is row-normalized, W j n ι n = ι n , we have Y j n c = S j n 1 c j n ι n + ϵ j n = ψ j n c j n ι n + S j n 1 ϵ j n = ψ j n c j n ι n + Y j n . So (57) could be rewritten as:
Y 1 n = α c ψ 1 n c 1 n + j = 2 m ψ j n c j n β j n c ι n + Y 1 , n β c + V n c .
Compare (58) and (42), OLS estimation of β c is the same as that of β and thus the corresponding statistics would be the same.

3.4. “Spurious” Regression with Equal Weights

Baltagi and Liu [38] show that under the special case where the spatial weight matrix is row-normalized and with equal spatial weights, i.e.,
w n , i j = 0 , if i = j , 1 n 1 , if i j .
spurious regression will not occur. This spatial weight matrix “is naturally suggested if all units are neighbors to each other and there is no other natural or observable measure of distance [39]”, such as interactions between students in a class or workers in a firm, etc. Without loss of generality, the DGP is assumed to be
Y n = λ 1 W n Y n + ϵ 1 n , X n = λ 2 W n X n + ϵ 2 n .
Consider the regression
Y n = α ι n + X n β + V n .
By the Frisch–Waugh–Lovell Theorem, the OLS estimation of β is β ^ = X n E n X n 1 X n E n Y n , where E n = I n J n n and J n is an n × n matrix of ones. Kelejian and Prucha [39] show the inverse of the matrix S n k = I n λ k W n , k = 1 , 2 , is
S n k 1 = δ k 1 J n δ k 2 I n ,
where δ k 1 = λ k n 1 + λ k 1 λ k and δ k 2 = n 1 n 1 + λ k . Using the fact that E n J n = 0 , Baltagi and Liu [38] show 1 n X n E n X n p σ 2 2 and 1 n X n E n Y n d N 0 , σ 1 2 σ 2 2 , so the asymptotic distribution of β ^ is given by
n β ^ = 1 n X n E n X n 1 1 n X n E n Y n d N 0 , σ 1 2 / σ 2 2 .
The asymptotic distribution of β ^ does not depend on λ 1 and λ 2 and is n consistent (compared with (49)), which means that the spurious regression does not occur.

4. Estimation and Inference

The spurious regression studied in the previous section suggests that OLS is not a good estimator in the SAR model with spatial unit roots. QMLE and 2SLS are alternative methods of estimation. This section briefly reviews these estimators and their performance under near unit roots.
When the error term of the DGP has a spatial autoregressive structure, OLS and (feasible) generalized least squares ((F)GLS) estimators are consistent. The efficiency of the OLS estimator was considered by Krämer and Donninger [40], Tilke [41] for the symmetric spatial weight matrices, and generalized by Krämer and Baltagi [42] with a broader covariance matrix. But the symmetry of the weight matrix is too restrictive to be used in practice, so Martellosio [43] generalizes this to nonsymmetric weights matrices. The efficiency of the OLS estimator is defined as η : = tr var X β ^ G L S tr var X β ^ O L S . But these papers generally focus on the relationship between W n and X, for example, when the column space of W n is contained by the column space of X. Following Lee and Yu [34], Baltagi et al. [37] derive the asymptotic properties of OLS, (F)GLS of β and point out important differences from conventional theory based on stationary spatial error.
A special SAR model has also been considered in regular lattices under spatial unit roots with the form
Y k , l = i = 0 p 1 j = 0 p 2 α i , j Y k i , l j + ϵ k , l , α 0 , 0 = 0 .
The simplest case of this special SAR model is the doubly geometric spatial autoregressive process:
Y k , l = α Y k 1 , l + β Y k , l 1 α β Y k 1 , l 1 + ϵ k , l .
It is called “unilateral” because only the previous units have effects on the latter ones and have a lower triangular weight matrix. It can be considered as the combination of two autoregressive (AR) models [44]. This model has been widely used in the area of image processing, agriculture trials, digital filtering, and other different fields. Model (64) is unstable when either | α | 1 or | β | 1 because of the existence of spatial unit roots [45,46].
A more complicated special case of the unilateral model is
Y k , l = α Y k 1 , l + β Y k , l 1 + γ Y k 1 , l 1 + ϵ k , l .
Because of its simple form, the estimation and inference of (64) can be derived without too many assumptions, as will be seen shortly. Moreover, when spatial unit roots exist, the limit of the variance of Y k , l are analytically obtained in Baran [47] when the parameters ( α , β , γ ) are located on the interior, on the edges, and on the vertices of the domain of stability. Paulauskas [48] explicitly shows that the growth rate of the variance of Y is different in dimensions d = 2 , 3 and d = 4 . Though this approach studies spatial unit roots from a different angle than that of Fingleton [14], it points to some similarity as in a recent paper [49] that will be discussed in Section 6.
Another possible way to remedy the problem caused by spatial unit roots in the SAR, SEM model is to relax the compactness assumption. When some parameters approach the boundary of the parameter space, consistency of extremum estimators could be obtained with compact parameter spaces. Thus the compactness assumption is standard in spatial econometrics because of the existence existence of the singular point 1 ρ m a x , see proofs in Kelejian and Prucha [6] Lee [5], Gupta [50]. But such an assumption is also restrictive in the sense that if we choose an arbitrary compact set on the open parameter space, the true global optimizer may be exclusive, especially for near unit root cases. A recent paper by Liu et al. [51] generalizes Theorem 2.7 in Newey and McFadden [52] (p. 2133), which relaxes compactness when the objective function of an extremum estimator is concave and allows the non-stochastic objective functions to depend on the sample size n. This generalization is suitable for spatial econometrics models because the sample observations are usually in triangular arrays. (A triangular array is a doubly indexed sequence of numbers or polynomials. Each row of the array is only as long as the row’s index. For example, the ith row contains only i elements.) The consistency of the QMLE of the SAR model and the MLE of the SAR Tobit model are investigated. But a closed-form solution is not obtained. On the other hand, Gupta [53] proposes a Newton-step computational algorithm of QMLE for a large-parameter-size SAR model, which is free from the compactness assumptions. Under the normality assumption, it has the same asymptotic efficiency as MLE, but has a closed-form solution and is computationally simple.

4.1. QMLE and 2SLS Methods for the (Mixed) SAR Model

4.1.1. Quasi-Maximum Likelihood Estimation Method

Lee [5] investigates the asymptotic distribution of the QMLE estimator of the mixed SAR model, which is the starting point for further analysis when spatial unit roots exist. This analysis is based on the discussion of the singularity of the information matrix of the log-likelihood function. Especially when the information matrix is singular, a scaling factor 1 h n is needed, where 1 h n is the order of the elements of the spatial weight matrix W n and thus 1 h n is the order elements of G n = W n S n 1 . We have seen in (39) that the order of elements of G n is ψ n in the near unit roots case; thus, a similar scaling factor will be needed.
The (mixed) SAR model under consideration is
Y n = X n β + λ W n Y n + V n ,
with its reduced form
Y n = S n 1 X n β 0 + V n = X n β 0 + λ 0 G n X n β 0 + S n 1 V n ,
since I n + λ 0 G n = S n 1 .
Also, Lee [5] imposes a weaker assumption about the spatial weight matrix and derives the information matrix.
Assumption 7. 
The elements w n , i j of W n are at most of order h n 1 , denoted by O 1 / h n , uniformly in all i , j where the rate sequence h n can be bounded or divergent. The ratio h n / n 0 as n goes to infinity.
Let V n ( δ ) = Y n X n β λ W n Y n , where δ = β , λ , then the log-likelihood function of (67) is
ln L n ( θ ) = n 2 ln ( 2 π ) n 2 ln σ 2 + ln S n ( λ ) 1 2 σ 2 V n ( δ ) V n ( δ ) ,
where θ = β , λ , σ 2 . The information matrix is
E 1 n ln L n θ 0 θ 1 n ln L n θ 0 θ = E 1 n 2 ln L n θ 0 θ θ + Ω θ , n ,
where
E 1 n 2 ln L n θ 0 θ θ = 1 σ 0 2 n X n X n 1 σ 0 2 n X n G n X n β 0 0 1 σ 0 2 n G n X n β 0 X n 1 σ 0 2 n G n X n β 0 G n X n β 0 + 1 n tr G n s G n 1 σ 0 2 n tr G n 0 1 σ 0 2 n tr G n 1 2 σ 0 4 .
with G n s = G n + G n . The existence of the extra Ω θ , n is because V n is not necessarily normally distributed.
To ensure the asymptotic distribution of QMLE θ ^ exists, Σ θ = lim n E 1 n 2 ln   L n θ 0 θ θ must be well defined. Lee [5] proves the nonsingularity of Σ θ can be guaranteed by the fact that there does not exist a nonzero vector λ = λ 1 , λ 2 , λ 3 such that a linear combination of columns of Σ θ is 0. This condition could be simplified as: there does not exist a λ 2 0 , such that
lim n 1 n σ 0 2 G n X n β 0 M n G n X n β 0 + lim n 1 n tr G n G n + tr G n 2 2 n tr 2 G n λ 2 = 0 .
Since each term in (72) is greater or equal to 0 (the first term is non-negative because it is symmetric; for the second term, tr G n G n + tr G n 2 2 n tr 2 G n = 1 2 tr C n + C n C n + C n = 1 2 t r ( C n s C n s ) 0 where C n = G n t r ( G n ) n I n and C n s = C n + C n ), Lee [5] studies the singularity of the information matrix in terms of these two terms, respectively. For the first term, by the partition matrix formula, lim n 1 n X n , G n X n β 0 X n , G n X n β 0 is nonsingular if and only if lim n 1 n X n X n and lim n 1 n G n X n β 0 M n G n X n β 0 are nonsingular. Moreover, under Assumption 7, if G n X n β 0 and X n are independent, one sufficient condition for the nonsingularity of Σ θ could be:
Assumption 8. 
lim n 1 n X n , G n X n β 0 X n , G n X n β 0 exists and is nonsingular.
However, Lee [5] states that if G n X n β 0 and X n are linearly dependent, for example, when W n is row-normalized and the relevant regressor is only the constant term, Assumption 8 should be adjusted to guarantee the second term in (72) is greater than 0:
Assumption 9. 
lim n 1 n G n X n β 0 M n G n X n β 0 = 0 and the h n is a bounded sequence and, for any λ λ 0 ,
lim n 1 n ln σ 0 2 S n 1 S n 1 1 n ln σ n 2 ( λ ) S n 1 ( λ ) S n 1 ( λ ) 0 .
Then under Assumption 7 and either 8 or 9, the asymptotic distribution of QMLE θ ^ will be
n θ ^ n θ 0 D N 0 , Σ θ 1 + Σ θ 1 Ω θ Σ θ 1 ,
where Ω θ = lim n Ω θ , n , Σ θ = lim n E 1 n 2 ln   L n θ 0 θ θ .
The above results are based on Σ θ being invertible of which a necessary condition is that { h n } is a bounded sequence. However, when lim n h n = , Σ θ will become singular because 1 n tr C n + C n C n + C n = O 1 h n = o ( 1 ) . Also, the singularity of the information matrix implies that the score function will be too flat to be useful and thus an adjustment of the rate should be imposed as in Lee [5]
Assumption 10. 
The h n is a divergent sequence, elements of M n G n X n β 0 have the uniform order O 1 h n , and lim n h n n G n X n β 0 M n G n X n β 0 = c with 0 c < . Under this situation, either (a) c > 0 , or (b) if c = 0
lim n h n n ln σ 0 2 S n 1 S n 1 h n n ln σ n 2 ( λ ) S n 1 ( λ ) S n 1 ( λ ) 0 ,
whenever λ λ 0 .
Lee [5] gives the asymptotic distributions of the QMLE under this rate-adjusted assumption:
n h n λ ^ n λ 0 D N 0 , σ λ 2 , n h n β ^ n β 0 D N 0 , σ λ 2 lim n X n X n 1 X n G n X n β 0 G n X n β 0 X n X n X n 1 , n σ ^ n 2 σ 0 2 = 1 n i = 1 n v i 2 σ 0 2 + o P ( 1 ) D N 0 , μ 4 σ 0 4 .

4.1.2. Generalized Spatial Two-Stage Least Squares (GS2SLS) Method

Kelejian and Prucha [7] proposed the GS2SLS method for the “cross-sectional (first-order) autoregressive spatial model with (first-order) autoregressive disturbances”:
Y n = X n β + λ W n Y n + u n | λ | < 1 = Z n δ + u n , u n = ρ M n u n + ϵ n , | ρ | < 1 ,
where Z n = X n , W n Y n and δ = ( β , λ ) . Since W n Y n is endogenous, it should be instrumented. Assume ρ is known (a consistent estimator is given by Kelejian and Prucha [6]), a Cochrane–Orcutt (CO)-type transformation applied to (77) yield the transformed regression
Y n * = Z n * δ + ϵ n = I n ρ M n Z n δ + ϵ n = Z n ρ M n Z n δ + ϵ n .
M n Z n should be instrumented. The ideal instruments are of course
E Z n = X n , W n E Y n and E M n Z n = M n X n , M n W n E Y n .
Then, by (10), E Y n = I n λ W n 1 X n β = i = 0 λ i W n i X n β , the ideal instruments are:
H n = X n , W n 1 X n , W n 2 X n , , M n X n , M n W n 1 X n , M n W n 2 X n , .
In practice, H n = X n , W n 1 X n , W n 2 X n , M n X n , M n W n 1 X n , M n W n 2 X n is used.
Assumptions about the instrument matrix are made to ensure that lim n 1 1 n H n Z n exists and is of full rank. In the near unit roots case, similar to the QMLE method, scaling factors are needed to guarantee this property as we will see later. With these assumptions, the GS2SLS procedure has three steps, as in Kelejian and Prucha [7], as follows:
  • Run 2SLS on Y n = Z n δ + u n with instruments H n . This yields δ ˜ n = Z ^ n Z ^ n 1 Z ^ n Y n , where Z ^ n = P H n Z n = X n , W n Y n ^ , P H n is the projection matrix of H n , and δ ˜ n = δ + O p ( n 1 2 ) .
  • Estimate ρ by Kelejian and Prucha [6] according to the GMM system: ρ n = G n λ + V n where λ = ρ , ρ 2 , σ 2 , then solve λ ˜ = G n g n , or by λ = argmin ρ , σ 2 V n V n . Both λ ˜ and λ are consistent, but λ is more efficient.
  • Assuming ρ is known, run 2SLS on the CO transformed regression (78) with instruments H n yields δ ^ n = Z ^ n * Z ^ n * 1 Z ^ n * Y n * , where
    Z ^ n * = P H n Z n * = P H n X n ρ M n X n , W n Y n ρ M n W n Y n = X n ρ M n X n , W n Y n ρ M n W n Y n ^ .
    By replacing ρ by its consistent estimation ρ ^ n (in Step 2). The feasible 2SLS estimator is
    δ ^ F , n = Z ^ n * ( ρ ^ n ) Z ^ n * ( ρ ^ n ) 1 Z ^ n * ( ρ ^ n ) Y n ( ρ ^ n ) .
Obviously, this procedure is for a SARAR model. If it is a SAR model, only step 1 is needed and H n = X n , W n 1 X n , W n 2 X n , .

4.1.3. Best Generalized Spatial Two Stage Least Squares (BGS2SLS) Estimators

Lee [8] does not drop the higher-order terms in (80) but use the fact that W n E Y n = W n I n λ W n 1 X n β = G n X n β by the definition of G n and proposes the best instrument:
H n * = I n ρ M n X n , W n I n λ W n 1 X n β = I n ρ M n X n , G n X n β .
With the corresponding simple instrumental variable estimator, the BGS2SLS is
δ ^ B , n = H n * Z n 1 H n * Y n .
If there is no SAR structure in the disturbance term, the best instrument is
H n * = X n , G n X n β ,
and β in (85) could be replaced by any consistent estimator such as the KP-GS2SLS estimator.
Compared with (80), H n * does not drop off the higher-order terms and is numerically equivalent to the ideal instrument, which in turn yields asymptotically optimal instrumental variable estimators.

4.2. Near Unit Roots in the SAR Model

Lee and Yu [15] study the asymptotic distribution of QMLE and 2SLS estiamtors of the SAR model by decomposing Y n (see Section 3.3.1). The model is given as
Y n = λ n 0 W n Y n + X n β 0 + ϵ n ,
where λ n 0 = 1 1 ψ n . And the reduced form, again, is
Y n = S n 1 X n β 0 + V n = X n β 0 + λ 0 G n X n β 0 + S n 1 V n .

4.2.1. QMLE

Obviously, the generated regressor G n X n β 0 is explosive because of ψ n in G n , see (39). This is very similar to the case in Assumption 10 (a) with lim n h n = in Section 4.1.1. This implies that the convergence rate of estimators of ρ n 0 is not the usual n case as in (76).
Additional assumptions are made in Lee and Yu [15] to control the magnitude of the unstable part of W n (Assumption 11), and specify the identification condition (Assumption 12), which are adjusted from Assumption 8:
Assumption 11. 
(1) lim n m n n > 0 ; (2) lim n n 1 tr W n u W n u > 0 ; (3) for any finite constant c , lim n n 1 tr I n c W n u I n c W n u > 0 .
Assumption 12. 
β 0 lim n n 1 X n W n u M n W n u X n β 0 > 0 holds.
Assumption 11 (1) (2) guarantee that m n is not too small compared to n and (3) implies that it is not too large. Assumption 12 is equivalent to (see Lee and Yu [15] (p. 338, Lemma 1 (7))) lim n 1 n ψ n 2 G n X n β 0 M n G n X n β 0 , which is a rate-adjusted version of Assumption 8, that ensures the identification uniqueness and implies nonsingularity of Σ θ . For a detailed discussion of these two assumptions, see Baltagi et al. [37] (p. 6). Also, since the adjusted rate is 1 ψ n 2 , QMLE λ ^ n would be n ψ n -consistent, see (76) derived under Assumption 10.
The information matrix will be the same as in (70). In Section 3.3.2, when studying the spurious regression of OLS, we mentioned that the scaling factor Y m is needed in terms of the order of X n . Here, for the QMLE, similar scaling factor will also be introduced because elements of Σ n and Ω n have different orders as the existence of G n . The second column and row of Σ n and Ω n , which are the derivatives with respect to the spatial coefficient λ n contain G n , thus, they have to be scaled by a factor 1 ψ n . Specifically, the ( 2 , 2 ) element should be scaled by a factor 1 ψ n 2 . This can be done by a left and right multiplying matrix Y θ , n 1 , where
Y θ , n = Y δ , n 0 ( k + 1 ) × 1 0 1 × ( k + 1 ) 1 and Y δ , n = I k 0 k × 1 0 1 × k ψ n .
Thus, Lee and Yu [15] give
n Y θ , n θ ^ n θ n 0 = Y θ , n 1 1 n 2 ln L n θ ˜ n θ θ Y θ , n 1 1 Y θ , n 1 1 n ln L n θ n 0 θ .
Let Σ = lim n Y θ , n 1 Σ n Y θ , n 1 and Ω = lim n Y θ , n 1 Ω n Y θ , n 1 and assume they exist, the asymptotic distribution of θ ^ n is
n Y θ , n θ ^ n θ n 0 D N 0 , Σ 1 + Σ 1 Ω Σ 1 .
Recently, Rossi and Lieberman [54] combine the near unit roots with a similarity-based weighted matrix and study the consistency of the QMLE estimator when the spatial coefficient 0 λ < 1 and λ = 1 , by allowing uncentered units. The element of similarity-based weighted matrix is w i , j = s x i , x j ; w 0 j i s x i , x j ; w 0 , where s x i , x j ; w 0 is some function that measures the similarity between unit i and j according to some parameter w 0 . The parameters they are most interested in are θ 2 = ( λ , w ) . They establish the connection between the λ and the order of uniform absolute row-sum norm of S 0 1 , S 0 1 = S 1 θ 20 = O n γ . This means that λ 0 = 1 γ = 1 [54] (p. 11, Proposition 1). Recall that Var ( Y n ) = σ ϵ 2 S n 1 S n 1 , so when γ = 0 , the variance is independent of n, corresponding to the standard SAR setup ( λ < 1 and fixed). In the case 0 < γ < 1 , the variance increases with the sample size but with a lower speed, which is the case that we have seen when studying near unit root; but when γ = 1 , λ = 1 , the variance increases so fast that the non-standard limit distribution of θ = ( λ , w ) has to be established on a case-by-case basis, according to the resulting S n 1 . Their result is much more complicated than that of Lee and Yu [15] because of the introduction of similarity structure in the weight matrix, but are much more flexible since now W n is no longer fixed but data-driven, and is potentially more useful in empirical work.

4.2.2. GS2SLS and BGS2SLS

Lee and Yu [15] derive the 2SLS estimators and their asymptotic distributions using the procedures mentioned above. Using instruments defined in (80), the GS2SLS estimator of δ n 0 is
δ ^ n , 2 s l s = Z n P H n Z n 1 Z n P H n Y n .
Since Z n = X n , W n Y n contains W n , the adjustment by Y δ , n is needed. The asymptotic properties of the GS2SLS estimators of λ n 0 and β 0 are obtained as follows:
Y δ , n n δ ^ n , 2 s l s δ n 0 = n β ^ n , 2 s l s β 0 n ψ n λ ^ n , 2 s l s λ n 0 d N 0 , Φ 2 s l s
where Φ 2 s l s = σ 0 2 lim n 1 n Y δ , n 1 Z n P H n Z n Y δ , n 1 1 . So the GS2SLS estimator λ ^ n , 2 s l s of λ n 0 is n ψ n -consistent, which is higher than the usual n rate in Kelejian and Prucha [7] and Lee [8], but β ^ n , 2 s l s has the usual n rate of convergence.
Choosing the instrument H n * = X n , G n X n β 0 as in (85), the BGS2SLS estimator is δ ^ n , b 2 s l s = H n * Z n 1 H n * Y n with asymptotic distribution:
Y δ , n n δ ^ n , b 2 s l s δ n 0 d N 0 , Φ b 2 s l s
where Φ b 2 s l s = σ 0 2 lim n 1 n Y δ , n 1 H n * H n * Y δ , n 1 1 . Since Φ b 2 s l s Φ 2 s l s is negative semidefinite, the BGS2SLS estimator is more efficient.
The above result is based on the fact that G n X n β and X n are independent, which makes sure that the instrument matrix is of full rank, otherwise the 2SLS estimator will be inconsistent [9]. Liu [55] shows that even though G n X n β and X n are linearly dependent, i.e., G n X n β = X n c n , where c n is a nonzero vector, we still have ρ ^ 2 S L S ρ = O p 1 ψ n , β ^ k , 2 S L S β k = O p c k n ψ n + O p 1 n under near unit roots case; and ρ ^ 2 S L S ρ = O p ( 1 ) , β ^ k , 2 S L S β k = O p c k n + O p 1 n under the regular case, as long as c n = o ( 1 ) . This is equivalent to G n X n β and X n , which are asymptotically independent.
To provide guidelines for empirical studies, Lee and Yu [15] conduct simulations to compare the performance of QMLE and 2SLS methods. QMLEs are relatively robust whether the error term is normally distributed or not. Moreover, as n increases (and the spatial coefficient is closer to the spatial unit roots), the QMLEs perform better than the 2SLS estimators because of smaller variances. One interesting phenomenon is that the best 2SLS estimators are even worse than the regular 2SLS estimators in some cases, which violates the theoretical result as shown in (93). One possible reason for this is that the best 2SLS estimator requires an initial consistent estimator by construction (see (85)) and under spatial unit roots, such an initial estimator may not be accurately calculated.

4.3. Near Unit Roots in the SEM Model

Baltagi et al. [37] extend the study of near unit roots from SAR model to SEM model by considering the OLS, GLS and FGLS estimation and properties of the corresponding statistics. The model is given as
Y n = X n β 0 + u n , u n = λ n 0 W n u n + ϵ n ,
with λ n 0 = 1 1 ψ n . Similar to Lee and Yu [15], u n could be decomposed to
u n = S n 1 ϵ n = ψ n W n u ϵ n + I n λ n 0 W ˜ n 1 ϵ n .
One more assumption that they impose is
Assumption 13. 
The elements of X n are nonstochastic and bounded, uniformly in n, lim n n 1 X n X n exists and is nonsingular. lim n n 1 X n W n u W n u X n exists. Furthermore, lim n n 1 X n S n S n X n exists and is nonsingular.
The OLS estimator β ^ OLS = X n X n 1 X n Y n has the asymptotic distribution when ψ n n 0
n ψ n β ^ OLS β 0 d N 0 , σ 0 2 lim n 1 n X n X n 1 lim n 1 n X n W n u W n u X n lim n 1 n X n X n 1 ,
and when ψ n n c < , β ^ OLS β 0 = O p ( 1 ) . Thus β ^ OLS = β 0 + O p ψ n n , which is n ψ n -consistent and is slower than the stationary error term case.
Baltagi et al. [37] also study the asymptotic properties of the GLS and FGLS estimators. If λ n 0 is known, β ^ G L S = X n S n S n X n 1 X n S n S n Y n , and
n β ^ G L S β 0 d N 0 , σ 0 2 lim n 1 n X n S n S n X n 1 ,
which implies that β ^ G L S is robust for the near unit roots in the error term because it has n rate of convergence. The feasible GLS (FGLS) could be achieved by replacing λ by a consistent estimator λ ^ n , which yields
β ^ F G L S = X n S ^ n S ^ n X n 1 X n S ^ n S ^ n Y n ,
where S ^ n = I n λ ^ n W n . It can be seen that FGLS is identical to the QMLE: concentrated log likelihood function of (94) with respect to λ is
ln L n ( λ ) = n 2 ( ln 2 π + 1 ) n 2 ln σ ^ n 2 ( λ ) + ln S n ( λ ) ,
where
σ ^ n 2 ( λ ) = n 1 u n S n ( λ ) P ¯ n ( λ ) S n ( λ ) u n ,
P ¯ n ( λ ) = I n S n ( λ ) X n X n S n ( λ ) S n ( λ ) X n 1 X n S n ( λ ) .
The QMLE is of order ψ n λ ^ n λ n 0 = o p ( 1 ) , which is ψ n -consistent and
β ^ n λ ^ n = X n S n λ ^ n S n λ ^ n X n 1 X n S n λ ^ n S n λ ^ n Y n , σ ^ n 2 λ ^ n = 1 n u n S n λ ^ n P ¯ n λ ^ n S n λ ^ n u n .
Comparing with (98), β ^ n λ ^ n is a FGLS of β using λ ^ n . Thus, the QMLE β ^ M L and the infeasible GLS estimator β ^ G L S have the same asymptotic distribution as shown before. Next, Baltagi et al. [37] consider the Wald test statistic for the null hypothesis H 0 : R β = r for OLS, GLS and FGLS, where R is a q × k matrix of rank q < k and r is q × 1 . For OLS,
W OLS : = R β ^ O L S r σ ^ O L S 2 R X n X n 1 R 1 R β ^ O L S r d ξ σ 0 2 lim n 1 n tr W n u W n u R lim n 1 n X n X n 1 R 1 ξ ,
where
ξ N 0 , σ 0 2 R ( lim n n 1 X n X n ) 1 ( lim n n 1 X n W n u W n u X n ) ( lim n n 1 X n X n ) 1 R ,
does not have a standard χ 2 distribution, which is similar to the F-statistic shown above. However, the GLS Wald statistic
W G L S = R β ^ G L S r σ 0 2 R X S n S n X 1 R 1 R β ^ G L S r d χ k 2 ,
has a chi-squared limiting distribution.
Baltagi et al. [37] conduct extensive simulations. Using the root mean squared error (RMSE) as the evaluation criteria, the QML (FGLS) estimators perform uniformly better than the OLS estimator. In particular, when the spatial coefficient is sufficiently close to 1 and the sample size n increases, the RMSE of the OLS estimator grows dramatically. Together with the fact that the Wald test statistic based on the QML method has a standard Chi-squared distribution, QMLE is recommended when near spatial unit roots exist in the spatial error model.

4.4. Doubly Geometric Spatial Autoregressive Process

The main difference between the SAR and the doubly geometric spatial autoregressive models is that the spatial dependence form of the latter is clearly specified. However, for the SAR model, such dependence relies on the specification of W n whose explicit form varies in different situations.
For model (65), based on the observation X k , : 1 k m and 1 n , Baran [47] shows that the asymptotic normality of the estimators α ^ m , n , β ^ m , n is m n α ^ m , n α β ^ m , n β D N 0 , Σ α , β in the stable case ( | α | < 1 and | β | < 1 ), with some covariance matrix Σ α , β . For the unstable case ( α n 1 , β n 1 ), using the martingale central limit theorem, Bhattacharyya et al. [56,57] show “one step Gauss-Newton” estimators are asymptotically normal with convergence rate n 3 / 2 . This is different from the classical time series A R ( 1 ) , where the OLS estimator converges to a fraction of functionals of the standard Brownian motion: T a O L S 1 W ( 1 ) 2 1 2 0 1 W ( r ) 2 d r [58] (p. 281).
Baran and Pap [59] consider the more complicated model as in (66). The model is stable if and only if ( α , β , γ ) S , where S is the open tetrahedron with vertices V : = { ( 1 , 1 , 1 ) , ( 1 , 1 , 1 ) , ( 1 , 1 , 1 ) , ( 1 , 1 , 1 ) } . They also prove that the OLS estimator is asymptotically normally distributed with the convergence rate n when the model is stable, and n 3 / 2 otherwise. (The simpler model Y k , l = α Y k 1 , l + β Y k , l 1 + ϵ k , l , with possibly α = β was investigated in Baran et al. [60], Baran et al. [44], Baran and Pap [61] under stable and unstable cases. Under different settings, the limiting distribution of the OLS estimator is normal but has different rates of convergence.)
Roknossadati and Zarepour [62,63] study the limiting behavior of M-estimation for the near unit roots of model (65). The M-estimator α ^ n , β ^ n of α n , β n is defined to minimize of the objective function
g α n , β n = i = 2 n j = 2 n λ Y i j α n Y i 1 , j β n Y i , j 1 + α n β n Y i 1 , j 1 ,
for some convex function λ ( · ) . Roknossadati and Zarepour [62] show that the self-normalized M-estimators are asymptotically normal, and when the series is stable, the convergence rate of M-estimators is still n 3 / 2 , same as in Bhattacharyya et al. [56,57]. But if it is unstable, i.e., when the model has infinite variance innovations, the M-estimates have a higher consistency rate.

5. Tests for Spatial Unit Roots and Nonstationarity

Recognizing the possible consequences of spatial unit roots, it is necessary to test for it. In fact, in nonstationary cases, the estimator is inconsistent and diverges [64]. If the series contains spatial unit roots, one may employ the spatial first difference as recommended by Fingleton [14]: after the first-order difference, such series will be converted to a stationary one, otherwise it is over-differenced and spatial correlation still exists. Based on this idea, Lauridsen and Kosfeld [17,18] propose two-stage LM tests to check for spatial unit roots. However, such LM tests have a high power function because of L M > L R > W a l d in finite samples and are not useful for spatial cointegration since they mis-specify the regression in the second stage. A Wald test is proposed by Lauridsen and Kosfeld [65] but it does not have a usual χ 2 distribution so simulation has to be conducted before each test to obtain the critical values. A different approach introduced by Beenstock et al. [19] uses the fact that when spatial unit roots exist, the variance explodes and the spatial impulse does not die out as distance increases, so they iterate on the parameter space to find out the value of the unit roots (for irregular lattice) and then generate nonstationary series to conduct interval estimation.
Martellosio [66] derives the power properties of invariant tests, for example, u ^ Q u ^ u ^ u ^ , where u ^ is the OLS residuals and Q is a fixed matrix. When Q = W n , we obtain the Cliff–Ord test. When the regression contains only a constant, the Cliff–Ord test reduces to Moran’s test as introduced before, which is best locally invariant as shown by King [67]. It has been shown that for the SEM model, as λ 1 ρ m a x , the test power vanishes. For the SAR model, as λ 1 ρ m a x , the limiting power is either 0 or 1. Krämer [68] shows similar conclusions but focuses on the symmetric weight matrix. Martellosio [69] further shows the power of any test vanishes as spatial correlation increases for a set of regression spaces. Heteroskedasticity robust tests have been studied. For example, Born and Breitung [70], Baltagi and Yang [71] design diagnostic tests for SEM and SAR employing the outer product of gradients (OPG) variant of the LM test which are robust against heteroskedastic (and non-normal) errors. But these tests suffer from the same deficiency as in Martellosio [66] because such test is asymptotically equivalent to Moran’s I. Baltagi and Yang [72] have also shown that the standard LM test undergoes finite sample distortion when spatial dependence is heavy in both spatial and panel data settings. Recently Preinerstorfer [73] suggests some modified tests to avoid this “zero-power trap” phenomenon, which works well for small spatial autocorrelation, but still has limiting power smaller than 1 (only 0.619 by simulation). Thus, the invariant test of I ( 0 ) null hypothesis is not satisfactory when spatial unit roots exist, and methodologies to determine it (Tests of the I ( 1 ) null hypothesis) deserve more attention.

5.1. Two-Stage LM Test for the Sources of Spurious Spatial Regression

Lauridsen and Kosfeld [17] develop a two-stage LM test to distinguish between two possible sources for spurious regression. The first one is the existence of spatial (near) unit roots in the regressand and/or regressors as in Fingleton [14], Lee and Yu [34]; the second is that the spatial error term itself is nonstationary. So the LM tests are essentially testing if the spatial process is stable or not. The idea originates from the fact that Fingleton [14] suggests a high value of Moran’s I statistic as an indicator for both spatial nonstationarity and spurious regression, but we cannot distinguish between them or even distinguish between the nonstationarity and the positive spatial correlation among the error terms, which by definition imply a high value of Moran’s I. Specifically, we are trying to distinguish if (i) the X n is a SAR process and we regressed Y n on X n or (ii) the model itself is SEM as Y n = X n β + ϵ n , where ϵ n = λ ϵ W n ϵ n + μ n , because both (i) and (ii) can cause spatial autocorrelation.
There are at least three advantages of the LM tests [17]. First is that compared with Wald or LR, LM is usually simpler to compute because it is constructed under H 0 . Second is that with the LM test, it is possible to control for some omitted model features such as heterogeneity and autoregression, as in Anselin [74], which will be discussed later. The last one is that, other statistics may not have a standard asymptotic distribution, such as the OLS Wald type statistic as in Baltagi et al. [37]. The proposed two-stage LM test is based on the SEM model and all four possible results are summarized in Lauridsen and Kosfeld [17] (Table 1):
  • Under H 0 : λ ϵ = 0 , the LM error statistic (LME) developed by Anselin [74] (p. 11, Equation (35)) is
    L M E = e W n e / σ 2 2 tr W n 2 + W n W n χ 2 ( 1 ) .
    Thus large values of L M E reject the null hypothesis, which implies either 0 < λ ϵ < 1 or λ ϵ = 1 .
  • The next step is to test if H 0 : λ ϵ = 1 . This could be carried out by using the spatial differencing we introduced before. Under H 0 , Δ ϵ = μ , thus the first order difference on the regression, Δ Y n = Δ X n β + μ n , yields i.i.d. error μ n , which means the value of differenced LME (DLME) should be close 0 under H 0 . But if λ ϵ < 1 , Δ represents overdifferencing, i.e., Δ ϵ n = ( I n W n ) ( I n λ ϵ W n ) 1 μ n , so spatial correlation in the error term still exists, and we cannot reject H 0 .
Similar procedures could be used to investigate whether Y n or any X n are spatially nonstationary, as the case in Lee and Yu [34]. Letting Z n be one of Y n , X n 1 , X n 2 , , Lauridsen and Kosfeld [17] suggest using the regression
Z n = α ι n + ϵ n , Δ Z n = α Δ ι n + Δ ϵ n = Δ ϵ n ,
to obtain LME and DLME, respectively. Z n is regressed on a constant term because there is no meaningful regressor but we still need the residuals.
Spatial cointegration could also be tested using this LM test. Thus, after determining Y n and X n are nonstationary, regress Y n on X n and Δ Y n on Δ X n to obtain LME and DLME. The cointegration relation exists if LME is 0; and non-cointegration if LME is positive and DLME is 0; the limiting case of “near cointegration” occurs if LME and DLME are positive.
Moreover, Lauridsen and Kosfeld [18] generalize their two-stage LM test to account for unobserved heteroskedasticity. They specify the covariance matrix of ϵ n , Ω = diag { σ 1 2 , , σ n 2 } , have the diagonal element σ i 2 = f ( Z i , λ Z ) , where Z i is P × 1 vector of observations of exogenous variables for region i, related to σ i 2 via the P × 1 vector of parameters λ Z . So the statistic in (107) should be adjusted as in Anselin [74] (p. 9, Equation (29)):
L M E H = e W n e / σ 2 2 tr W n 2 + W n W n + f Z Z Z 1 Z f 2 χ 2 ( P + 1 ) ,
with f i = e i 2 σ 2 1 and Z as the n × P matrix containing the Z vectors that cause heteroskedasticity.
However, the Lauridsen and Kosfeld test procedure is not without problems. Beenstock et al. [19] point out that this procedure is not suitable for testing spatial cointegration since the second stage is misspecified. To see it more clearly, regress Δ Y n = β Δ X n + v n . The LM procedure asserts that if v n is not spatially correlated, then X n and Y n are spatially integrated; and if v n is spatially correlated, X n and Y n are not spatially integrated because of overdifferencing. Nevertheless, regressing Δ Y n on Δ X n is equivalent to regressing two white noise series, ϵ Y and ϵ X ; because Y and X are both I(1). Hence, the corresponding residuals must be not spatially correlated as long as ϵ Y and ϵ X are independent, regardless of whether v n is spatially correlated, nonstationary, or not.

5.2. A Wald Test for Spatial Nonstationarity

Lauridsen and Kosfeld [65] suggest a Wald post-test statistic. Based on MLE, under H 0 : R θ = q , the general form of the Wald test is W = ( R θ q ) R V R 1 ( R θ q ) , where V = I 1 is the inverse of information matrix. If we specify the null hypothesis as λ = 1 , then with R = 0 , 1 , 0 and q = 1 , we have W = ( λ 1 ) 2 V λ , where V λ is the diagonal element of V corresponding to λ . However, as mentioned before, Wald statistics may not have a standard distribution so simulations are conducted. Unlike Fingleton [14], to generate SAR series with spatial unit roots, Lauridsen and Kosfeld [65] do not introduce the noncircular matrix. Thus λ = 1 is a singular point of ( I n λ W n ) so the inverse does not exist. To solve this issue, they use the Moore–Penrose pseudoinverse.
According to Monte Carlo simulation, they find the critical limit of the Wald test under spatial nonstationarity is higher than the χ 2 ( 1 ) distribution, especially for the 5th and 10th percentile.

5.3. Test Unit Roots and Cointegration in the Sense of Spatial Impulses

Beenstock et al. [19] come up with an innovative method to test spatial unit roots and spatial cointegration by considering the behavior of the variance and the spatial impulse. Also, they do not assume unconnected spatial units or row normalize the spatial weight matrix either. Thus, based on the topology of the unit neighborhood, the spatial unit roots are λ * = 1 n in the regular lattices where n is the maximum and general number of neighbors of each unit. For example, n = 2 , for bilateral space, n = 4 for rook lattice and n = 8 for queen lattice, with λ * = 1 2 , 1 4 , 1 8 respectively. (“The weight matrix with first-order contiguity according to the rook criterion has the cells immediately above, below, to the right, and to the left, for a total of four neighboring cells. The weight matrix with first-order contiguity according to the queen criterion is eight cells immediately surrounding the central cell” [75] (p. 131). For the introduction of other types of the spatial weight matrices, see Kelejian and Robinson [26] (pp. 94–95).) And with spatial unit roots, ( I n λ * W n ) 1 is still well-defined because of the existence of the edge effect, that is, there exist some units having fewer neighbors than n. However, even though ( I n λ * W n ) 1 exists, the variance tends to explode even in finite sample space, which provides us with a way to determine the spatial unit roots for any arbitrary irregular lattices.

5.3.1. Spatial Impulse

The spatial impulse response is essentially the consequence of the shocks from one location to another. Intuitively, shocks should have no effect on the remote units if the spatial data are stationary. Beenstock et al. [19] first consider the simplest SAR model in lateral space:
Y n , j = λ Y n , j + 1 + Y n , j 1 + u n , j , j = , , , 1 λ 1 L + L 2 Y n , j = λ 1 u n , j 1 ,
where L denote a spatial lag operator such that L i Y n , j = Y n , j i . The auxiliary equation is
x 2 λ 1 x + 1 = 0 .
When the discriminant of the above equation is greater than 0, 0 < λ < 1 2 , and there are two different solutions, x 1 < 1 < x 2 , by Vieta’s formula. Hence, Beenstock et al. [19] express Y n , j as
Y n , j = λ 1 x 1 1 x 1 i = 1 x 1 i u n , j i + i = 0 x 1 i u n , j + i , Y n , j u n , j i = Y n , j u n , j + i = λ 1 x 1 i x 1 1 x 1 ,
where (112) is known as the Wold representation that expresses Y n , j in terms of the shocks. The impulse from location j i to i tends to 0 because x 1 < 1 . Also, x 1 varies with λ . When λ = λ * = 1 2 , x 1 = x 2 = 1 , so the impulse does not die out with distance and explodes. This fact can also be seen from
Var ( Y n , j ) = λ 2 1 + x 1 2 1 x 1 2 x 1 1 x 1 2 σ u 2 .
If 0 < λ < 1 2 , Var ( Y n , j ) is finite and independent of j; if λ = 1 2 , x 1 = 1 , and Var ( Y n , j ) is infinite.
For the bilateral space case, λ * = 1 n are the spatial unit roots as shown before. Because of the edge effect, the singular point is strictly greater but approaches λ * . This fact shows a downside of the row-normalized spatial weight matrix: it overstates the true weight of the unit at the edge of the lattice. For example, in a rook lattice, the units have three neighbors with weight 1 3 at the edge, and four neighbors with weight 1 4 in the center. The row-normalized procedure assigns a higher weight to the neighbors of edge units. This weight assignment is not necessarily reasonable and makes the spatial unit roots the same as the singular points. Moreover, without row-normalization, the edging units play the role of the unconnected unit as in Fingleton [14]. The general SAR model in bilateral space is Y n = λ W n Y n + u n and the Wold representation is Y n = A n u n , where A n = ( I n λ W n ) 1 = I n + i = 1 λ i W n i . Let the spatial impulse response be defined as d Y n , j d u n , j = a j j and d Y n , j d u n , i = a j i . Analytical solutions of spatial impulse response in bilateral and higher dimension lattices are not obtained, but Beenstock et al. [19] expect a j j to be positively related to the number of spatial units because of the larger spillover effect and a j i varies inversely with the distance between i and j in the stationary case. If spatial unit roots exist, as in the lateral case, the impulse a j i would not die out as the distance increases. This is supported by the simulation, though only the finite sample case could be simulated, see Beenstock and Felsenstein [76] (Figures 5.2 and 5.4). Compared with λ < λ * , when λ = λ * , it obviously shows a qualitative difference in the persistence of spatial impulses, as well as in the tendency for the explosion of the variance.
In the irregular lattices, the number of neighbors for the unit is undetermined generally, so the spatial unit roots, λ * , cannot be calculated as the reciprocal of n. However, since the nonstationarity implies that spatial impulses do not disappear, one can find the empirical spatial unit roots by simulation.
The simulation method in Beenstock et al. [19] to calculate the critical value is pretty flexible and can be adapted to different models. For example, when both dynamic and spatial terms are included, Beenstock and Felsenstein [23] develop a similar procedure for testing cointegration in nonstationary panel data when estimating the spatial spillover effect in housing construction for Israel.

5.3.2. Spatial Unit Roots and Cointegration Tests

Knowing the spatial unit roots λ * , Beenstock et al. [19] conduct Monte Carlo simulations to generate the artificial SAR series and use the MLE method to estimate SAC to obtain the corresponding distributions under different topologies (different sample size, criteria, etc.). Results show that the empirical distribution of SAC could be used to construct interval estimation and critical values for statistics that test spatial unit roots. For the spatial cointegration test, a similar procedure applies, but OLS estimation is used.

5.4. Some Applications

Kosfeld and Lauridsen [77] offer an application of the two-stage LM test in Section 5.1 to the income and productivity convergence in the German regional labor market. They find highly significant LME and DLME statistics (refer to formulas above like (107)) for all variables, which means the spatial unit roots are rejected. Yesilyurt and Elhorst [20] estimate the spatial interaction effects of inflation in Turkey. Because the regional inflation rates have a high tendency to co-move over time, they question whether the inflation rates of different regions are stationary in space. Using the two-stage LM statistics from (108), they find that the inflation curve is stationary in space. Olejnik [21] studies the income process of the extended European 25 based on the augmented Solow model taking into account the spatial autocorrelation effect. The stationarity of the error term as well as all variables in the model are investigated. No problem of spurious regression is found. Machado et al. [22] examine the spatial correlation of traffic accidents of vulnerable road users (such as pedestrians and cyclists) in big cities and detect the factors that contribute to these accidents. Because their study covers several cities, the model specifications may vary across different locations. Thus they use the two-stage LM statistic to choose the best model, see Machado et al. [22] (Table 4). Though the Wald post-test in Lauridsen and Kosfeld [65] is asymptotically equivalent to the LM test, “It is generally recommended to choose among these alternatives on the basis of computational ease [78] (p. 94)”.

6. Related Topics

Spatial panel data have been studied extensively. The spatial dependence is incorporated in the error component [12,75,79] or by spatial lag dependence [11,80]. See Baltagi [2] for a textbook discussion. Also, the panel data model can have time lagged dependent variables. If the panel data model includes both spatial and dynamic features, it is named as spatial dynamic panel data (SDPD) model by Yu et al. [10]. Yu et al. [10], Yu and Lee [16] and Yu et al. [36] study the QMLE estimator of the SDPD model under stable, unit roots, and spatial cointegration respectively. The concept of unit roots under the SDPD model is a combination of the spatial and dynamic one. To see this more clearly, Yu et al. [10] specify the model as
Y n t = λ 0 W n Y n t + γ 0 Y n , t 1 + ρ 0 W n Y n , t 1 + X n t β 0 + c n 0 + V n t ,
where Y n t = y 1 t , y 2 t , , y n t and V n t = v 1 t , v 2 t , , v n t are n × 1 column vectors. Since S n ( λ ) = I n λ W n , assuming S n is invertable, the reduced form is
Y n t = A n Y n , t 1 + S n 1 X n t β 0 + S n 1 c n 0 + α t 0 S n 1 ι n + S n 1 V n t ,
where A n = S n 1 γ 0 I n + ρ 0 W n . If the infinite sums are well-defined, then by continuous substitution
Y n t = h = 0 A n h S n 1 c n 0 + X n , t h β 0 + V n , t h .
So instead of focusing on the singular points of S n , A n should be considered, which contains λ , γ and ρ : the parameter of contemporaneous spatial effect, time lagged variable and time–spatial effect. A similar process as in Section 3.3.1, letting w ¯ n = diag w ¯ n 1 , w ¯ n 2 , , w ¯ n n be the eigenvalue matrix of W n , Yu et al. [10] show that the eigenvalue matrix of A n is D n = I n λ 0 w ¯ n 1 γ 0 I n + ρ 0 w ¯ n , which can be decomposed as D n = γ 0 + ρ 0 1 λ 0 J n + D ˜ n . The power matrix of A n follows as A n h = γ 0 + ρ 0 1 λ 0 h R n J n R n 1 + B n h with B n h = R n D ˜ n h R n 1 since the eigenvector matrix R n is orthogonal and J n D ˜ n = 0 . Thus, whether Y n t is stable or not depends on the value of γ 0 + ρ 0 1 λ 0 compared with 1. Consequently, the decomposition of Y n t , which is a generalization of (40), can be expressed as
Y n t = Y n t u + Y n t s + Y n t α ,
where Y n t s is a possible stable part, Y n t u is a possible unstable part, and Y n t α is the time effect part, see Lee and Yu [81] for details. A data transformation procedure is imposed by them to eliminate both the time effects and the possible unstable term. Based on their analysis, the eigenvalues of A n , the asymptotic properties of QMLE and bias are derived. When eigenvalues of A n are all less than 1 ( γ 0 + λ 0 + ρ 0 < 1 ), or some equal to 1 ( γ 0 + λ 0 + ρ 0 = 1 and γ 0 1 ), or all equal to 1 ( γ 0 + λ 0 + ρ 0 = 1 and γ 0 = 1 ), the information matrix has different properties, see Yu and Lee [16] (Table 5).
Thus, the test for the unit eigenvalues of A n is of great importance. Most attention is paid to the unit roots in the time dimension, i.e., γ 0 = 1 and equivalently, if λ 0 + ρ 0 = 0 . Unit root tests in panel data under spatial dependence have been extensively studied, see Baltagi [2] (Section 12.3) for a summary. Also, the performance of different tests has been considered in Baltagi et al. [24]. The test for H 0 : γ 0 + ρ 0 + λ 0 = 1 has been investigated in Lee and Yu [81] (Section 14.3.4). Such a test works well when λ 0 + γ 0 + ρ 0 < 1 . However, when λ 0 + γ 0 + ρ 0 > 1 and T is small, it is not reliable. Thus, further study of the unit root test for the SDPD model should be investigated.
Recently, another approach to describe strong spatial dependence has been proposed by Müller and Watson [49]. Since the spatial units are not neatly arranged, i.e., irregular lattice, they do not model the spatial dependence by SAR model but “posit a continuous parameter model of spatial variation [49]”. They use the Lévy–Brownian motion to define the spatial I ( 1 ) process, L ( s ) , s R d , which is a generalization of the Wiener process that is widely used in time series in d dimensions (when d = 2 , it could be regarded as a random walk on the plane). The advantage of the Lévy–Brownian motion is that such a process is isotropic, which means the relative variance between two locations is determined by the distance but not the orientation (see Anselin [1] (p. 42)). The functional central limit theorem (FCLT) is established to measure the asymptotic behavior of such a process. When regressing two independent I ( 1 ) processes, spurious regression also occurs since classical, HAC-corrected and clustered standard errors F statistics diverge to infinity, which is similar to that in Fingleton [14]. To remedy this situation, “difference” regression is again considered. Unlike time series, Müller and Watson [49] introduce the isotropic differences that treat all directions symmetrically. That is, they regard the weighted values of the neighborhood as the “average” value of the current location, just as W n Y n in the SAR model. And the difference transformation is defined as y l * = 1 n l κ b λ n 1 s s l y y l , where κ b is some weighting function. Their simulations show regressions using isotropic differences do not suffer from spurious regression problems and valid inference can be conducted. Some test procedures for H 0 : I ( 1 ) and H 0 : I ( 0 ) are also suggested.

7. Conclusions

This paper briefly surveys spatial unit roots in spatial models. First, some fundamental concepts in spatial econometrics are introduced. Spatial unit roots in SAR and SEM models may lead to spurious regression. For the estimation and inference in the presence of spatial unit roots, QMLE and 2SLS methods are generally used and have satisfactory properties after scaling. The compactness assumption has been recently relaxed in spatial econometrics which potentially makes the spatial unit roots no longer a concern but its implication to concepts like stationarity, and spatial cointegration should be further investigated. The doubly geometric spatial autoregressive process, has been widely used in some scientific fields of which the most concern is about regular lattice. Similar to time series, exact orders of convergence for different estimators are obtained because of its simple specification. But this limits its application in economics where irregular lattice and different types of weight matrices are applied.
To detect possible spatial unit roots, as well as spatial cointegration, several test procedures have been proposed. Their applications are rather limited and depend heavily on simulations to obtain critical values. This could be explained by the fact that statistics under spatial scenarios generally do not have standard asymptotic distributions, not to mention the irregular lattice.
Lastly, some related topics were introduced. The idea of singular points is generalized in SDPD model because such a model includes the time lagged variable that is based on the traditional SAR model. However, the existing literature focuses on the temporal unit roots in the SDPD model. Recently, an innovative way to study spatial unit roots describes the underlying spatial process using Lévy–Brownian motion, which is a generalization and spatial analogy to the time series counterpart. The limitations of different approaches and further research were also discussed.

Author Contributions

Conceptualization, B.H.B. and J.S.; methodology, B.H.B. and J.S.; formal analysis, B.H.B. and J.S.; writing—original draft preparation, B.H.B. and J.S.; writing—review and editing, B.H.B. and J.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments

We would like to thank the editors and four anonymous referees for their valuable comments and suggestions.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Anselin, L. Spatial Econometrics: Methods and Models; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1988. [Google Scholar]
  2. Baltagi, B.H. Econometric Analysis of Panel Data, 6th ed.; Springer: Cham, Switzerland, 2021. [Google Scholar]
  3. Elhorst, J.P. Spatial Econometrics: From Cross-Sectional Data to Spatial Panels; Springer: Heidelberg, Germany, 2014; Volume 479. [Google Scholar]
  4. Ord, K. Estimation methods for models of spatial interaction. J. Am. Stat. Assoc. 1975, 70, 120–126. [Google Scholar] [CrossRef]
  5. Lee, L.F. Asymptotic distribution of quasi-maximum likelihood estimators for spatial autoregressive models. Econometrica 2004, 72, 1899–1925. [Google Scholar] [CrossRef]
  6. Kelejian, H.H.; Prucha, I.R. A generalized moments estimator for the autoregressive parameter in a spatial model. Int. Econ. Rev. 1999, 40, 509–533. [Google Scholar] [CrossRef]
  7. Kelejian, H.H.; Prucha, I.R. A generalized spatial two-stage least squares procedure for estimating a spatial autoregressive model with autoregressive disturbances. J. Real Estate Financ. Econ. 1998, 17, 99–121. [Google Scholar] [CrossRef]
  8. Lee, L.F. Best spatial two-stage least squares estimators for a spatial autoregressive model with autoregressive disturbances. Econom. Rev. 2003, 22, 307–335. [Google Scholar] [CrossRef]
  9. Lee, L.F. GMM and 2SLS estimation of mixed regressive, spatial autoregressive models. J. Econom. 2007, 137, 489–514. [Google Scholar] [CrossRef]
  10. Yu, J.; De Jong, R.; Lee, L.F. Quasi-maximum likelihood estimators for spatial dynamic panel data with fixed effects when both n and T are large. J. Econom. 2008, 146, 118–134. [Google Scholar] [CrossRef]
  11. Baltagi, B.H.; Liu, L. Instrumental variable estimation of a spatial autoregressive panel model with random effects. Econ. Lett. 2011, 111, 135–137. [Google Scholar] [CrossRef]
  12. Kapoor, M.; Kelejian, H.H.; Prucha, I.R. Panel data models with spatially correlated error components. J. Econom. 2007, 140, 97–130. [Google Scholar] [CrossRef]
  13. Keller, W.; Shiue, C.H. The origin of spatial interaction. J. Econom. 2007, 140, 304–332. [Google Scholar] [CrossRef]
  14. Fingleton, B. Spurious spatial regression: Some Monte Carlo results with a spatial unit root and spatial cointegration. J. Reg. Sci. 1999, 39, 1–19. [Google Scholar] [CrossRef]
  15. Lee, L.F.; Yu, J. Near unit root in the spatial autoregressive model. Spat. Econ. Anal. 2013, 8, 314–351. [Google Scholar] [CrossRef]
  16. Yu, J.; Lee, L.F. Estimation of unit root spatial dynamic panel data models. Econom. Theory 2010, 26, 1332–1362. [Google Scholar] [CrossRef]
  17. Lauridsen, J.; Kosfeld, R. A test strategy for spurious spatial regression, spatial nonstationarity, and spatial cointegration. Pap. Reg. Sci. 2006, 85, 363–377. [Google Scholar] [CrossRef]
  18. Lauridsen, J.; Kosfeld, R. Spatial cointegration and heteroscedasticity. J. Geogr. Syst. 2007, 9, 253–265. [Google Scholar] [CrossRef]
  19. Beenstock, M.; Feldman, D.; Felsenstein, D. Testing for unit roots and cointegration in spatial cross-section data. Spat. Econ. Anal. 2012, 7, 203–222. [Google Scholar] [CrossRef]
  20. Yesilyurt, F.; Elhorst, J.P. A regional analysis of inflation dynamics in Turkey. Ann. Reg. Sci. 2014, 52, 1–17. [Google Scholar] [CrossRef]
  21. Olejnik, A. Using the spatial autoregressively distributed lag model in assessing the regional convergence of per-capita income in the EU25. Pap. Reg. Sci. 2008, 87, 371–385. [Google Scholar] [CrossRef]
  22. Machado, C.A.S.; Giannotti, M.A.; Chiaravalloti Neto, F.; Tripodi, A.; Persia, L.; Quintanilha, J.A. Characterization of black spot zones for vulnerable road users in São Paulo (Brazil) and Rome (Italy). ISPRS Int. J. Geo-Inf. 2015, 4, 858–882. [Google Scholar] [CrossRef]
  23. Beenstock, M.; Felsenstein, D. Estimating spatial spillover in housing construction with nonstationary panel data. J. Hous. Econ. 2015, 28, 42–58. [Google Scholar] [CrossRef]
  24. Baltagi, B.H.; Bresson, G.; Pirotte, A. Panel unit root tests and spatial dependence. J. Appl. Econom. 2007, 22, 339–360. [Google Scholar] [CrossRef]
  25. Horn, R.A.; Johnson, C.R. Matrix Analysis, 2nd ed.; Cambridge University Press: New York, NY, USA, 2012. [Google Scholar]
  26. Kelejian, H.H.; Robinson, D.P. Spatial correlation: A suggested alternative to the autoregressive model. In New Directions in Spatial Econometrics; Anselin, L., Florax, R.J.G.M., Eds.; Springer: Berlin/Heidelberg, Germany, 1995; pp. 75–95. [Google Scholar]
  27. Griffith, D.A. Advanced Spatial Statistics: Special Topics in the Exploration of Quantitative Spatial Data Series; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012; Volume 12. [Google Scholar]
  28. Hamilton, J.D. Time Series Analysis; Princeton University Press: Princeton, NJ, USA, 1994. [Google Scholar]
  29. Kelejian, H.H.; Prucha, I.R. On the asymptotic distribution of the Moran I test statistic with applications. J. Econom. 2001, 104, 219–257. [Google Scholar] [CrossRef]
  30. Kelejian, H.H.; Prucha, I.R. Specification and estimation of spatial autoregressive models with autoregressive and heteroskedastic disturbances. J. Econom. 2010, 157, 53–67. [Google Scholar] [CrossRef] [PubMed]
  31. Lee, L.F.; Yang, C.; Yu, J. QML and efficient GMM estimation of spatial autoregressive models with dominant (popular) units. J. Bus. Econ. Stat. 2023, 41, 550–562. [Google Scholar] [CrossRef]
  32. Mur, J.; Trívez, F.J. Unit roots and deterministic trends in spatial econometric models. Int. Reg. Sci. Rev. 2003, 26, 289–312. [Google Scholar] [CrossRef]
  33. Robert, E.F. Co-integration and error correction: Representation, estimation, and testing. Econometrica 1987, 55, 251–276. [Google Scholar]
  34. Lee, L.F.; Yu, J. Spatial nonstationarity and spurious regression: The case with a row-normalized spatial weights matrix. Spat. Econ. Anal. 2009, 4, 301–327. [Google Scholar] [CrossRef]
  35. Lee, L.F. Consistency and efficiency of least squares estimation for mixed regressive, spatial autoregressive models. Econom. Theory 2002, 18, 252–277. [Google Scholar] [CrossRef]
  36. Yu, J.; de Jong, R.; Lee, L.F. Estimation for spatial dynamic panel data with fixed effects: The case of spatial cointegration. J. Econom. 2012, 167, 16–37. [Google Scholar] [CrossRef]
  37. Baltagi, B.H.; Kao, C.; Liu, L. The estimation and testing of a linear regression with near unit root in the spatial autoregressive error term. Spat. Econ. Anal. 2013, 8, 241–270. [Google Scholar] [CrossRef]
  38. Baltagi, B.H.; Liu, L. Spurious spatial regression with equal weights. Stat. Probab. Lett. 2010, 80, 1640–1642. [Google Scholar] [CrossRef]
  39. Kelejian, H.H.; Prucha, I.R. 2SLS and OLS in a spatial autoregressive model with equal spatial weights. Reg. Sci. Urban Econ. 2002, 32, 691–707. [Google Scholar] [CrossRef]
  40. Krämer, W.; Donninger, C. Spatial autocorrelation among errors and the relative efficiency of OLS in the linear regression model. J. Am. Stat. Assoc. 1987, 82, 577–579. [Google Scholar]
  41. Tilke, C. The relative efficiency of OLS in the linear regression model with spatially autocorrelated errors. Stat. Pap. 1993, 34, 263–270. [Google Scholar] [CrossRef]
  42. Krämer, W.; Baltagi, B. A general condition for an optimal limiting efficiency of OLS in the general linear regression model. Econ. Lett. 1996, 50, 13–17. [Google Scholar] [CrossRef]
  43. Martellosio, F. Efficiency of the OLS estimator in the vicinity of a spatial unit root. Stat. Probab. Lett. 2011, 81, 1285–1291. [Google Scholar] [CrossRef]
  44. Baran, S.; Pap, G.; van Zuijlen, M.C. Asymptotic inference for unit roots in spatial triangular autoregression. Acta Appl. Math. 2007, 96, 17–42. [Google Scholar] [CrossRef]
  45. Basu, S.; Reinsel, G.C. A note on properties of spatial Yule-Walker estimators. J. Stat. Comput. Simul. 1992, 41, 243–255. [Google Scholar] [CrossRef]
  46. Basu, S.; Reinsel, G.C. Properties of the spatial unilateral first-order ARMA model. Adv. Appl. Probab. 1993, 25, 631–648. [Google Scholar] [CrossRef]
  47. Baran, S. On the variances of a spatial unit root model. Lith. Math. J. 2011, 51, 122–140. [Google Scholar] [CrossRef]
  48. Paulauskas, V. On unit roots for spatial autoregressive models. J. Multivar. Anal. 2007, 98, 209–226. [Google Scholar] [CrossRef]
  49. Müller, U.K.; Watson, M.W. Spatial Unit Roots; Princeton University: Princeton, NJ, USA, 2022. [Google Scholar]
  50. Gupta, A. Estimation of spatial autoregressions with stochastic weight matrices. Econom. Theory 2019, 35, 417–463. [Google Scholar] [CrossRef]
  51. Liu, T.; Xu, X.; Lee, L.F. Consistency without compactness of the parameter space in spatial econometrics. Econ. Lett. 2022, 210, 110224. [Google Scholar] [CrossRef]
  52. Newey, W.K.; McFadden, D. Large sample estimation and hypothesis testing. In Handbook of Econometrics; Elsevier: Amsterdam, The Netherlands, 1994; Volume 4, Chapter 35; pp. 2111–2245. [Google Scholar]
  53. Gupta, A. Efficient closed-form estimation of large spatial autoregressions. J. Econom. 2023, 232, 148–167. [Google Scholar] [CrossRef]
  54. Rossi, F.; Lieberman, O. Spatial autoregressions with an extended parameter space and similarity-based weights. J. Econom. 2023, 235, 1770–1798. [Google Scholar] [CrossRef]
  55. Liu, L. A note on 2SLS estimation of the mixed regressive spatial autoregressive model. Econ. Lett. 2015, 134, 49–52. [Google Scholar] [CrossRef]
  56. Bhattacharyya, B.; Richardson, G.; Franklin, L. Asymptotic inference for near unit roots in spatial autoregression. Ann. Stat. 1997, 25, 1709–1724. [Google Scholar] [CrossRef]
  57. Bhattacharyya, B.; Khalil, T.; Richardson, G. Gauss-Newton estimation of parameters for a spatial autoregression model. Stat. Probab. Lett. 1996, 28, 173–179. [Google Scholar] [CrossRef]
  58. Phillips, P.C. Time series regression with a unit root. Econometrica 1987, 55, 277–301. [Google Scholar] [CrossRef]
  59. Baran, S.; Pap, G. Parameter estimation in a spatial unilateral unit root autoregressive model. J. Multivar. Anal. 2012, 107, 282–305. [Google Scholar] [CrossRef]
  60. Baran, S.; Pap, G.; Van Zuijlen, M.C. Asymptotic inference for a nearly unstable sequence of stationary spatial AR models. Stat. Probab. Lett. 2004, 69, 53–61. [Google Scholar] [CrossRef]
  61. Baran, S.; Pap, G. On the least squares estimator in a nearly unstable sequence of stationary spatial AR models. J. Multivar. Anal. 2009, 100, 686–698. [Google Scholar] [CrossRef]
  62. Roknossadati, S.; Zarepour, M. M-estimation for a spatial unilateral autoregressive model with infinite variance innovations. Econom. Theory 2010, 26, 1663–1682. [Google Scholar] [CrossRef]
  63. Roknossadati, S.; Zarepour, M. M-estimation for near unit roots in spatial autoregression with infinite variance. Statistics 2011, 45, 337–348. [Google Scholar] [CrossRef]
  64. Ahlgren, N.; Gerkman, L. Inference in unilateral spatial econometric models. Bull. Int. Stat. Inst. 2007, 56, 1–44. [Google Scholar]
  65. Lauridsen, J.; Kosfeld, R. A Wald test for spatial nonstationarity. Estud. Econ. Apl. 2004, 22, 475–486. [Google Scholar]
  66. Martellosio, F. Power properties of invariant tests for spatial autocorrelation in linear regression. Econom. Theory 2010, 26, 152–186. [Google Scholar] [CrossRef]
  67. King, M.L. A small sample property of the Cliff-Ord test for spatial correlation. J. R. Stat. Soc. Ser. B (Methodological) 1981, 43, 263–264. [Google Scholar] [CrossRef]
  68. Krämer, W. Finite sample power of Cliff–Ord-type tests for spatial disturbance correlation in linear regression. J. Stat. Plan. Inference 2005, 128, 489–496. [Google Scholar] [CrossRef]
  69. Martellosio, F. Testing for spatial autocorrelation: The regressors that make the power disappear. Econom. Rev. 2012, 31, 215–240. [Google Scholar] [CrossRef]
  70. Born, B.; Breitung, J. Simple regression-based tests for spatial dependence. Econom. J. 2011, 14, 330–342. [Google Scholar] [CrossRef]
  71. Baltagi, B.H.; Yang, Z. Heteroskedasticity and non-normality robust LM tests for spatial dependence. Reg. Sci. Urban Econ. 2013, 43, 725–739. [Google Scholar] [CrossRef]
  72. Baltagi, B.H.; Yang, Z. Standardized LM tests for spatial error dependence in linear or panel regressions. Econom. J. 2013, 16, 103–134. [Google Scholar] [CrossRef]
  73. Preinerstorfer, D. How to avoid the zero-power trap in testing for correlation. Econom. Theory 2021, 39, 1292–1324. [Google Scholar] [CrossRef]
  74. Anselin, L. Lagrange multiplier test diagnostics for spatial dependence and spatial heterogeneity. Geogr. Anal. 1988, 20, 1–17. [Google Scholar] [CrossRef]
  75. Baltagi, B.H.; Song, S.H.; Koh, W. Testing panel data regression models with spatial error correlation. J. Econom. 2003, 117, 123–150. [Google Scholar] [CrossRef]
  76. Beenstock, M.; Felsenstein, D. Unit root and cointegration tests in spatial cross-section data. In The Econometric Analysis of Non-Stationary Spatial Panel Data; Advances in Spatial Science; Springer: Berlin/Heidelberg, Germany, 2019; Chapter 5; pp. 97–127. [Google Scholar]
  77. Kosfeld, R.; Lauridsen, J. Dynamic spatial modelling of regional convergence processes. Empir. Econ. 2004, 29, 705–722. [Google Scholar] [CrossRef]
  78. Vaona, A. Spatial autocorrelation and the sensitivity of RESET: A simulation study. J. Geogr. Syst. 2010, 12, 89–103. [Google Scholar] [CrossRef]
  79. Fingleton, B. A generalized method of moments estimator for a spatial model with moving average errors, with application to real estate prices. Empir. Econ. 2008, 34, 35–57. [Google Scholar] [CrossRef]
  80. Baltagi, B.H.; Liu, L. Testing for random effects and spatial lag dependence in panel data models. Stat. Probab. Lett. 2008, 78, 3304–3306. [Google Scholar] [CrossRef]
  81. Lee, L.F.; Yu, J. A unified transformation approach for the estimation of spatial dynamic panel data models: Stability, spatial cointegration and explosive roots. In Handbook on Empirical Economics and Finance; Ullah, A., Giles, D., Eds.; Taylor & Francis Group: New York, NY, USA, 2011; Chapter 13; pp. 397–434. [Google Scholar]
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Baltagi, B.H.; Shu, J. A Survey of Spatial Unit Roots. Mathematics 2024, 12, 1052. https://doi.org/10.3390/math12071052

AMA Style

Baltagi BH, Shu J. A Survey of Spatial Unit Roots. Mathematics. 2024; 12(7):1052. https://doi.org/10.3390/math12071052

Chicago/Turabian Style

Baltagi, Badi H., and Junjie Shu. 2024. "A Survey of Spatial Unit Roots" Mathematics 12, no. 7: 1052. https://doi.org/10.3390/math12071052

APA Style

Baltagi, B. H., & Shu, J. (2024). A Survey of Spatial Unit Roots. Mathematics, 12(7), 1052. https://doi.org/10.3390/math12071052

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop